# Using COSMO-SAC models in the mixing rules of the translated-consistent Peng-Robinson equation-of-state: a rigorous benchmarking across a broad and diverse database

Francisco Paes, Baptiste Sirjean, Jean-Noël Jaubert, Romain Privat

► To cite this version:

Francisco Paes, Baptiste Sirjean, Jean-Noël Jaubert, Romain Privat. Using COSMO-SAC models in the mixing rules of the translated-consistent Peng-Robinson equation-of-state: a rigorous benchmarking across a broad and diverse database. Chemical Engineering Science, 2026, 324, pp.123300. 10.1016/j.ces.2026.123300. hal-05458297

HAL Id: hal-05458297

https://hal.univ-lorraine.fr/hal-05458297v1

Submitted on 14 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Using COSMO-SAC models in the mixing rules of the translated-consistent Peng-Robinson equation-of-state: a rigorous benchmarking across a broad and diverse database

Francisco Paes ID, Baptiste Sirjean *, Jean-Noël Jaubert * ID, Romain Privat * ID

Université de Lorraine, CNRS, LRGP, F-54000 Nancy, France

## ARTICLE INFO

**Keywords:**

COSMO-SAC
Equation of State
Mixing Rules
Phase equilibria

## ABSTRACT

Among the various continuum solvation models available, the COSMO-type ones are the most accurate option as they incorporate a robust statistical thermodynamic formulation combined with quantum-based calculations, ensuring a good level of accuracy and a high predictive power. Therefore, COSMO-type models offer an interesting alternative to group contribution-based approaches like UNIFAC, whose applicability relies on the availability of functional groups and their associated interaction parameters. When it comes to calculating thermodynamic data at low and high pressures, UNIFAC is typically used as the excess Gibbs energy model within the mixing rules of cubic equations of state, such as PSRK, VTPR, and UMR-PRU. In view of the advantages associated with the COSMO-based approach, this work aims to systematically analyze the use of COSMO-based models within mixing rules of cubic equations of state. For this purpose, three versions of the COSMO-SAC model (COSMO-SAC 2002, 2010, and dsp) were implemented (and re-parameterized) alongside the translated-consistent Peng-Robinson equation of state (tc-PR EoS) for a comprehensive benchmarking study based on a large and diverse experimental database of 160 binary mixtures. This database comprises data on vapor-liquid equilibria at both low and high pressures, as well as critical points and mixing properties, such as enthalpy and heat capacity. By combining the tc-PR EoS with more advanced COSMO-SAC models (2010 and dsp), it was possible to effectively represent the phase behavior of a wide range of mixtures, consistently reproducing complex phase diagrams at least from a qualitative point of view and, in many cases, also quantitatively.

## 1. Introduction

Predictive thermodynamic models are essential for estimating the behavior of mixtures without relying on extensive experimental data. This capability is key in the context of sustainable development, where designing new products and processes plays a major role (De Hemptinne et al., 2022). These models are useful for predicting phase equilibria and the thermodynamic properties of pure compounds and mixtures. They are important tools for process screening, design, and optimization, especially in the early stages of industrial projects.

Kontogeorgis et al. recently conducted a survey on industrial demand for thermodynamic modeling (Kontogeorgis et al., 2021). Regarding predictive models for mixtures, the study revealed a growing interest in COSMO-type models over the past decade, mainly due to their distinct advantages over traditional predictive models based on group

contributions, which include:

* No need for an extensive experimental database to fit Binary Interaction Parameters (BIPs) between functional groups, since COSMO-based models only depend on pure compound information obtained through quantum-chemistry calculations.

* It is not limited to the availability of predefined functional groups, increasing its scope of application to a wider range of chemical species.

* It can handle different isomeric forms, conformations, and structural arrangements, capturing the effect of geometric variations on the thermodynamic behavior of molecules.

COSMO-type models are notably versatile in handling diverse chemical species; the only prerequisite is the knowledge of their atomic

<footnote>* Corresponding authors.</footnote>

<footnote>E-mail addresses: francisco.paes@univ-lorraine.fr (F. Paes), baptiste.sirjean@univ-lorraine.fr (B. Sirjean), jean-noel.jaubert@univ-lorraine.fr (J.-N. Jaubert), romain.privat@univ-lorraine.fr (R. Privat).</footnote>
coordinates, which are used as input data to perform quantum-based COSMO calculations in molecular cavities to predict liquid phase effects (Klamt, 1995; Mullins et al., 2006). Despite their flexibility, challenges persist. For instance, it is well-known that COSMO-based models are originally adapted to describe incompressible liquids far from their critical region.

Many research studies have focused on extending the application of COSMO-based models to high temperatures and high pressures. For example, Costa et al. (Costa et al., 2016a, 2016b) developed the σ-MTC equation of state, which integrates the charge density profiles from COSMO computations with the generalized van der Waals theory (Vera and Prausnitz, 1972). More recently, Soares et al. (de P. Soares et al., 2019) proposed a novel COSMO-based equation of state (EoS), referred to as COSMO-SAC-phi, which couples the COSMO-SAC model (Lin and Sandler, 2002) with the lattice-fluid theory by incorporating free volume in the form of holes.

It has also been thoroughly demonstrated that COSMO-type models can be successfully combined with cubic EoSs through the so-called advanced mixing rules. One of the first investigations to explore this combination was conducted by Constantinescu et al., in 2005 (Constantinescu et al., 2005). The authors combined the COSMO-RS (Eckert and Klamt, 2002; Klamt, 1995) and COSMO-SAC 2002 (Lin and Sandler, 2002) models with the Soave-Redlich-Kwong (SRK) EoS using the Huron-Vidal mixing rule (Huron and Vidal, 1979). Encouraging results were obtained, with average relative deviations of less than 5 % in their predictions of the saturation pressure predictions of six binary mixtures. A year later, Shimoyama et al. (Shimoyama et al., 2006) performed a similar study coupling the SRK EoS with UNIQUAC, whose BIP's were fitted to generic data generated with COSMO-RS. This study focused on vapor–liquid and liquid–liquid equilibria of a few binary and ternary mixtures involving alcohols and hydrocarbons. Similarly, Lee and Lin (Lee and Lin, 2007) investigated the coupling of COSMO-SAC with the Peng-Robinson EoS through the Wong-Sandler mixing rule, analyzing 12 binary systems containing alkanes, alcohols, water, and acetone. In that study, the authors found that using only the residual part of the COSMO-SAC model in the mixing rule yielded superior results, especially for asymmetric mixtures at high temperatures.

Leonhard et al. (Leonhard et al., 2009) conducted a more detailed study to evaluate the performance of COSMO-RS in combination with the Peng-Robinson EoS using three different mixing rules: Wong-Sandler, MHV1 (first order Modified-Huron-Vidal) and Ahlers-Gmehling mixing rules. Their analysis was conducted on liquid-vapor equilibrium data from 31 binary systems. The authors concluded that incorporating a COSMO-type model into the mixing rule of a cubic EoS does not improve the representation of a system if the COSMO-type model alone (i.e., when used in a $γ-ω$ approach) is already inaccurate. In other words, the authors confirmed that incorporating COSMO-RS within a mixing rule extends the method to high-pressure conditions rather than improving its accuracy. A similar study was carried out by Hsieh and Lin (Hsieh and Lin, 2012), but with emphasis on pressure effects on liquid-liquid equilibrium (LLE). The authors showed that adopting a mixing rule with an infinite reference pressure improves the representation of the pressure dependence of LLE. On the other hand, using a mixing rule with zero reference pressure produces results nearly identical to those of the standalone COSMO-SAC model, reducing the influence of using a cubic equation of state.

Building on the work of Mambo-Lomba and Paricaud (Mambo-Lomba and Paricaud, 2021), Wang et al. (Wang et al., 2023) used an advanced version of the COSMO-SAC model – COSMO-SAC-dsp (Hsieh et al., 2014) – to predict the phase equilibria of 233 binary systems involving refrigerant fluids, mostly halogenated hydrocarbon mixtures). For this type of system, the authors found that the Peng-Robinson EoS based on COSMO-SAC-dsp outperformed the PSRK EoS. However, this improvement was mainly in predictive power, since the UNIFAC method applied to PSRK lacks several binary interaction parameters involving halogenated groups.

Recently, Prino and Voutsas (Prinos and Voutsas, 2025) introduced a revised version of the COSMO-SAC 2010 model, featuring a more detailed description of screening charge density profiles through a specific treatment for water molecules, along with a more consistent description of the combinatorial term. The authors tested the modified COSMO-SAC 2010 model at high pressures using the Peng-Robinson EoS, relying on the universal mixing rule previously proposed for the UMR-PRU EoS (Voutsas et al., 2006). Compared to the original COSMO-SAC 2010 version, the new model improved the description of azeotropic mixtures due to the more detailed description of hydrogen bonds (which implies a greater number of fitted parameters). Additionally, size-asymmetric systems, especially those involving gases, were better represented due to the new formulation of the combinatorial term. To support their conclusions, the authors evaluated high-pressure vapor-liquid equilibrium (VLE) data for 12 binary systems in total.

While many studies have examined the combination of cubic equations of state with the COSMO-based models, a more comprehensive analysis would contribute meaningfully to the literature. Such an analysis should involve comparisons across an extensive and diverse database encompassing a broad range of conditions and species. To fill this gap, this study presents such an analysis by benchmarking the combination of a cubic EoS with COSMO-type models across 160 binary systems comprising low and high pressure data, using a database developed to evaluate the performance of EoSs (Jaubert et al., 2020).

For this purpose, the translated-consistent Peng-Robinson (tc-PR (Le Guennec et al., 2016b)) EoS was combined with COSMO-SAC models in versions 2002 (Lin and Sandler, 2002), 2010 (Hsieh et al., 2010), and 2014 (Hsieh et al., 2014) – also known as COSMO-SAC-dsp – through the recently developed Unified Huron Vidal Michelsen (UHVM) mixing rule (Privat et al., 2025). The tc-PR EoS was selected because it provides reliable predictions for pure-component properties, combining a thermodynamically consistent $\alpha$-function with an appropriate volume-translation constant. Regarding the UHVM mixing rule, it offers the advantage of being physically sound, as it incorporates only the residual contribution to the excess Gibbs energy. This approach therefore avoids inconsistencies between the combinatorial term derived from the equation of state and that of the COSMO-SAC model.

Each COSMO-SAC variant was re-parameterized (i.e., its universal constants were refitted) to fully exploit the tc-PR/COSMO-SAC EoS. These versions were selected to represent the most widely used COSMO-SAC formulations: the 2002 model as a baseline, the 2010 model to illustrate the effect of a more detailed $\sigma$-profile description, and the dp model to highlight the impact of including a specific dispersion term, in contrast to the fully electrostatic residual contribution obtained from the other versions.

Finally, to support the evaluation of the results, we applied the scoring system proposed by Jaubert et al. (Jaubert et al., 2020), which enables a quick and effective comparison of the performance of each tc-PR/COSMO-SAC variant (and with other models).

## 2. Modeling

This section presents the formulation of the tc-PR cubic EoS combined with the “Unified Huron-Vidal Michelsen” (UHVM) mixing rule (Privat et al., 2025), together with the COSMO-SAC model variants considered in this work (COSMO-SAC 2002, COSMO-SAC 2010, and COSMO-SAC-dsp). Together, these formulations form the basis for the phase-equilibria and mixing property calculations discussed in Section 3.

### 2.1. The tc-PR equation of state with advanced mixing rules

#### 2.1.1. EoS formulation

For a given species i, the formulation of the tc-PR EoS is introduced in Eq. (1) below.
$$ P = \frac{RT}{(V_m + c_i - b_i)} - \frac{a_i(T)}{(V_m + c_i)(V_m + c_i + b_i) + b_i(V_m + c_i - b_i)} $$ 
 where: 
 $$ \left\{ \begin{array}{l} a_i(T) = 0.45724 \frac{(RT_{c,exp,i})^2}{P_{c,exp,i}} \alpha_i^{Twu-91}(T_{r,i}) \\ \alpha_i^{Twu-91}(T_{r,i}) = T_{r,i}^{N(M_i-1)} \exp[L_i(1-T_{r,i}^{M_iN_i})] \\ \text{with: } T_{r,i} = \frac{T}{T_{c,exp,i}} \\ b_i = 0.0778 \frac{RT_{c,exp,i}}{P_{c,exp,i}} \\ c_i = V_{m,non\ translated\ EoS}(T_{r,i} = 0.8) - V_{m,exp}(T_{r,i} = 0.8) \end{array} \right. \quad (1) $$

Where T is the temperature, P is the pressure, and Vm the molar volume. The tc-PR EoS has three parameters: the attractive parameter (aᵢ), the co-volume of the non-translated EoS (bᵢ) and a constant correction called volume translation (cᵢ). The parameter aᵢ characterizes the attractive forces between molecules, the parameter bᵢ is the volume occupied by one mole of the substance at infinite pressure (when the EoS is not translated), and cᵢ is the volume translation parameter, which has the role of improving the description of volumetric properties (Jaubert et al., 2016; Peneloux et al., 1989).

As shown in Eq. (1), both aᵢ and bᵢ are derived using the experimental values of the critical temperature (T<sub>c,exp,i</sub>) and critical pressure (P<sub>c,exp,i</sub>) of the substance i. The temperature dependence of aᵢ is calculated using the Twu-91 α-function, which ensures thermodynamic consistency across a broad temperature range (Le Guennec et al., 2016a). This α-function employs three parameters (Lᵢ, Mᵢ, and Nᵢ), either fitted to experimental data or determined through generalized correlations. The volume translation constant (cᵢ) is determined to exactly reproduce the experimental value of the saturated-liquid molar volume of the pure component i at a reduced temperature T<sub>r,i</sub> = 0.8. Optimal values of Lᵢ, Mᵢ, Nᵢ and cᵢ for 1,800 compounds can be found elsewhere (Le Guennec et al., 2016b; Magoulas and Tassios, 1990; Pina-Martinez et al., 2022).

### 2.1.2. Mixing rules

Advanced mixing rules are expressions for the mixture parameters a<sub>mix</sub> and b<sub>mix</sub> obtained by equating the expression of a given excess Gibbs energy model to the expression of the excess Gibbs energy stemming from the EoS (g<sup>E,EoS</sup> = g<sup>E,y</sup>). Deriving these expressions for a<sub>mix</sub> and b<sub>mix</sub> in this way provides a more accurate description of complex polar and associating mixtures than the classical Van der Waals mixing rules do. Following a recent discussion on the strengths and weaknesses of advanced mixing rules for cubic EoS (Privat et al., 2025), using the equality of residual excess Helmholtz energies (a<sub>res</sub><sup>E,EoS</sup> = a<sub>res</sub><sup>E,y</sup>) rather than g<sup>E,EoS</sup> = g<sup>E,y</sup> as fundamental equation (sometimes named matching equation) to derive the expressions of a<sub>mix</sub> and b<sub>mix</sub> eliminates most of the limitations of this approach and, particularly the potential inconsistencies introduced by the combinatorial part of activity coefficient models (for more details, the reader is referred to (Privat et al., 2025)). The result of this identification is the UHVM (Unified Huron Vidal Michelsen) mixing rule, and is given by Eq. (2).

$$ \frac{a_{\text{mix}}(T, \mathbf{x})}{b_{\text{mix}}(\mathbf{x})} = \frac{a_{\text{res}}^{E,y}(T, \mathbf{x})}{\Lambda_{\text{EoS}}} + \sum_{i} x_i \frac{a_i(T)}{b_i} $$ 
 with 
 $$ \Lambda_{\text{EoS}} = -0.53 \quad (2) $$

Note that applying the matching equation a<sub>res</sub><sup>E,EoS</sup> = a<sub>res</sub><sup>E,y</sup> results in the same formulation of the mixing rule for a<sub>mix</sub>/b<sub>mix</sub>, regardless of the

reference pressure considered (0 or +∞). Despite the similarity with the mixing rule proposed by Huron and Vidal, this formulation allows for a quadratic mixing rule for b<sub>mix</sub> (as shown in Eq. (3)), rather than the linear combination of pure-component co-volumes required in the Huron-Vidal mixing rule, which often yields poorer results.

$$ \begin{array}{l} b_{\text{mix}}(\mathbf{x}) = \sum_{i} \sum_{j} x_i x_j b_{ij} \\ \text{with: } \left\{ \begin{array}{l} b_{ij} = \left( \frac{b_i^{1/s} + b_j^{1/s}}{2} \right)^s \\ s = 3/2 \end{array} \right. \end{array} \quad (3) $$

In agreement with theory (Privat et al., 2016), the c<sub>mix</sub> coefficient is often assigned a linear mixing rule.

$$ c_{\text{mix}}(\mathbf{x}) = \sum_{i} x_i c_i \quad (4) $$

In this study, the computation of the activity coefficients (γ<sub>i</sub><sup>res</sup>) that appear in Eq. (2) was performed using three different versions of the COSMO-SAC model (COSMO-SAC 2002, COSMO-SAC 2010, and COSMO-SAC-dsp), which are briefly presented below.

## 2.2. Calculating activity coefficients with COSMO-SAC models

In the COSMO-SAC groundwork, the logarithm of the activity coefficient lnγ<sub>i,S</sub> of a solute i in a mixture S(the solvent) is classically calculated by summing the residual (lnγ<sub>i,S</sub><sup>res</sup>) and the combinatorial (lnγ<sub>i,S</sub><sup>res</sup>) terms (in other words, the excess chemical potential of i results from attractive – residual – contributions and repulsive – combinatorial – contributions). The residual term can itself be divided into a dispersion contribution (lnγ<sub>i,S</sub><sup>res,disp</sup>) and an electrostatic contribution (lnγ<sub>i,S</sub><sup>res,elec</sup>).

$$ \ln\gamma_{i,S} = \frac{\ln(\gamma_{i,S}^{\text{res,elec}}) + \ln(\gamma_{i,S}^{\text{res,disp}}) + \ln(\gamma_{i,S}^{\text{comb}})}{\ln(\gamma_{i,S}^{\text{res}})} \quad (5) $$

The combinatorial part is usually represented by a Staverman-Guggenheim term. This accounts for the non-ideal contribution of a mixture arising from differences in the size and shape of its components. As can be seen from Eq. (2), this term is not involved in the mixing rule for the quantity a<sub>mix</sub>/b<sub>mix</sub>. Conversely, both the dispersive and electrostatic residual contributions share the same enthalpic origin, arising from molecular interactions that affect the energy of the system: the electrostatic term describes the permanent electrostatic interactions between molecules in the mixture while the dispersive term reflects non-electrostatic enthalpic effects (Hsieh et al., 2014). Therefore, both contributions belong to the family of Van der Waals forces and must be involved in the UHVM mixing rule introduced in Eq. (2). It should be noted that the COSMO-SAC 2002 and 2010 models do not include an explicit dispersive contribution. However, the COSMO-SAC-dsp model, an improved version of COSMO-SAC 2010, incorporates a dispersion term. These aspects will be examined in greater detail in the following paragraphs.

### 2.2.1. Screening charges as the descriptor of residual interactions

The COSMO-SAC models use the screening charge densities of surface segments on molecular cavities as descriptors for calculating the residual contribution to activity coefficients. These very localized charge densities, obtained from DFT-based COSMO calculations, are incorporated into a lattice-based framework analogous to Guggenheim's Quasi-Chemical Theory (as employed with functional groups in models like UNIFAC). The three-dimensional surface charge information is simplified into histograms known as σ-profiles, which serve as a key input for evaluating segment interactions and ultimately computing activity
<figure><img src="image_2.png" /><figcaption>Fig. 1. σ-profiles of acetic acid within the COSMO-SAC 2002 and COSMO-SAC 2010 frameworks.</figcaption></figure>

coefficients of species in liquid mixtures.

The 2002 version of COSMO-SAC employs a single σ-profile for representing molecular surface charge distributions. In contrast, the 2010 version introduces three distinct σ-profiles to classify the surface segments into the following categories:

* OH: (1) The segment belongs to an oxygen atom bound to a hydrogen atom and its charge density is positive, or, (2) the segment belongs to a hydrogen atom bound to an oxygen atom and its charge is negative.

* OT: (1) The segment belongs to an oxygen atom (bound to an atom other than hydrogen), a nitrogen or a fluorine atom, and its charge density is positive, or, (2) the segment belongs to a hydrogen atom bound to a nitrogen or fluorine atom, and its charge density is negative.

* nhb: otherwise (nhb = non-hydrogen bond).

Fig. 1 illustrates the differences in σ-profiles between the 2002 and 2010 versions of COSMO-SAC, using acetic acid as an example. An overview of how the σ-profiles are calculated across the different versions of COSMO-SAC is given in the Supplementary Material (refer to Appendix A for the details about the organization of the supplementary material).

### 2.2.2. Electrostatic residual contributions from COSMO-SAC 2002

As mentioned above, the electrostatic residual contribution in the COSMO-SAC 2002 model is determined using a single σ-profile, as illustrated in Eq. (6).

$$ \ln(Y_{i,S}^{\text{res.elec}}) = n_i \sum_{\sigma_m} p_i(\sigma_m) [\ln(\Gamma_S(\sigma_m)) - \ln(\Gamma_i(\sigma_m))] \quad (6) $$

Where $\Gamma_S$ and $\Gamma_i$ represent the activity coefficients of a given segment $m$ in the mixture $S$ and in pure $i$, respectively. Moreover, $n_i$ is an estimate of the number of segments of a molecule $i$ ($n_i = A_i/a_{eff}$). The activity coefficients of the segments, either in a mixture or in a pure compound, are calculated by Eq. (7). In the case of pure compound $i$, $p_S(\sigma)$ and $\Gamma_S$ are replaced by $p_i(\sigma)$ and $\Gamma_i$.

$$ \ln(\Gamma_S(\sigma_m)) = - \ln \left\{ \sum_{\sigma_n} p_S(\sigma_n) \Gamma_S(\sigma_n) \exp \left[ - \frac{\Delta W(\sigma_m, \sigma_n)}{RT} \right] \right\} \quad (7) $$

Here, the parameter $\Delta W(\sigma_m, \sigma_n)$ represents the contact energy between two segments $\sigma_m$ and $\sigma_n$, incorporating both electrostatic and hydrogen bonding effects, as shown in Eqs (8).

$$ \Delta W(\sigma_m, \sigma_n) = \underbrace{\frac{\alpha'}{2}(\sigma_m + \sigma_n)^2}_{\text{electrostatic}} + \underbrace{C_{hb}(T) \times \min\{0, \min[0, \sigma_{acc} + \sigma_{hb}]\} \max[0, \sigma_{don} - \sigma_{hb}]}_{\text{Hydrogen - bonding}} \quad (8) $$
with:
$$ \begin{cases} \sigma_{acc} = \max\{\sigma_m; \sigma_n\} \\ \sigma_{don} = \min\{\sigma_m; \sigma_n\} \end{cases} $$
<figure><img src="image_3.png" /><figcaption>Fig. 2. In (a) the temperature-dependent function of the electrostatic (or misfit) constant ($C_{ES}$) proposed by Hsieh et al. (Hsieh et al., 2010). In (b) a phase diagram for the system argon + methane at 90.68 K generated by the models tc-PR/COSMO-SAC 2002 and tc-PR/COSMO-SAC 2010. Note the discrepancy between the results obtained with tc-PR/COSMO-SAC 2010 and the experimental data due to the use of the $C_{ES}(T)$ function at low temperatures.</figcaption></figure>

Where, $\alpha'$ is the electrostatic interaction parameter, $\sigma_{hb}$ represents a cut-off value for hydrogen-bonding, and $C_{hb}$ denotes the hydrogen-bonding interaction constant. A temperature-dependent function was introduced for $C_{hb}$ to account for the reduction in entropy that goes along with hydrogen bonding (Klamt, 2005).

$$f_{hb}(T) = T\ln\left[1 + 0.005\exp\left(\frac{4.78 \text{ kcal/mol}}{RT}\right)\right] \quad (9)$$

$$C_{hb}(T) = C_{hb} \times \frac{f_{hb}(T)}{f_{hb}(298.15 \text{ K})} \quad (10)$$

This function is important for better capturing the temperature dependence of hydrogen-bonding interactions, which progressively weaken at higher temperatures.

### 2.2.3. Electrostatic residual contributions from COSMO-SAC 2010

Based on the separation of $\sigma$-profiles discussed earlier, Hsieh et al. (Hsieh et al., 2010) introduced an enhanced version of the COSMO-SAC model for calculating the electrostatic residual contribution to activity coefficients. This improved model is presented in Eq. (11).

$$\ln(\gamma_{i,S}^{\text{res,elec}}) = n_i \sum_{t}^{nhb,OH,OT} \sum_{\sigma_m} p_t^s(\sigma_m) [\ln(\Gamma_S^t(\sigma_m)) - \ln(\Gamma_i^t(\sigma_m))] \quad (11)$$

Similar to the model proposed by Lin and Sandler in 2002, $\Gamma_S^t$ and $\Gamma_i^t$ also denote the activity coefficient of a given segment $m$ in the mixture $S$ and in the pure compound $i$, but, in this case, associated to one of the t-type $\sigma$-profiles ($t = nhb, OH, or OT$). The calculation of the segment activity coefficients is performed by Eq. (12).

$$\ln(\Gamma_S^t(\sigma_m^t)) = -\ln\left\{\sum_{s}^{nhb,OH,OT} \sum_{\sigma_n} p_S^s(\sigma_n^s) \Gamma_S^s(\sigma_n^s) \exp\left[-\frac{\Delta W(\sigma_m^t, \sigma_n^s)}{RT}\right]\right\} \quad (12)$$

In this version of the COSMO-SAC model, the contact energy between the segments $\sigma_m^t$ and $\sigma_n^s$ is calculated by Eq. (13).

$$\Delta W(\sigma_m^t, \sigma_n^s) = \underbrace{\text{electrostatic}}_{\text{C}_{ES}(T) \times (\sigma_m^t + \sigma_n^s)^2} - \underbrace{\text{hydrogen-bonding}}_{C_{hb}(\sigma_m^t, \sigma_n^s) \times (\sigma_m^t - \sigma_n^s)^2} \quad (13)$$

Note that the electrostatic part does not distinguish between $\sigma$-profiles, since the same electrostatic parameter ($C_{ES}$) is used for all types of interaction. Hsieh et al. (Hsieh et al., 2010) made this parameter an empirical function of temperature, as follows ($C_{ES}$ is given in kcal Å$^4$ mol$^{-1}$ e$^{-2}$):

$$C_{ES}(T) = 6525.69 + \frac{1.4859 \times 10^8}{T^2} \quad (14)$$

However, Fig. 2(a) shows that this temperature-dependent function produces inconsistent results at low temperatures (below 250 K) due to its divergence towards infinity at 0 K. This divergence problem does not occur at high temperatures, as $C_{ES}$ asymptotically approaches a reasonable constant value Fig. 2(b) shows the Pxy projection of the phase diagram at 91 K for the system argon + methane using both the tc-PR/COSMO-SAC 2002 and tc-PR/COSMO-SAC 2010 models.

In order to avoid this divergence issue, we opted to keep the electrostatic parameter constant, as in the COSMO-SAC 2002 version:

$$C_{ES} = \frac{\alpha'}{2} \quad (15)$$

In COSMO-SAC 2010, the hydrogen bond interaction is calculated differently based on the nature of segments $t$ and $s$, as the interaction constant $C_{hb}$ assumes different values depending on the $t/s$ combination, as shown below. It should be noted that the same temperature-dependent function introduced in Eq. (10) was used in this COSMO-SAC 2010 implementation.

$$C_{hb}(\sigma_m^t, \sigma_n^s) = \begin{cases} c_{OH-OH}, & \text{if } s = t = \text{OH} \text{ and } \sigma_m^t \times \sigma_n^s < 0 \\ c_{OT-OT}, & \text{if } s = t = \text{OT} \text{ and } \sigma_m^t \times \sigma_n^s < 0 \\ c_{OH-OT}, & \text{if } s = \text{OH}, t = \text{OT} \text{ and } \sigma_m^t \times \sigma_n^s < 0 \\ 0, & \text{otherwise} \end{cases} \quad (16)$$

The parameters $c_{OH-OH}$, $c_{OH-OT}$ and $c_{OT-OT}$ are universal constants fitted to experimental data. This formulation of hydrogen bonds provides more detail than the COSMO-SAC 2002 version proposed by Lin and Sandler in 2002, since this older version – as well as the different versions of the COSMO-RS model (Gerlach et al., 2022; Klamt et al., 2010; Klamt and Eckert, 2000) – uses a single universal constant to represent all types of hydrogen bonding interactions, regardless their nature. For instance, O–H···O or N–H···N hydrogen bond interactions are calculated using the same universal constant, $C_{hb}$. This often results in less accurate phase-equilibrium predictions for systems in which hydrogen bonding plays a dominant role, especially mixtures containing self-associating molecules, such as water, alcohols, carboxylic acids or amines. The COSMO-SAC 2010 model captures variations in association strength accompanying hydrogen bonding more accurately by handling association with greater specificity.
### 2.2.4. Dispersion contributions from COSMO-SAC-dsp

In 2014, Hsieh and collaborators (Hsieh et al., 2014) proposed an enhanced version of the COSMO-SAC 2010 model by introducing a dispersion term (the so-called COSMO-SAC-dsp model), while retaining the same electrostatic residual contribution formulation (and its universal constants) as in their model from 2010. A dispersive term is added to the COSMO-SAC 2010 version through a simplified model based on the 1-parameter Margules equation. For a binary system, the dispersive term is written as follows:

$$ 
\left\{ 
\begin{aligned}
\ln(\gamma_{1,S}^{\text{res,disp}}) &= A_{\text{disp}}x_2^2 \\
\ln(\gamma_{2,S}^{\text{res,disp}}) &= A_{\text{disp}}x_1^2
\end{aligned}
\right. 
\qquad (17) 
$$

Where the parameter $A_{\text{disp}}$ is calculated based on the dispersion parameters ($\epsilon_1/k_b$) of each component of the binary system through Eq. (18) (Bell et al., 2020; Hsieh et al., 2014).

$$ A_{\text{disp}} = w \left[ \frac{1}{2} \left( \frac{\epsilon_1}{k_b} + \frac{\epsilon_2}{k_b} \right) - \sqrt{\frac{\epsilon_1}{k_b} \times \frac{\epsilon_2}{k_b}} \right] \qquad (18) $$

An atomic additivity model is used to calculate the dispersion parameters of each molecule, according to Eq. (19).

$$ \frac{\epsilon_i}{k_b} = \frac{1}{N_{\text{atoms}}} \sum_{\text{atom } j} \frac{\epsilon_{\text{atom } j}}{k_b} \qquad (19) $$

Here, $\epsilon_{\text{atom } j}/k_b$ represents the contribution of an atom j in the calculation of the dispersion parameter of a given molecule i. Further, $N_{\text{atoms}}$ denotes the total number of atoms for which $|\epsilon_{\text{atom } j}/k_b| > 0$. Moreover, the constant $w$ that appears in Eq. (18) takes different values depending on the binary system considered, as shown Eq. (20).

$$ w = \begin{cases} -w_0 \quad (\text{in K}^{-1}), & \text{if water} + \text{hb only - acceptor} \\ -w_0 \quad (\text{in K}^{-1}), & \text{if COOH} + \text{nhb or hb donor - acceptor} \\ -w_0 \quad (\text{in K}^{-1}), & \text{if water} + \text{COOH} \\ +w_0 \quad (\text{in K}^{-1}), & \text{otherwise} \end{cases} \qquad (20) $$

Where $w_0$ represents a universal constant of the model, also fitted to experimental data.

## 3. Technical details

In this section, we detail how the pure compound inputs were obtained in order to perform calculations in mixtures using the tc-PR/COSMO-SAC models. In addition, the procedure for optimizing the universal constants of each version of the COSMO-SAC models is detailed, along with an overview of the methodology to perform phase-equilibrium calculations and score the models according to the benchmarking metrics proposed by Jaubert et al. (Jaubert et al., 2020).

### 3.1. Pure compound inputs of the tc-PR/COSMO-SAC models

As mentioned above, one of the key advantages of COSMO-type models is their ability to describe phase equilibria and thermodynamic properties in mixtures, relying only on pure component information.

In this implementation of the COSMO-SAC models coupled with the tc-PR EoS, the resulting model requires the critical temperature and pressure of each component of the mixture, along with their respective Twu-91 $\alpha$-function and volume translation parameters. The critical properties were sourced from the DIPPR database, while optimized values for the Twu-91 $\alpha$-function and volume translation parameters were obtained from the work of Piña-Martinez et al. (Pina-Martinez et al., 2022).

For the activity coefficients calculations, the COSMO files were sourced from the publication of Bell et al. (Bell et al., 2020). These COSMO files were generated using DMol3 with the PBE functional

<caption>Table 1</caption>

Binary association code (BAC) and data distribution for a total of 32,727 experimental datapoints.

<table><thead><tr><td>Binary association code (BAC)</td><td>Data distribution</td><td>Description</td></tr></thead><tbody><tr><td>BAC 1 (NA + NA)</td><td>20.34 %</td><td rowspan="2">Mixtures in which association does not take place</td></tr><tr><td>BAC 2 (NA + HA)</td><td>14.77 %</td></tr><tr><td>BAC 3 (NA + HD)</td><td>6.21 %</td><td></td></tr><tr><td>BAC 4 (HD + HD) or (HA + HA)</td><td>8.13 %</td><td></td></tr><tr><td>BAC 5 (NA + SA)</td><td>7.93 %</td><td>Mixtures in which only self-association takes place</td></tr><tr><td>BAC 6 (HD + HA)</td><td>8.81 %</td><td>Mixtures in which only cross-association takes place</td></tr><tr><td>BAC 7 (HD + SA)</td><td>4.25 %</td><td rowspan="3">Mixtures in which both cross and self-association take place</td></tr><tr><td>BAC 8 (HA + SA)</td><td>14.09 %</td></tr><tr><td>BAC 9 (SA + SA)</td><td>15.46 %</td></tr></tbody></table>

(Perdew-Burke-Ernzerhof) for exchange-correlation within the density functional theory (DFT) framework, combined with the DNP basis set. All $\sigma$-profiles and dispersion parameters were calculated from these COSMO files. Additional details are given in the Supplementary Material.

### 3.2. Benchmarking database from Jaubert et al. (Jaubert et al., 2020)

The experimental data used for the optimization of the COSMO-SAC models and evaluation of their performance within the UHVM mixing rules were obtained from the databank published by Jaubert et al. (Jaubert et al., 2020). This comprehensive dataset contains information on 200 binary systems involving various chemical species and experimental thermodynamic data, including phase-equilibrium data at both low and high pressures (such as VLE, LLE, VLLE, azeotropic and critical points), and energetic properties like enthalpy ($h^M$) and heat capacity ($c_p^M$) changes on mixing.

The original database of 200 binary systems was not used in its entirety. Instead, 160 systems were retained in this work for the following reasons:

* Availability of the $\sigma$-profiles for both components of the binary mixture;

* Only components containing C, H, O, N, F, and Cl were considered, and systems including other elements were excluded. This is because the dispersion parameters were calibrated with a particular focus on these atoms.

In any case, the final dataset is uniformly distributed according to the association behavior of each mixture constituent, allowing for thorough model evaluation. Basically, a molecule can be classified as non-associating (NA), as exemplified by n-alkanes; hydrogen-donor (HD), including chlorinated hydrocarbons; hydrogen-acceptor (HA), such as aldehydes and ketones; or self-associating (SA), which has both hydrogen donor and hydrogen acceptor sites, as observed in water, alcohols and carboxylic acids.

The different combinations of these molecules characterize the so-called binary association codes (BACs), as denoted in Table 1. This table also shows the distribution of experimental data across the different BACs.

### 3.3. Optimizing the COSMO-SAC universal constants on a subset

A subset of 45 systems (5 per BAC) was employed to optimize the universal constants associated with the COSMO-SAC models. This database included selected data on liquid-vapor equilibrium at constant temperature. The optimization objective function (OF) was defined
based on the calculations of bubble pressures and vapor-phase compositions $y_1$ (for given temperature $T$ and liquid-phase composition $x_1$):

The coefficients of this expansion ($C_k^{RK}$) were determined using a polynomial regression technique on a set of data generated by the COSMO-SAC model. The dataset initially consisted of at least 5 points

$$ OF = \sum_{\substack{\text{exp data, } \\ 1 \le k \le N}} \sqrt{\left[ \frac{P_{\text{bubble}}^{\text{k,exp}} - P_{\text{bubble}}^{\text{calc}}(T_{\text{exp},k}, x_{1,\text{exp},k})}{P_{\text{bubble}}^{\text{k,exp}}} \right]^2 + \left[ y_1^{\text{k,exp}} - y_1^{\text{calc}}(T_{\text{exp},k}, x_{1,\text{exp},k}) \right]^2} \quad (21) $$

Where, $P_{\text{bubble}}^{\text{k,exp}}$ and $P_{\text{bubble}}^{\text{calc}}(T_{\text{exp},k}, x_{1,\text{exp},k})$ are the experimental and calculated bubble pressures of a given mixture with known temperature and composition (the one of the liquid phase), whereas $y_1^{\text{k,exp}}$ and $y_1^{\text{calc}}(T_{\text{exp},k}, x_{1,\text{exp},k})$ are the experimental and calculated compositions (in terms of molar fraction) in the vapor phase. The bubble pressures of the binary mixtures were calculated using a bubble point algorithm (Michelsen and Mollerup, 2007) implemented in in-house Matlab codes (r2024b version). The minimization of the objective function OF was performed with a Nelder-Mead optimizer. Initially, only the residual term parameters of the COSMO-SAC 2002 and 2010 versions were optimized (i.e., $\alpha^*$, $C_{hb}$, $\sigma_{hb}$, $c_{OH-OH}$, $c_{OH-OT}$ and $c_{OT-OT}$). It should be noted that the effective contact area ($a_{\text{eff}}$) was maintained identical to that used in generating the $\sigma$-profiles (refer to Supplementary Materials for more details).

Afterwards, the dispersion parameters ($w_0$ and $\varepsilon_{\text{atom}}/k_b$ of each atom) of the tc-PR/COSMO-SAC-dsp model were optimized, keeping the same constants of the residual term from the COSMO-SAC 2010 version fitted in the previous step. Originally, the $\varepsilon_{\text{atom}}/k_b$ parameters cover the atoms C, N, O, F, Cl, and H. In this study, additional parameters were included to better represent gas molecules such as CO$_2$, CO, N$_2$, and NH$_3$. For this reason, additional phase equilibrium data involving these gas molecules were incorporated to adjust the dispersion parameters specific to the corresponding gas compounds.

### 3.4. Calculating phase equilibrium data

A bubble point algorithm was previously chosen due to its relatively lower computational cost to calculate bubble pressures, making it a good option for optimization tasks. This is particularly important in the framework of COSMO-based models since they are computationally demanding due to the implicit solution of Eqs. (7) and (12). Nonetheless, after the optimization procedure, a flash-type algorithm was used to perform the phase-equilibrium calculations for all the data contained in the complete database. This type of phase-equilibrium calculation was selected as it offers greater robustness, particularly when dealing with more complex phase diagrams (which cannot be simulated using a simplified bubble point algorithm, especially when one component is in a supercritical state).

A flash-type algorithm often requires a greater number of iterations (especially near azeotropic or critical points), as well as the temperature and composition derivatives of the excess property ($g^E$) used in the mixing rule, which incorporates the COSMO-based model. These derivatives were calculated as follows:

* For the temperature derivatives of the excess function, numerical differentiation was applied using a central difference scheme.

* In order to accelerate the calculation of compositional derivatives (with respect to molar fraction vector, $x$) at a constant temperature, the Redlich-Kister expansion was considered:

$$ g^E(\mathbf{x}, T) = x_1 x_2 \sum_{k=0}^{m} C_k^{RK} (x_1 - x_2)^k \quad (22) $$

($x_1, g^E(T, x_1)$), with compositions $x_1$ evenly distributed between 0 and 1, and a second-order expansion ($m=2$) was used as a starting point. The fit was only accepted if the correlation coefficient ($R^2$) exceeded 0.99. If this threshold was not met, the number of datapoints from COSMO-SAC was increased, up to a maximum of 50 points. If this limit was reached without achieving an $R^2$ of 0.99, the order of the polynomial expansion ($m$) was increased. Note that increasing the order of the expansion is performed in a second step to avoid overfitting when a high-order polynomial is considered, which could negatively impact the accuracy of the compositional derivatives.

Overall, using the methodology described above, TP-flash calculations are accelerated by a factor of around 1,000 compared to implementations in which the COSMO-SAC model is directly used to compute activity coefficients at each iteration. This gain is especially significant when Newton-based solvers are employed, as they require composition derivatives.

### 3.5. Rating the results

As previously mentioned, we employed the Mark20 evaluation scale to benchmark the model performance, which was first proposed by Jaubert et al. (Jaubert et al., 2020) and later adopted in subsequent studies (Nikolaidis et al., 2023; Pina-Martinez et al., 2021). This scale assesses the model accuracy for different types of phase equilibrium and mixing properties data, assigning each of them a score between 0 and 20. The final score is obtained by averaging these individual evaluations, providing a global measure of the model performance. A step-by-step guide of this procedure is provided below:

<figure><img src="image_4.png" /><figcaption>Fig. 3. Illustration of ‘outside-the-model’ datapoints on a phase diagram at constant temperature.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 4. Qualitative analysis of scores within the MARK<sub>20</sub> scale.</figcaption></figure>

**Step 1:** we calculate the mean average percentage error (MAPE) between calculations and experimental data across the 9 BACs. This calculation is performed for each property considered: liquid-phase composition (x), vapor-phase composition (y), 3-phase pressure (P<sub>LLV</sub>), 3-phase composition (z<sub>LLV</sub>), critical pressure (P<sub>c</sub>), critical composition (x<sub>c</sub>), azeotropic pressure (P<sub>azeo</sub>), azeotropic composition (x<sub>azeo</sub>), enthalpy of mixing (h<sup>M</sup>), and heat capacity of mixing (c<sub>p</sub><sup>M</sup>). In the equations below, “prop” refers to one of the properties listed above.

$$MAPE_{prop}^{BAC} = \frac{1}{N_{\text{data}}} \sum_{i=1}^{N_{\text{data}}} \left| \frac{prop_i^{\text{exp}} - prop_i^{\text{calc}}}{prop_i^{\text{exp}}} \right| \quad (23)$$

**Step 2:** We assign a Mark value for each property within a BAC using the Equation (24). The value of α<sub>prop</sub> parameter that appears in this equation is 0.5 for all properties except critical points, enthalpies of mixing, and heat capacity of mixing, for which α<sub>prop</sub> is assigned values of 0.75, 0.25, and 0.1, respectively.

$$ \left\{ \begin{array}{l} \text{Mark}_{prop}^{BAC} = (20 - \alpha_{prop}) \text{MAPE}_{prop}^{BAC} \times \text{SR}_{prop}^{BAC} \\ \text{with:} \quad \text{SR}_{prop}^{BAC} = 1 - \frac{N_{\text{otm}} \text{ points}}{N_{\text{all}} \text{ points}} \end{array} \right. \quad (24) $$

Note that SR<sub>prop</sub><sup>BAC</sup> stands for the 'success ratio', which is a measure of the model to actually reproduce the experimental datapoints. In Fig. 3, experimental points for which the relative deviation cannot be calculated are classified as 'outside the model' (otm). As shown in Fig. 3, this situation may occur when the pressures at which experimental data are available exceed the azeotropic pressure predicted by the thermodynamic model. While these points are excluded from the MAPE calculation, they influence the final score through the 'success ratio,' defined as the proportion of out-of-model points relative to the total number of experimental datapoints.

**Step 3:** the MARK value assigned to an entire BAC is calculated by averaging the Mark values of all the properties. Note that properties that do not appear in the BAC considered are removed from the average calculation.

**Step 4:** the Mark20 is then calculated by category of association, as shown in Eq. (25). In the notation used, “N-A” refers to non-associative mixtures, “S-A” to mixtures in which only self-association occurs, “C-A” to mixtures in which only cross-association occurs, and “S-A + C-A” to mixture showing both types of association.

<caption>**Table 2**</caption>

Universal constants of COSMO-SAC models. The other parameters of the COSMO-SAC-dsp are equivalent to those of COSMO-SAC 2010. Note that the a<sub>eff</sub> (Å<sup>2</sup>) parameter was not re-optimized and remains the same as that used to generate the σ-profiles. The pairwise segment interaction constants (α', c<sub>hb</sub>, c<sub>OH-OH</sub>, c<sub>OH-OT</sub>, and c<sub>OT-OT</sub>) have units of kcal Å<sup>4</sup> mol<sup>-1</sup> e<sup>-2</sup>, while σ<sub>hb</sub> is given in e Å<sup>-2</sup>. Moreover, the dispersion constants (ε<sub>atom</sub>/k<sub>hb</sub>) are given in K.

<table><thead><tr><th colspan="2">COSMO-SAC 2002 (Lin and Sandler, 2002)</th></tr><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>a<sub>eff</sub></td><td>7.5</td></tr><tr><td>α'</td><td>18,168.61</td></tr><tr><td>C<sub>hb</sub></td><td>93,921.84</td></tr><tr><td>σ<sub>hb</sub></td><td>0.01041</td></tr></tbody><thead><tr><th colspan="2">COSMO-SAC 2010 (Hsieh et al., 2010)</th></tr><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>a<sub>eff</sub></td><td>7.25</td></tr><tr><td>α'</td><td>17,277.85</td></tr><tr><td>c<sub>OH-OH</sub></td><td>3,879.91</td></tr><tr><td>c<sub>OH-OT</sub></td><td>2,837.95</td></tr><tr><td>c<sub>OT-OT</sub></td><td>941.85</td></tr></tbody><thead><tr><th colspan="2">COSMO-SAC-dsp (Hsieh et al., 2014)</th></tr><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>w<sub>0</sub></td><td>0.5102</td></tr><tr><td>C(sp<sub>3</sub>)</td><td>338.49</td></tr><tr><td>C(sp<sub>2</sub>)</td><td>355.42</td></tr><tr><td>C(sp)</td><td>348.23</td></tr><tr><td>N(sp<sub>3</sub>)</td><td>247.28</td></tr><tr><td>N(sp<sub>2</sub>)</td><td>223.50</td></tr><tr><td>N(sp)</td><td>288.03</td></tr><tr><td>-O-</td><td>351.88</td></tr><tr><td>=O</td><td>211.75</td></tr><tr><td>F</td><td>254.44</td></tr><tr><td>Cl</td><td>352.14</td></tr><tr><td>H (H<sub>2</sub>O)</td><td>306.28</td></tr><tr><td>H (OH)</td><td>226.49</td></tr><tr><td>H (NH)</td><td>389.10</td></tr><tr><td>H (other)</td><td>0</td></tr><tr><td>=O (CO<sub>2</sub>) *</td><td>249.34</td></tr><tr><td>=[O + ] *</td><td>303.73</td></tr><tr><td>=[N] *</td><td>293.85</td></tr><tr><td>H (NH<sub>3</sub>) *</td><td>346.08</td></tr></tbody></table>

$$ \left\{ \begin{array}{l} \text{Mark}_{(N-A)} = \frac{1}{4} \sum_{i=1}^{4} \text{Mark}_{\text{BAC } i} \\ \text{Mark}_{(S-A)} = \text{Mark}_{\text{BAC } 5} \\ \text{Mark}_{(C-A)} = \text{Mark}_{\text{BAC } 6} \\ \text{Mark}_{(S-A+C-A)} = \frac{1}{3} \sum_{i=7}^{9} \text{Mark}_{\text{BAC } i} \end{array} \right. \quad (25) $$

**Step 5:** the model is assigned a global mark, considering all the BACs and association categories:

$$ \text{Mark}_{\text{all}} = \frac{\text{Mark}_{(N-A)} + \text{Mark}_{(S-A)} + \text{Mark}_{(C-A)} + \text{Mark}_{(S-A+C-A)}}{4} \quad (26) $$

In summary, the methodology assigns an overall score to thermodynamic models based on mean absolute percentage error (MAPE). This approach is useful because defining a single MAPE across all properties can be misleading; for instance, the MAPE of mixing molar heat capacities may appear large simply because the values are small, making small absolute deviations look disproportionately high. The Mark20 methodology mitigates these distortions when deriving an overall performance metric. In addition, in order to provide a qualitative

<footnote>* New parameters added in this work.</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 5. Global Mark₂₀ results for tc-PR/COSMO-SAC 2002, tc-PR/COSMO-SAC 2010 and tc-PR/COSMO-SAC-dsp. The results are compared with correlative equations of state, such as tc-PR/Wilson, tc-PR/NRTL, tc-PR/UNIQUAC, along with PC-SAFT. All models were evaluated on binary mixtures data sourced from (Jaubert et al., 2020).</figcaption></figure>

interpretation of the MARK20 results, we established a performance ranking for the models (Fig. 4) based on their scores. In general, correlative equations of state—i.e., those using binary interaction parameters fitted to experimental data—fall within the ‘average’ to ‘good’ range.

# 4. Results

## 4.1. General results

### 4.1.1. Parametrization of the COSMO-SAC models

Table 2 lists the optimized values of the universal constants associated with each version of the COSMO-SAC model considered in this study. The parameters used to initialize the optimization algorithm were taken from their respective original publications.

It is important to note that the parameters reported in Table 2 were obtained for use in the UHVM mixing rule within the PR EoS, and are not optimal for a $\gamma-\varphi$ framework. For $\gamma-\varphi$ calculations, the original constants of the respective COSMO-SAC models should be used. With respect to the choice between the $\gamma-\varphi$ and $\varphi-\varphi$ approaches, it is worth noting that both yield very similar results when the calculations are performed under subcritical conditions for all components of the mixture. A more detailed discussion of this point is provided in Supplementary Material. Given that this work focuses on the $\varphi-\varphi$ approach, the subsequent sections address this formulation exclusively.

### 4.1.2. Mark₂₀ results

After parameterizing each COSMO-SAC model version, the thermodynamic calculations were carried out for the full dataset of 160 systems. The general results, classified according to the Mark20 methodology described in section 3.3, are presented in Fig. 5, while more detailed results are given in the Supplementary Material (please refer to Appendix A for information on the organization of the supplementary material). Fig. 5 compares the tc-PR/COSMO-SAC models with correlative equations of state for which experimental data were used to fit binary interaction parameters. Among the equations of state considered are the tc-PR EoS, combined with various correlative activity coefficient models (Wilson, NRTL, and UNIQUAC) (Pina-Martinez et al., 2021), and the PC-SAFT EoS (Nikolaidis et al., 2023, 2021). Note that these correlative EoSs were also assessed on the database proposed by Jaubert et al. (Jaubert et al., 2020).

As expected, the cubic EoSs combined with COSMO-SAC models generally performed worse than the correlative models with per-system parameters fitted to experimental data. It is important to note, however,

<figure><img src="image_7.png" /><figcaption>Fig. 6. Mark₂₀ results for tc-PR/COSMO-SAC 2002, tc-PR/COSMO-SAC 2010 and tc-PR/COSMO-SAC-dsp across each binary association code (see Table 1). ‘NA’ stands for non-associating mixtures, ‘SA’ concerns to mixtures for which only self-association take place, ‘CA’ represents mixtures for which only cross-association take place, while ‘SA + CA’ stands for mixtures where both self- and cross-association take place.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Mark₂₀ results for tc-PR/COSMO-SAC 2002, tc-PR/COSMO-SAC 2010 and tc-PR/COSMO-SAC-dsp across each property considered: liquid-phase composition (x), vapor-phase composition (y), 3-phase pressure (pllv), 3-phase composition (xllv), critical pressure (pcrit), critical composition (xcrit), azeotropic pressure (paz), azeotropic composition (xaz), enthalpy of mixing (hM), and heat capacity of mixing (cpm).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Calculation of phase envelopes for BAC 5 systems using the tc-PR/COSMO-SAC 2002 and tc-PR/COSMO-SAC 2010 models.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 9. Calculation of phase envelopes for BAC 6 systems using the tc-PR/ COSMO-SAC 2002 and tc-PR/COSMO-SAC 2010 models. In (a), the phase envelope for the mixture chlorodifluoromethane + dimethyl ether at 323 K, and in (b), the phase envelope for the mixture acetone + chloroform at 323 K.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Calculation of phase envelopes for acetone + alcohol and acetone + water systems using the tc-PR/COSMO-SAC 2002 and tc-PR/COSMO-SAC 2010 models.</figcaption></figure>

that the tc-PR/COSMO-SAC models are based on predictive mixing rules and their performance depends on COSMO files derived from DFT calculations. This makes the COSMO-based equations of state applicable to a much wider variety of mixtures, since experimental data is not required to evaluate model parameters.

Delving deeper into the results, Fig. 6 and Fig. 7 show that the

performance of the tc-PR EoS across the different thermodynamic properties and equilibrium data improves with the use of more advanced COSMO-SAC versions (2010 and dp5), as these models provide a more refined treatment of COSMO files for the calculation of activity coefficients. However, greater attention must be given to whether a dispersion term should be included for BAC1 mixtures, as these systems
<figure><img src="image_12.png" /><figcaption>Fig. 11. Calculation of phase envelopes for aqueous BAC 9 systems using the tc-PR/COSMO-SAC 2002 and tc-PR/COSMO-SAC 2010 models.</figcaption></figure>

often behave close to ideality, with activity coefficients near unity. This helps explain why the dpst model performs less accurately than the other versions for this BAC.

Overall, phase diagrams (even those exhibiting azeotropic behavior), are well described by the tc-PR/COSMO-SAC EoSs, particularly when using COSMO-SAC 2010 and its dpst variant. Nonetheless, all models exhibited rather unsatisfactory performance regarding the calculation of critical points, composition of three-phase lines, and energetic properties such as enthalpy and heat capacity changes on mixing.

For critical points and three-phase lines, the poor performance can be attributed to the low success ratio obtained across most BACs, generally below 50% (Refer to Supplementary Material). With regard to energetic properties, the low grades on the Mark₂₀ scale can be attributed to the model formulation, which lacks functional forms to consistently describe the temperature-dependence of segment-segment interactions. Additionally, the entire dataset covers a wide temperature range, varying from a few kelvins up to nearly 1000 K, depending on the binary system under consideration. Maintaining accuracy under such wide conditions is challenging for a predictive model based on only a few universal constants.

## 4.2. Phase equilibria

### 4.2.1. From COSMO-SAC 2002 to COSMO-SAC-2010

For non-associating systems in BACs 1 to 4, the tc-PR/COSMO-SAC 2002 and tc-PR/COSMO-SAC 2010 versions exhibit similar performance (Fig. 6), as their formulations are very similar. For such systems, the observed differences in the results arise primarily from the values of the universal parameters. However, for systems involving hydrogen bonding (BACs 5 to 9), the use of the COSMO-SAC-2010 model in the mixing rule offers a clear advantage.

To illustrate this, let us first consider BAC 5, which consists of mixtures of non-associating and self-associating molecules, such as hydrocarbons and alcohols. The hydrogen-bonding dynamics in these systems is complex because the self-associating molecules form hydrogen bonds with each other that are disrupted by the presence of the non-associating molecules. This explains why these mixtures are often partially miscible. In general, COSMO-type models encounter challenges in representing this type of behavior, particularly in their more simplified formulations, as illustrated in Fig. 8. For example, the figure shows a clear underestimation of the deviation from ideality when the tc-PR/COSMO-SAC 2002 version is used to describe alkane + alcohol mixtures.

The enhanced treatment of hydrogen bonds highlights the superiority of the tc-PR/COSMO-SAC 2010 model for BAC 6. This category
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Calculation of isothermal phase diagrams for the system carbon dioxide + ethane using the tc-PR/COSMO-SAC 2010 and tc-PR/COSMO-SAC-dsp models. Note that, in the absence of the dispersion term, the tc-PR/COSMO-SAC 2010 model fails to predict the existence of azeotropic points.</figcaption></figure>

refers to mixtures where the components do not form hydrogen bonds in their pure state, but form cross-hydrogen bonds when mixed. These cross-hydrogen bonds increase intermolecular attraction, resulting in negative deviations from ideality. Qualitatively, these negative deviations are well captured by the tc-PR/COSMO-SAC models, as illustrated in Fig. 9. Although both models tend to overestimate the negative deviations from ideality, the tc-PR/COSMO-SAC 2010 model clearly outperforms when used with the UHVM mixing rules (Fig. 10).

A general improvement was also observed in the associating systems comprising BACs 7 to 9 when the COSMO-SAC 2010 version is used. Considering the acetone + alcohol and acetone + water mixtures as illustrative examples, it can be observed that the updated formulation proposed by Hsieh et al. (Hsieh et al., 2010) improves the results considerably.

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Calculation of isothermal phase diagrams for the systems acetone + benzene (BAC 2) and difluoromethane + propane (BAC 3) using the tc-PR/COSMO-SAC 2010 and tc-PR/COSMO-SAC-dsp models.</figcaption></figure>

As COSMO-SAC 2010 introduces hydrogen bond constants (c<sub>OT-OT</sub>, c<sub>OH-OH</sub>, and c<sub>OH-OT</sub>) tailored to the specific interaction between the different σ-profiles (nhb, OH, and OT), the acetone-acetone (OT-OT), alcohol-alcohol (OH-OH), and acetone-alcohol (OT-OH) interactions vary their intensity in a more realistic way. Using a single hydrogen bonding constant maintains the hydrogen bond interactions between molecules of the same species and between different species at a more uniform level. This explains why the bubble curve obtained with the tc-PR/COSMO-SAC 2002 model approaches the behavior of an ideal solution (i.e., when the bubble line is the straight line connecting the saturation points of the two pure components, meaning that the interactions between molecules of different species have similar magnitude to the interactions in pure state).

Systems with complex hydrogen-bond networks (especially BAC 9) are also challenging to be described by predictive models, such as COSMO-type ones. Aqueous mixtures, in particular, represent a well-known limitation, especially when water is combined with another self-associating molecule such as an alcohol, nitrile, or amine. In such cases, complex hydrogen-bond networks arise, involving both self- and cross-association. Fig. 11 illustrates this behavior with four representative examples.

It is worth noting that at higher temperatures the influence of hydrogen bonding diminishes, which explains the good agreement observed in Fig. 11, panel (d). This is partially due to the use of a simplified temperature-dependent function on hydrogen-bond interactions.

### 4.2.2. Effect of the dispersion term

As discussed above, Hsieh and co-authors improved the COSMO-SAC-2010 model by adding a dispersion term in the form of a simplified 1-parameter Margules model where the constant is determined from the molecular dispersion parameter of the components, which, in turn, is calculated using an atomic additivity model (Hsieh et al., 2014). Comparing with our previous results, an overall improvement in the tc-PR/COSMO-SAC model across all BACs is observed.

From a modeling standpoint, the electrostatic-based formulation of COSMO-type models is often insufficient to describe the ideality deviations of systems containing at least one nonpolar (or nearly nonpolar) compound, such as mixtures belonging to BACs 1 to 4. An example is presented in Fig. 12 for the binary mixture carbon dioxide + ethane, which exhibits significant deviations from ideality attributable to dispersion effects. Additional examples to reinforce this point are given in Fig. 13.

For mixtures from BAC 6, in particular, the dispersion term attenuates the negative deviations from ideality, thereby aligning the model

<figure><img src="image_15.png" /></figure>
<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 14. Calculation of phase envelopes for BAC 6 systems using the tc-PR/COSMO-SAC 2010 and tc-PR/COSMO-SAC-dsp models.</figcaption></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 15. Residual activity coefficient ($\gamma$) as a function of the (composition in terms of molar fraction, $x_1$) for the systems methanol + water and carbon dioxide + ethane. The calculation of the residual activity coefficients was performed using COSMO-SAC 2010, while the dispersion-based activity coefficients were calculated with the 1-parameter Margules model used in the COSMO-SAC-dsp implementation.</figcaption></figure>

more closely with experimental data. In fact, the dispersion term reduces systematic errors observed, for example, in the underprediction of azeotropic pressures, as illustrated by Fig. 14.

For BACs 5, 7, 8, and 9, which exhibit more strongly connected hydrogen-bond networks (which can be attributed to the presence of self-associating molecules), the contribution of the dispersive term be-

comes less significant. As a result, the residual activity coefficient takes values significantly different from unity, whereas the dispersive contribution remains closer to one. This explains why the MARK₂₀ scores for these BACs remained nearly unchanged when moving from the tc-PR/COSMO-SAC 2010 to the tc-PR/COSMO-SAC-dsp version. To illustrate this, Fig. 15 presents the residual and dispersive activity
<figure><img src="image_22.png" /><figcaption>**Fig. 16.** Enthalpy changes of mixing for nine binary systems calculated with the tc-PR/COSMO-SAC 2002, tc-PR/COSMO-SAC 2010, and tc-PR/COSMO-SAC-dsp models.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. 17.** Isobaric heat capacities of mixing for three binary mixtures calculated with the tc-PR/COSMO-SAC 2002, tc-PR/COSMO-SAC 2010, and tc-PR/COSMO-SAC-dsp models.</figcaption></figure>

coefficients calculated within the COSMO-SAC-dsp framework. For an associating system from BAC 9 (methanol + water), the absolute value of the residual ln(γ) is considerably larger than that of the dispersive contribution. In contrast, for a non-associating system such as CO₂ + ethane, the dispersive term becomes more relevant.

### 4.3. Energetic properties (enthalpy and isobaric heat capacity of mixing)

Regarding the dependence of the COSMO models on temperature, all tc-PR/COSMO-SAC models, overall, showed a poor correlation with experimental data for the calculation of enthalpies and heat capacity of mixing, with most BACs scoring below 5/20 for these energetic properties. This highlights a well-known limitation of COSMO models: their limited precision in capturing temperature effects. As previously
reported in the literature (Chang and Siepmann, 2024; Constantinescu et al., 2009), this limitation directly impacts the calculation accuracy of derivative properties.

For illustration purposes, Fig. 16 reports the enthalpy of mixing results for nine representative binary mixtures, while Fig. 17 presents the outcomes for the calculation of heat capacities of mixing involving 3 binary systems.

In most cases, the predictions obtained for these properties were only qualitative, indicating that the tc-PR/COSMO-SAC models were able to capture general trends but not to produce accurate predictions from a quantitative viewpoint. For mixtures where hydrogen bonding does not take place (e.g., hydrocarbons with halogenated compounds, ketones, ethers and esters), the tc-PR/COSMO-SAC 2002 and tc-PR/COSMO-SAC-dsp models yield better results. On the other hand, concerning systems containing self-associating molecules such as water, alcohols, and carboxylic acids, the 2002 version performed poorly when compared to newer versions.

It should be noted that achieving simultaneously accurate descriptions of both phase-equilibrium and energy-related properties with a single equation of state represents a well-known challenge in thermodynamics, especially across wide temperature ranges (Jaubert et al., 2022). Even correlative models relying on binary interaction parameters specifically fitted to experimental data, often face limitations in this regard. Naturally, the complexity of this task increases further when the mixing rule is based on pure predictions of activity coefficients calculated with COSMO-SAC models, which generally do not account for this type of energetic property during the optimization of their few universal constants.

In this context, the results reported here reflect the overall performance in predicting energetic properties rather than a strict quantitative validation, while also pointing out to possible improvements through the incorporation of consistent temperature-dependent functions for electrostatic, hydrogen-bonding, and dispersive interactions. This implies the introduction of new universal constants, which should be simultaneously adjusted against both phase-equilibrium data (at low and high temperatures) and derivative property data.

## 5. Conclusions

This paper presents a systematic analysis of predictive EoS obtained by combining the tc-PR EoS with COSMO-SAC models through advanced (UHVM) mixing rules with the view to predict phase equilibria and mixing properties. It was observed that the tc-PR/COSMO-SAC model successfully reproduced experimental data, at least qualitatively, and in several cases quantitatively too, even for complex diagrams involving, azeotropic points, triphasic lines, and critical points. However, when it comes to the energetic properties of mixing, such as enthalpies and isobaric heat capacities, the tc-PR/COSMO-SAC models performed poorly for many of the binary systems studied.

To improve the COSMO-type approach within mixing rules, we recommend analyzing the following aspects in future research:

* It is recommended to use, at least, an additional σ-profile to represent interactions with water molecules (nhb, OT, OH, w), rather than relying solely on the three σ-profiles used in the COSMO-SAC 2010 formulation (nhb, OT, OH). As Prinos and Voutsas have suggested (Prinos and Voutsas, 2025), this slight modification of the model will enhance its accuracy in capturing water interactions, while enabling a more focused parametrization of OH-OH and OH-OT interactions in the absence of aqueous systems.

* It is also recommended that a detailed analysis of the functions used to model the temperature dependence of electrostatic, hydrogen bonding, and dispersion terms is conducted, with the aim of enhancing the model accuracy, especially for predicting energetic properties of mixing. Although temperature-dependent functions have been proposed in the literature for electrostatic (Hsieh et al.,

2010), hydrogen-bonding, and dispersion parameters (Klamt and Eckert, 2000), these implementations have primarily focused on activity coefficients and phase equilibrium, with energetic properties generally not addressed.

* Additionally, given the important role of dispersion, implementing a more comprehensive formulation for dispersive interactions within the contact energy term, as described in references (Grigorash et al., 2025, 2024), would further improve the overall performance of the tc-PR/COSMO-SAC model.

Such improvements could significantly enhance the use of COSMO-type models within a φ-φ framework, even at low pressures. This advantage is reinforced by the availability of accurate machine learning models for estimating the input parameters of cubic equations of state (including critical temperatures, critical pressures, and acentric factors (Biswas et al., 2023; Bounaceur et al., 2025; Carande et al., 2015)), which renders the approach purely predictive.

## Funding sources

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Project: 101003318 – BioSCoPe). High-performance computing resources were provided by IDRIS under the allocation AD010812434R4 made by GENCI and also by the EXPLOR center hosted by the University of Lorraine.

## CRediT authorship contribution statement

**Francisco Paes:** Writing – review & editing, Writing – original draft, Visualization, Methodology, Investigation, Conceptualization. **Baptiste Sirjean:** Writing – review & editing, Validation, Supervision, Methodology, Conceptualization. **Jean-Noël Jaubert:** Writing – review & editing, Validation, Supervision, Methodology, Conceptualization. **Romain Privat:** Writing – review & editing, Validation, Supervision, Methodology, Conceptualization.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.ces.2026.123300.

## Data availability

Data will be made available on request.

## References

Bell, I.H., Mickoleit, E., Hsieh, C.-M., Lin, S.-T., Vrabec, J., Breitkopf, C., Jäger, A., 2020. A benchmark open-source implementation of COSMO-SAC. J. Chem. Theory Comput. 16, 2635–2646. https://doi.org/10.1021/acs.jctc.9b01016.

Biswas, S., Chung, Y., Ramirez, J., Wu, H., Green, W.H., 2023. Predicting critical properties and acentric factors of fluids using multitask machine learning. J. Chem. Inf. Model. 63, 4574–4588. https://doi.org/10.1021/acs.jcim.3c00546.

Bounaceur, R., Paes, F., Privat, R., Jaubert, J.-N., 2025. AI-powered prediction of critical properties and boiling points: a hybrid ensemble learning and QSPR approach. J. Cheminformatics 17, 132. https://doi.org/10.1186/s13321-025-01062-9.

Carande, W.H., Kazakov, A., Muzny, C., Frenkel, M., 2015. Quantitative structure–property relationship predictions of critical properties and acentric factors for pure compounds. J. Chem. Eng. Data 60, 1377–1387. https://doi.org/10.1021/je501093v.

Chang, C.-K., Siepmann, J.I., 2024. Enthalpy of Mixing and liquid structure for mixtures of primary and secondary alcohols: molecular simulations and COSMO-SAC
calculations. J. Chem. Eng. Data 69, 509-525. https://doi.org/10.1021/acs.jced.3c00415.

Constantinescu, D., Klamt, A., Geană, D., 2005. Vapor-liquid equilibrium prediction at high pressures using activity coefficients at infinite dilution from COSMO-RS-type methods. Fluid Phase Equilib. 231, 231-238. https://doi.org/10.1016/j.fluid.2005.01.014.

Constantinescu, D., Rarey, J., Gmehling, J., 2009. Application of COSMO-RS type models to the prediction of excess enthalpies. Ind. Eng. Chem. Res. 48, 8710-8725. https://doi.org/10.1021/ie900315p.

Costa, C.T.O.G., Tavares, F.W., Secchi, A.R., 2016a. Equation of state based on the hole-lattice theory and surface-charge density (COSMO): part A – pure compounds. Fluid Phase Equilib. 409, 472-481. https://doi.org/10.1016/j.fluid.2015.11.010.

Costa, C.T.O.G., Tavares, F.W., Secchi, A.R., 2016b. Equation of state based on the hole-lattice theory and surface-charge density (COSMO): part B – vapor-liquid equilibrium for mixtures. Fluid Phase Equilib. 419, 1-10. https://doi.org/10.1016/j.fluid.2016.03.005.

De Hemptinne, J.-C., Kontogeorgis, G.M., Dohrn, R., Economou, I.G., Ten Kate, A., Kuitunen, S., Fele Zilnik, L., De Angelis, M.G., Vesovic, V., 2022. A view on the future of applied thermodynamics. Ind. Eng. Chem. Res. 61, 14664-14680. https://doi.org/10.1021/acs.iecr.2c01906.

Eckert, F., Klamt, A., 2002. Fast solvent screening via quantum chemistry: COSMO-RS approach. AIChE J. 48, 369-385. https://doi.org/10.1002/aic.690480220.

Gerlach, T., Müller, S., de Castilla, A.G., Smirnova, I., 2021. An open source COSMO-RS implementation and parameterization supporting the efficient implementation of multiple segment descriptors. Fluid Phase Equilib. 560, 113472. https://doi.org/10.1016/j.fluid.2022.113472.

Grigorash, D., Müller, S., Paricaud, P., Stenby, E.H., Smirnova, I., Yan, W., 2024. A comprehensive approach to incorporating intermolecular dispersion into the openCOSMO-RS model. part 1: halocarbons. https://doi.org/10.48550/ARXIV.2406.05244.

Grigorash, D., Müller, S., Heid, E., Neese, F., Liakos, D., Riplinger, C., García-Ratés, M., Paricaud, P., Stenby, E.H., Smirnova, I., Yan, W., 2025. A comprehensive approach to incorporating intermolecular dispersion into the openCOSMO-RS model. Part 2: Atomic polarizabilities. https://doi.org/10.48550/ARXIV.2502.08520.

Hsieh, M.-K., Lin, S.-T., 2012. Effect of mixing rule boundary conditions on high pressure (liquid + liquid) equilibrium prediction. J. Chem. Thermodyn. 47, 33-41. https://doi.org/10.1016/j.jct.2011.09.023.

Hsieh, C.-M., Sandler, S.I., Lin, S.-T., 2010. Improvements of COSMO-SAC for vapor-liquid and liquid-liquid equilibrium predictions. Fluid Phase Equilib. 297, 90-97. https://doi.org/10.1016/j.fluid.2010.06.011.

Hsieh, C.-M., Lin, S.-T., Vrabec, J., 2014. Considering the dispersive interactions in the COSMO-SAC model for more accurate predictions of fluid phase behavior. Fluid Phase Equilib. 367, 109-116. https://doi.org/10.1016/j.fluid.2014.01.032.

Huron, M.-J., Vidal, J., 1979. New mixing rules in simple equations of state for representing vapour-liquid equilibria of strongly non-ideal mixtures. Fluid Phase Equilib. 3, 255-271. https://doi.org/10.1016/0378-3812(79)80001-1.

Jaubert, J.-N., Privat, R., Le Guennec, Y., Coniglio, L., 2016. Note on the properties altered by application of a Pénéeloux-type volume translation to an equation of state. Fluid Phase Equilib. 419, 88-95. https://doi.org/10.1016/j.fluid.2016.03.012.

Jaubert, J.-N., Le Guennec, Y., Pina-Martinez, A., Ramirez-Velez, N., Lasala, S., Schmid, B., Nikolaidis, I.K., Economou, I.G., Privat, R., 2020. Benchmark database containing binary-system-high-quality-certified data for cross-comparing thermodynamic models and assessing their accuracy. Ind. Eng. Chem. Res. 59, 14981-15027. https://doi.org/10.1021/acs.iecr.0c01734.

Jaubert, J.-N., Qian, J.-W., Lasala, S., Privat, R., 2022. The impressive impact of including enthalpy and heat capacity of mixing data when parameterising equations of state. application to the development of the E-PPR78 (Enhanced-Predictive-Peng-Robinson-78) model. Fluid Phase Equilib. 560, 113456. https://doi.org/10.1016/j.fluid.2022.113456.

Klamt, A., 1995. Conductor-like screening model for real solvents: a new approach to the quantitative calculation of solvation phenomena. J. Phys. Chem. 99, 2224-2235. https://doi.org/10.1021/j100007a062.

Klamt, A., 2005. COSMO-RS: from quantum chemistry to fluid phase thermodynamics and drug design, 1st ed. Elsevier, Amsterdam Boston.

Klamt, A., Eckert, F., 2000. COSMO-RS: a novel and efficient method for a priori prediction of thermophysical data of liquids. Fluid Phase Equilib. 172, 43-72. https://doi.org/10.1016/S0378-3812(00)00357-5.

Klamt, A., Eckert, F., Arlt, W., 2010. COSMO-RS: an alternative to simulation for calculating thermodynamic properties of liquid mixtures. Annu. Rev. Chem. Biomol. Eng. 1, 101-122. https://doi.org/10.1146/annurev-chembioeng-073009-100903.

Kontogeorgis, G.M., Dohrn, R., Economou, I.G., de Hemptinne, J.-C., ten Kate, A., Kuitunen, S., Mooijer, M., Zilnik, L.F., Vesovic, V., 2021. Industrial requirements for thermodynamic and transport properties: 2020. Ind. Eng. Chem. Res. 60, 4987-5013. https://doi.org/10.1021/acs.iecr.0c05356.

Lee, M.-T., Lin, S.-T., 2007. Prediction of mixture vapor-liquid equilibrium from the combined use of Peng-Robinson equation of state and COSMO-SAC activity coefficient model through the Wong-Sandler mixing rule. Fluid Phase Equilib. 254, 28-34. https://doi.org/10.1016/j.fluid.2007.02.012.

Le Guennec, Y., Lasala, S., Privat, R., Jaubert, J.-N., 2016a. A consistency test for α-functions of cubic equations of state. Fluid Phase Equilib. 427, 513-538. https://doi.org/10.1016/j.fluid.2016.07.026.

Le Guennec, Y., Privat, R., Jaubert, J.-N., 2016b. Development of the translated-consistent tc-PR and tc-RK cubic equations of state for a safe and accurate prediction of volumetric, energetic and saturation properties of pure compounds in the sub- and super-critical domains. Fluid Phase Equilib. 429, 301-312. https://doi.org/10.1016/j.fluid.2016.09.003.

Leonhard, K., Veverka, J., Lucas, K., 2009. A comparison of mixing rules for the combination of COSMO-RS and the Peng-Robinson equation of state. Fluid Phase Equilib. 275, 105-115. https://doi.org/10.1016/j.fluid.2008.09.016.

Lin, S.-T., Sandler, S.I., 2002. A priori phase equilibrium prediction from a segment contribution solvation model. Ind. Eng. Chem. Res. 41, 899-913. https://doi.org/10.1021/ie001047w.

Magoulas, K., Tassios, D., 1990. Thermophysical properties of n-Alkanes from C1 to C20 and their prediction for higher ones. Fluid Phase Equilib. 56, 119-140. https://doi.org/10.1016/0378-3812(90)85098-U.

Mambo-Lomba, D., Paricaud, P., 2021. Predictions of thermodynamic properties and phase equilibria of refrigerant systems with COSMO approaches. Int. J. Refrig. 124, 52-63. https://doi.org/10.1016/j.ijrefrig.2020.11.005.

Mikelsen, M.L., Mollerup, J.M., 2007. Thermodynamic models: fundamentals & computational aspects, 2. ed. Tie-Line Publications, Holte.

Mullins, E., Oldland, R., Liu, Y.A., Wang, S., Sandler, S.I., Chen, C.-C., Zwolak, M., Seavey, K.C., 2006. Sigma-profile database for using COSMO-based thermodynamic methods. Ind. Eng. Chem. Res. 45, 4389-4415. https://doi.org/10.1021/ie06370h.

Nikolaidis, I.K., Privat, R., Jaubert, J.-N., Economou, I.G., 2021. Assessment of the perturbed chain-statistical associating fluid theory equation of state against a benchmark database of high-quality binary-system data. Ind. Eng. Chem. Res. 60, 8935-8946. https://doi.org/10.1021/acs.iecr.1c01234.

Nikolaidis, I.K., Privat, R., Jaubert, J.-N., Economou, I.G., 2023. The effect of binary interaction parameters and induced association on the performance of the PC-SAPT equation of state: validation against a benchmark database of high-quality binary-system data. J. Chem. Eng. Data. https://doi.org/10.1021/acs.jced.2c00781.

Peneloux, A., Abdoul, W., Rauzy, E., 1989. Excess functions and equations of state. Fluid Phase Equilib. 47, 115-132. https://doi.org/10.1016/0378-3812(89)80172-4.

Pina-Martinez, A., Privat, R., Nikolaidis, I.K., Economou, I.G., Jaubert, J.-N., 2021. What is the optimal activity coefficient model to be combined with the translated-consistent Peng-Robinson equation of state through advanced mixing rules? Cross-comparison and grading of the Wilson, UNIQUAC, and NRTL a^ models against a benchmark database involving 200 binary systems. Ind. Eng. Chem. Res. 60, 17228-17247. https://doi.org/10.1021/acs.iecr.1c03003.

Pina-Martinez, A., Privat, R., Jaubert, J., 2022. Use of 300,000 pseudo-experimental data over 1800 pure fluids to assess the performance of four cubic equations of state: SRK, PR, tc-RK, and tc-PR. AIChE J 68, e17518. https://doi.org/10.1002/aic.17518.

Pinos, N., Voutsas, E., 2025. Phase equilibrium calculations at low and high pressures with a modified COSMO-SAC model. Fluid Phase Equilib. 589, 114277. https://doi.org/10.1016/j.fluid.2024.114277.

Privat, R., Jaubert, J.-N., Le Guennec, Y., 2016. Incorporation of a volume translation in an equation of state for fluid mixtures: which combining rule? Which effect on properties of mixing? Fluid Phase Equilib. 427, 414-420. https://doi.org/10.1016/j.fluid.2016.07.035.

Privat, R., Jaubert, J.-N., Kontogeorgis, G.M., 2025. Let us rethink advanced mixing rules for cubic equations of state. Fluid Phase Equilib. 114455. https://doi.org/10.1016/j.fluid.2025.114455.

Shimoyama, Y., Iwai, Y., Takada, S., Arai, Y., Tsuji, T., Hiaki, T., 2006. Prediction of phase equilibria for mixtures containing water, hydrocarbons and alcohols at high temperatures and pressures by cubic equation of state with GE type mixing rule based on COSMO-RS. Fluid Phase Equilib. 243, 183-192. https://doi.org/10.1016/j.fluid.2006.03.007.

Soares, R.D.e.P., Baladão, L.F., Staudt, P.B., 2019. A pairwise surface contact equation of state: COSMO-SAC-Phi. Fluid Phase Equilib. 488, 13-26. https://doi.org/10.1016/j.fluid.2019.01.015.

Vera, J.H., Prausnitz, J.M., 1972. Generalized van der Waals theory for dense fluids. Chem. Eng. J. 3, 1-13. https://doi.org/10.1016/0300-9467(72)85001-9.

Voutsas, E., Louli, V., Boukouvalas, C., Magoulas, K., Tassios, D., 2006. Thermodynamic property calculations with the universal mixing rule for EoS/G^ models: results with the Peng-Robinson EoS and a UNIFAC model. Fluid Phase Equilib. 241, 216-228. https://doi.org/10.1016/j.fluid.2005.12.028.

Wang, J., Sa, T., Zhu, L., 2023. Vapor-liquid equilibrium predictions of refrigerant systems using COSMO based Gex-EoS methods. Fluid Phase Equilib. 563, 113584. https://doi.org/10.1016/j.fluid.2022.113584.

Further reading

Ahmed, N.C.-B., Negadi, L., Mokbel, I., Jose, J., 2011. Phase equilibrium properties of binary aqueous solutions containing ethanediamine, 1,2-diaminopropane, 1,3-diaminopropane, or 1,4-diaminobutane at several temperatures. J. Chem. Thermodyn. 43, 719-724. https://doi.org/10.1016/j.jct.2010.12.010.

Barr-David, F., Dodge, B.F., 1959. Vapor-liquid equilibrium at high pressures, the systems ethanol-water and 2-propanol-water. J. Chem. Eng. Data 4, 107-121. https://doi.org/10.1021/je60002a003.

Bayles, J.W., Letcher, T.M., 1971. Thermodynamics of some binary liquid mixtures containing aliphatic amines. J. Chem. Eng. Data 16, 266-271. https://doi.org/10.1021/je60050a026.

Beare, W.G., McVicar, G.A., Ferguson, J.B., 1930. The determination of vapor and liquid compositions in binary systems. J. Phys. Chem. 34, 1310-1318. https://doi.org/10.1021/j150312a015.

Belabbaci, A., Villamanań, R.M., Negadi, L., Martín, M.C., 2012. Vapor-liquid equilibria of binary mixtures containing 2-butanol and hydrocarbons at 313.15 K. J. Chem. Eng. Data 57, 982-987. https://doi.org/10.1021/je201284y.

Bendiab, H., Roux-Desgranges, G., Roux, A.H., Grolier, J.-P.-E., Patterson, D., 1994. Excess heat capacities of ternary systems containing chlorobenzene or
chloronaphthalene. J. Solut. Chem. 23, 307–323. https://doi.org/10.1007/BF00973552.

Bobbo, S., Camporese, R., Stryjek, R., 1998. (Vapour + liquid) equilibrium measurements and correlations of the refrigerant mixture dimethyl ether (RE170) + 1,1,1,3,3,3-hexafluoropropane (R236fa) at the temperatures (303.68 and 323.75) K. J. Chem. Thermodyn. 30, 1041–1046. https://doi.org/10.1006/jctc.1998.0369.

Chadha, R., Tripathi, A.D., 1995. Excess molar enthalpies of 1,1,2,2-tetrachloroethane + 2-methylfuran, + tetrahydrofuran, + 1,4-dioxane, and + cyclopentanone at 308.15 and 318.15 K. J. Chem. Eng. Data 40, 645–646. https://doi.org/10.1021/je00019a024.

Cochran, S.W., Holste, J.C., Marsh, K.N., Gammon, B.E., Hall, K.R., 1993. Enthalpies of water+ methanol mixtures between 180 K and 320 K. Fluid Phase Equilib. 88, 171–181.

Coquelet, C., Chareton, A., Valtz, A., Baba-Ahmed, A., Richon, D., 2003. Vapor–liquid equilibrium data for the azeotropic difluoromethane + propane system at temperatures from 294.83 to 343.26 K and pressures up to 5.4 MPa. J. Chem. Eng. Data 48, 317–323. https://doi.org/10.1021/je020115d.

Drobig, J., Lindenau, E., 1966. A new apparatus to determine the isothermal liquid-vapor equilibrium. Wiss Z Tech Hochsch Chem'carl Schorlemmer "leuna-Mersebg. 8, 144–147.

Fenclova, D., Dohnal, V., Perez-Casas, S., Frigolet, C., Costas, M., 1995. Excess heat capacities, excess enthalpies, and excess volumes of trichloromethane+ some oxygenated solvents or+ the corresponding homomorphic hydrocarbons at 298.15 K. ELDATA Int. Electron J. Phys.-Chem. Data 1, 217–238.

Fredenslund, A., Mollerup, J., 1974. Measurement and prediction of equilibrium ratios for the C2H6+ CO2 system. J. Chem. Soc. Faraday Trans. 1 Phys. Chem. Condens. Phases 70, 1653. https://doi.org/10.1039/f19747001653.

Gmehling, J., Meents, B., 1992. Excess enthalpies of some alkanol (C1-C4)+ alkane or+ 1-alkene (C6, C8) mixtures. Int. Data Ser. Sel. Data Mix. Ser. A 20, 153–160.

Griswold, J., 1952. Phase-equilibria of the acetone-methanol-water system from 100°C into the critical region. Chem. Eng. Progr. Symp. Syr. 48, 18–34.

Hirobe, H., 1926. Thermochemical studies. J. Fac. Sci. Univ. Tokyo Sect. 1A Math. 1, 155–222.

Hu, J., Tamura, K., Murakami, S., 1997. Excess thermodynamic properties of binary mixtures of ethyl formate with benzene, ethanol, and 2,2,2-trifluoroethan-1-ol at 298.15 K. Fluid Phase Equilib. 131, 197–212. https://doi.org/10.1016/S0378-3812 (96)03230-X.

Im, J., Lee, G., Lee, J., Kim, H., 2007. Vapor–liquid equilibria of the 1,1,1-trifluoroethane (HFC-143a) + dimethyl ether (DME) system. Fluid Phase Equilib. 251, 59–62. https://doi.org/10.1016/j.fluid.2006.10.017.

Löwen, B., Schulz, S., 1995. Excess molar enthalpies of acetone + water, cyclohexane, methanol, 1-propanol, 2-propanol, 1-butanol and 1-pentanol at 283.15, 298.15,

323.15, 343.15 and 363.15 K. Thermochim Acta 262, 69–82. https://doi.org/10.1016/0040-6031(95)02278-A.

McKinnon, I., Williamson, A., 1964. Heats of mixing of carbon tetrachloride with dioxan. Aust. J. Chem. 17, 1374–1378. https://doi.org/10.1071/CH9641374.

Nagata, I., Tamura, K., Tokuriki, S., 1981. Excess enthalpies and complex formation of acetonitrile with acetone, chloroform, and benzene. Thermochim Acta 47, 315–331. https://doi.org/10.1016/0040-6031(81)80110-4.

Noles, J.R., Zollweg, J.A., 1992. Vapor-liquid equilibrium for chlorodifluoromethane + dimethyl ether from 283 to 395 K at pressures to 5.0 MPa. J. Chem. Eng. Data 37, 306–310. https://doi.org/10.1021/je00007a008.

Oh, B.C., Kim, Y., Shin, H.Y., Kim, H., 2004. Vapor–liquid equilibria for the system 1-propanol+n-hexane near the critical region. Fluid Phase Equilib. 220, 41–46. https://doi.org/10.1016/j.fluid.2003.12.004.

Ortega Saavedra, J., Plácido Suarez, J., 1996. Excess enthalpies of 12 binary liquid mixtures of .alpha.,.omega.-dibromoalkanes (C2, C4, C6) + benzene, toluene, ethylbenzene, or butylbenzene at 298.15 K. ELDATA Int. Electron. J. Phys.-Chem. Data ISSN 1265-1532 V 2 P 121-130.

Röck, H., Schröder, W. (Eds.), 1957. Dampf-Flüssigkeits-Gleichgewichtsmessungen im System Azeton-Chloroform. Z. Für Phys. Chem. 11, 41–55. https://doi.org/10.1524/zpch.1957.11.1.2.041.

Seo, J., Lee, S., Kim, H., 2002. Measurement and correlation of vapor–liquid equilibria for the ethanol + *n*-heptane system near the critical region. J. Chem. Eng. Data 47, 974–977. https://doi.org/10.1021/je0200057.

Shatas, J.P., Abbott, M.M., Van Ness, H.C., 1975. Excess thermodynamic functions for ternary systems. II. Chloroform-ethanol-*n*-heptane at 50Deg. J. Chem. Eng. Data 20, 406–409. https://doi.org/10.1021/je60067a020.

Sugi, H., Katayama, T., 1978. Ternary liquid-liquid and miscible binary vapor-liquid equilibrium data for the two systems n-hexane-ethanol-acetonitrile and water-acetonitrile-ethyl acetate. J. Chem. Eng. Jpn 11, 167–172. https://doi.org/10.1252/jcej.11.167.

Teodorescu, M., Barhala, A., Dragoescu, D., 2006. Isothermal (vapour + liquid) equilibria for the binary (cyclopentanone or cyclohexanone with 1,1,2,2-tetra-chloroethane) systems at temperatures of (343.15, 353.15, and 363.15) K. J. Chem. Thermodyn. 38, 1432–1437. https://doi.org/10.1016/j.jctc.2006.01.010.

Vinichenko, I.G., Susarev, M.P., 1966. Investigation and calculation of liquid-vapor equilibrium in the system acetone - ethanol - *n*-hexane. J. Appl. Chem. USSR 39, 1475–1478.

Carroll, B.H., Mathews, J.H., 1924. A Calorimeter for Heats of Mixing at Elevated Temperatures. J. Am. Chem. Soc. 46, 30–36. https://doi.org/10.1021/ja01666a004.

Ewert, M., 1936. Studies on the Theory of Concentrated Solutions. XIII. Aqueous Solutions of Organic Compounds. Bull. Soc. Chim. Belg. 45, 493–515.