Assessment of the SM12, SM8, and SMD Solvation Models for Predicting Limiting Activity Coefficients at 298.15 K

Abstract

The SMx (x = 12, 8, or D) universal solvent models are implicit solvent models which using electronic structure calculations can compute solvation free energies at 298.15 K. While solvation free energy is an important thermophysical property, within the thermodynamic modeling of phase equilibrium, limiting (or infinite dilution) activity coefficients are preferred since they may be used to parameterize excess Gibbs free energy models to model phase equilibrium. Conveniently, the two quantities are related. Therefore the present study was performed to assess the ability to use the SMx universal solvent models to predict limiting activity coefficients. Two methods of calculating the limiting activity coefficient where compared: (1) the solvation free energy and self-solvation free energy were both predicted and (2) the self-solvation free energy was computed using readily available vapor pressure data. Overall the first method is preferred as it results in a cancellation of errors, specifically for the case in which water is a solute. The SM12 model was compared to both the Universal Quasichemical Functional-group Activity Coefficients (UNIFAC) and modified separation of cohesive energy density (MOSCED) models. MOSCED was the highest performer, yet had the smallest available compound inventory. UNIFAC and SM12 exhibited comparable performance. Therefore further exploration and research should be conducted into the viability of using the SMx models for phase equilibrium calculations.

1. Introduction

Limiting (or infinite dilution) activity coefficients are of particular interest both in their fundamental and practical applications. As will be shown, the log limiting activity coefficient for solute i infinitely dilute in solvent j is proportional to the transfer free energy of i from a solution of pure i to pure j (in which it is infinitely dilute). Therefore fundamentally, the limiting activity coefficient sheds insight into how favorable the interactions of i with j are relative to i with itself. This fundamental molecular-level understanding provides the framework that directly impacts design schemes for chemical processes [1].

The practical applications are vast for they are used in the design of separation equipment, solvent selection, the design of pharmaceuticals, stripping operations, leaching operations, and the extraction of dilute materials [2,3]. Previous studies have also shown that limiting activity coefficients can be used to predict the formation of an azeotrope [4,5,6,7,8]. The equilibrium composition of a dilute mixture can also be determined, thus having industrial applications in azeotropic distillation and high-purity extraction [9]. This can serve as a means for process intensification in the environmental sphere, since environmental engineering addresses the implications of a chemical very dilute in a phase, often requiring dilute separation processes [10]. For example, limiting activity coefficients can be used to determine partition coefficients for a dilute species between two phases, which is critical in predicting the bio-accumulation of long-lived chemicals in the environment and food chain [2].

Of particular interest, it has been demonstrated that the limiting activity coefficient of each component in a binary mixture can be used to parameterize a binary interaction excess Gibbs free energy model to make phase equilibrium calculations over an entire range of concentrations [8,11,12,13]. These parameterizations are particularly advantageous for maximizing the efficiency of separation processes, which comprise a large majority of chemical plant operations and costs. The underlying principle is that limiting activity coefficients are a direct measure of the interactions in solution, and correspond to maximum deviations from ideality. It is therefore reasonable to assume we can extrapolate to more concentrated solutions where the solution would be more ideal (i.e., closer to the pure component limit). Consider for example the use of Wilson’s equation [13]. For a binary system we can compute the composition dependent activity coefficient of component 1 and 2, at mole fraction composition $x_1$ and $x_2=1-x_1$ as:

$$\begin{array}{cc}\hfill ln{\gamma }_1& =-ln\left(x_1+{\mathsf{\Lambda }}_{12}{x}_2\right)+x_2\left(\frac{{\mathsf{\Lambda }}_{12}}{x_1+{\mathsf{\Lambda }}_{12}{x}_2}-\frac{{\mathsf{\Lambda }}_{21}}{{\mathsf{\Lambda }}_{21}{x}_1+x_2}\right)\hfill \\ \hfill ln{\gamma }_2& =-ln\left(x_2+{\mathsf{\Lambda }}_{21}{x}_1\right)-x_1\left(\frac{{\mathsf{\Lambda }}_{12}}{x_1+{\mathsf{\Lambda }}_{12}{x}_2}-\frac{{\mathsf{\Lambda }}_{21}}{{\mathsf{\Lambda }}_{21}{x}_1+x_2}\right)\hfill \end{array}$$

(1)

where ${\mathsf{\Lambda }}_{12}$ and ${\mathsf{\Lambda }}_{21}$ are adjustable parameters which may be related to the binary (intermolecular) interaction parameters (BIPs) of the system ($a_{12}$ and $a_{21}$):

$$\begin{array}{c}\hfill {\mathsf{\Lambda }}_{12}=\frac{v_2}{v_1}exp\left[-\frac{a_{12}}{RT}\right]\\ \hfill {\mathsf{\Lambda }}_{21}=\frac{v_1}{v_2}exp\left[-\frac{a_{21}}{RT}\right]\end{array}$$

(2)

where R is the molar gas constant, T is the absolute temperature, and $v_1$ and $v_2$ are the molar volume of pure component 1 and 2, respectively. At infinite dilution, Equation (1) reduces to:

$$\begin{array}{cc}\hfill ln{\gamma }_1^{\infty }& =-ln\left({\mathsf{\Lambda }}_{12}\right)+1-{\mathsf{\Lambda }}_{21}\hfill \\ \hfill ln{\gamma }_2^{\infty }& =-ln\left({\mathsf{\Lambda }}_{21}\right)+1-{\mathsf{\Lambda }}_{12}\hfill \end{array}$$

(3)

It follows that with knowledge of the limiting activity coefficient of each component in the binary mixture, Equation (3) reduces to a system of two equations with two unknowns that can be used to parameterize Wilson’s equation. This would allow one to calculate composition dependent activity coefficients which could be used, for example, to model vapor/liquid equilibrium. Moreover, this approach may readily be generalized to multi-component systems. Wilson’s equation has two binary parameters for each pair in a multi-components system, where limiting activity coefficients may be used to parameterize each binary pair [13].

Limiting activity coefficients can be experimentally determined through a variety of methods including: gas chromatography, gas stripping, ebulliometry, and the dew point method [9,14]. However when there is little or no experimental data present, predictive models and tools can be implemented. Some popular methods include group contributions methods, solubility parameter based methods, molecular simulations, and electronic structure calculations. The Universal Quasichemical Functional-group Activity Coefficients (UNIFAC) method is a group contribution-based binary interaction excess Gibbs free energy model designed by breaking molecules into functional groups and calculating the configuration energy of the molecule as the sum of group interaction energies [15,16]. Similarly, mod-UNIFAC (Dortmund) is a revision of the original UNIFAC model that accounts for the temperature dependence of activity coefficients by including temperature-dependent interaction parameters within the model [17,18]. Solubility parameter based methods such as the modified separation of cohesive energy density (MOSCED) [8,19,20] and Hansen solubility parameter (HSP) methods [21], while they limit the user to modeling systems for which solute parameters are available, they give molecular-level insight into the the behavior of a binary system. Molecular simulations and electronic structure calculations can be computationally rigorous, yet are advantageous for they only require the input of a molecular structure [22,23,24,25,26].

The pursuit of accurate, efficient predictive models is vital as it could reduce the experimental and design cost associated with chemical processes. This is of particular interest in separation processes, as they may account for 40–70 percent of the capital and operating cost of a chemical plant [27]. Predictive computational tools have become an essential asset in conducting feasibility and comparative studies, as they are significantly quicker and more cost-efficient than experimental determinations [28]. They are necessary when property data is unavailable, specifically in the case of dealing with new products and compounds. However within computational advancement, there is a risk in losing the intuitive molecular understanding of the processes that lead to observed outcomes. Thus, including thermodynamic parameters such as limiting activity coefficients in the computational process provides molecular insight that could greatly assist in resultant recommendations for industrial design [10].

This study seeks to assess the capability of the SMx (x = 12, 8, or D) universal solvent models for predicting limiting activity coefficients [29,30,31]. The SMx universal solvent models are implicit or continuum solvent models that can compute solvation free energies at 298.15 K using electronic structure calculations. The models have long been under active development with their accuracy continuously improving and are implemented in a number of software packages [29,32,33,34]. The only input needed from the user is a structure of the molecule; therefore the SMx models are set apart from other methods in the sheer simplicity for the user. Of specific interest to us is their potential future use to parameterize binary interaction excess Gibbs free energy models for phase equilibrium predictions. In addition, a major goal of this work is to determine how best to predict limiting activity coefficients using solvation free energies computed using the SMx universal solvent models. Specifically, as will be shown, the prediction of the limiting activity coefficient is dependent on both the solvation free energy and self-solvation free energy. While the self-solvation free energy may be computed directly, it may additionally be computed from readily available vapor pressure data.

Previous work to use the SMx universal solvent models to predict limiting activity coefficients is limited. Recently, Lisboa and Pliego [26] used SMD to compute 46 limiting activity coefficients. In the present study, we compute over 3000 limiting activity coefficients with comparison to reference data, for each SM12, SM8, and SMD. In addition, we compare two methods of calculating the limiting activity coefficient: (1) the solvation free energy and self-solvation free energy were both predicted and (2) the self-solvation free energy was computed using readily available vapor pressure data. To better understand the two methods, we compare the individual contributions of the solvation and self-solvation free energy to available reference data. In general, there is not a significant difference between the limiting activity coefficients computed using the two methods. However, for the cases in which water is a solute, the results are significantly better when both the solvation and self-solvation free energy are predicted using SMx. This results from a cancellation of errors in the two terms. We therefore recommend that in general SMx be used to predict both the solvation and self-solvation free energy when computing limiting activity coefficients. Comparison is additionally made to predictions using UNIFAC and MOSCED. SM12 and UNIFAC are found to exhibit comparable performance, motivating future exploration and research into the viability of using the SMx models for phase equilibrium calculations.

2. Method

Performing electronic structure calculations in the SMx (x = 12, 8, or D) universal solvent model, one calculates the solvation free energy of a solute (component i) at infinite dilution in a solvent (component j), $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$, where i = {1 or 2} and $j=$ {1 or 2} [29,30,31,35,36]. The solvation free energy in this context is defined as taking a solute from a non-interacting ideal gas state to solution at the same molecular density (or concentration). We have shown previously that the solvation free energy is readily related to the limiting activity coefficient (${\gamma }_{i,j}^{\infty }$) as [24,37,38,39,40,41,42]:

$$ln{\gamma }_{i,j}^{\infty }\left(T,P\right)=\frac{1}{RT}\left[\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}\left(T,P\right)-\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}\left(T,P\right)\right]+ln\frac{v_i\left(T,P\right)}{v_j\left(T,P\right)}$$

(4)

where $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ is the solute “self”-solvation free energy, R is the molar gas constant, and $v_i$ and $v_j$ correspond to the pure liquid molar volume of component i and j, respectively. The difference $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}-\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ is equivalent to the transfer free energy of i from a solution of pure i to pure j (in which it is infinitely dilute). When $i=j$ we obtain the correct limiting behavior that ${\gamma }_{i,j}^{\infty }=1$. The self-solvation free energy, $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$, may readily be calculated by performing a solvation free energy calculation for component i in itself. Alternatively, we can relate $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ to the pure liquid fugacity of component i, $f_i^0$, as [38]:

$$ln{f}_i^0\left(T,P\right)=\frac{1}{RT}\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}\left(T,P\right)+ln\frac{RT}{v_i\left(T,P\right)}$$

(5)

We can expand $f_i^0$ as [13]:

$$f_i^0\left(T,P\right)={\varphi }_i^{\mathrm{sat}}\left(T\right)P_i^{\mathrm{sat}}\left(T\right)exp\left\{{\int }_{P^{\mathrm{sat}}}^P\frac{v_i\left(T,P\right)}{RT}dP\right\}$$

(6)

where ${\varphi }_i^{\mathrm{sat}}$ and $P_i^{\mathrm{sat}}$ are the fugacity coefficient and vapor pressure of pure component i at saturation at T, and the term in brackets is the Poynting correction, and accounts for the change in fugacity in going from $P^{\mathrm{sat}}$ to P at constant T. If we assume that the vapor phase is an ideal gas and that the Poynting correction is negligible, we have [13,38,43,44]:

$$f_i^0\left(T,P\right)=P_i^{\mathrm{sat}}\left(T\right)$$

(7)

If we are at low pressures well removed from the critical point, and we have a non-self associating fluid (i.e., no carboxylic acid), use of this expression is generally reasonable. Often times, the pure component property $P_i^{\mathrm{sat}}$ is known. This presents an opportunity to use reference data to compute $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ as:

$$\frac{1}{RT}\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}\left(T,P\right)=ln\frac{P_i^{\mathrm{sat}}\left(T\right)v_i\left(T,P\right)}{RT}$$

(8)

and then using electronic structure calculations in the SMx universal solvent models to compute $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$, which is a property of the mixture.

3. Computational Details

All solvation free energy calculations were performed in the SMx (x = 12, 8, or D) universal solvent models using QChem 5.1.2 (supplementary materials) [45,46]. We selected systems by searching through the Minnesota Solvation Database version 2012 [47] for binary systems for which QChem had (SMx) solvent parameters for both components. The Minnesota Solvation Database was selected as our reference for comparison as it was used in the parameterization of SM12 [29], and earlier versions were used in the parameterization of SM8 and SMD [30,31]. The data in the Minnesota Solvation Database is restricted to 298.15 K, and it follows that calculations using the SMx solvent models are restricted to 298.15 K. In addition, the Minnesota Solvation Database contains only solvation and transfer free energy data. We computed reference limiting activity coefficients from the solvation free energies using Equation (4), where the self-solvation free energy was computed via Equation (8) using reference vapor pressure and molar volume data [48]. Note that iodine and carboxylic acid containing molecules were excluded from our search. We excluded the seven iodine hydrocarbons because SM8 and SMD were not parameterized to treat iodine containing solutes [30,31]. The six molecules containing a carboxylic acid functional group were excluded as they are known to dimerize in the vapor phase [49], and we wish to compare the use of Equations (5) and (7) when computing limiting activity coefficients. Excluding these 13 molecules, QChem contains SMx parameters for 148 additional molecules (solvents).

This resulted in a data set of solvation free energies (and hence limiting activities) for 1140 systems, consisting of 107 unique solutes and 73 unique solvents. The Minnesota Solvation Database contained self-solvation free energies at 298.15 K for only 62 of the 148 possible molecules. We therefore instead computed a set of reference self-solvation free energies for all 148 molecules using Equation (8) with reference vapor pressure and molar volume data [48]. For the 62 molecules for which values were available in the Minnesota Solvation Database, the computed values were found to be consistent. The inclusion of both solvation free energies and self-solvation free energies will allow us to look at the individual contributions to the limiting activity coefficient.

As already mentioned, the Minnesota Solvation Database was used in the parameterization of the SMx universal solvent models. Its use was motivated by the desire to assess the prediction of the individual contributions, $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$ and $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$, and the combined effect on $ln{\gamma }_{i,j}^{\infty }$. As a secondary test, we curated reference data from Parts 1 to 6 of DECHEMA’s “Activity Coefficients at Infinite Dilution” volumes [50,51,52,53,54,55]. We again restricted ourselves to binary systems for which both components were in our set of 148 molecules (solvents) for which QChem contains SMx parameters. This resulted in limiting activity coefficients for an additional 2163 systems, consisting of 117 unique solutes and 77 unique solvents. Note that as compared to the Minnesota Solvation Database wherein there was only a single entry for each unique system (i.e., solute-solvent pair), duplicates were present in DECHEMA. We retained all of the data, including duplicates, in the final set.

When curating reference data, we only included binary systems for which QChem had (SMx) solvent parameters for both components. The reason for this is two-fold. (1) First, if one were to compute both $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$ and $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ to predict $ln{\gamma }_{i,j}^{\infty }$ via Equation (4), a solvation free energy would be required when both i and j take on the role of solvent, requiring that SMx solvent parameters be available for both components. (2) On the other hand, if $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ was instead computed using reference vapor pressure data, then in theory solvent parameters for i would not be necessary. However, even if this is the case, we are interested in the potential future application of using SMx computed solvation free energies to parameterize a binary interaction excess Gibbs free energy model for phase equilibrium calculations. For this case, for a binary pair both ${\gamma }_{i,j}^{\infty }$ and ${\gamma }_{j,i}^{\infty }$ are required; no matter how the self-solvation free energy is computed, a solvation free energy calculation is required in both solvents, making necessary the availability of SMx solvent parameters for both components.

As we will discuss, we find that challenges occur in predicting limiting activity coefficients for aqueous systems. We therefore sourced an additional data set composed of organics in water and water in organics at 298.15 K from “Yaws’ Handbook of Properties for Aqueous Systems” [56]. Data was available for 121 systems where water is the solvent and 90 systems where water is the solute, with no duplicates.

The calculations proceeded by first obtaining a daylight SMILES [57,58] for each of the 148 molecules from PubChem [59]. We then generated 3-D structures from the SMILES using Open Babel 2.3.2 [60,61]. Open Babel was then used to perform a systematic conformational search to identify the lowest energy conformer, followed by geometry optimization. The conformational search and geometry optimization were conducted using the General Amber Force Field [62,63] with Gasteiger partial charges [64]. The final structures were then further refined by performing a geometry optimization in vacuum using QChem 5.1.2 at the M06-2X/cc-pVDZ level of theory/basis set [35,65,66]. This final structure was then used for all of the solvation free energy calculations.

To calculate the solvation (including self-solvation) free energies, single point energy calculations were performed in vacuum and in the SM12, SM8, and SMD universal solvent models for each solvent of interest at the M06-2X/6-31G(d) level of theory/basis set using QChem 5.1.2. Using QChem 5.1.2, the calculation of the solvation free energy is handled completely internally by the software. We emphasize that we used the same conformation for the solute in the ideal gas (vacuum) and solution phases, consistent with the SM12, SM8, and SMD parameterizations [29,30,31]. It is possible to consider different configurations for the two phases, and it is possible to incorporate multiple configurations into the calculation of the solvation free energy which may be important for some systems [67]. However, this was not pursued here for simplicity.

All of the SMx computed self-solvation free energies ($\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$), solvation free energies ($\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$), and limiting activity coefficients (${\gamma }_{i,j}^{\infty }$), along with the reference values and results of the reference calculations (UNIFAC and MOSCED) to which comparison is made, are tabulated in the supporting information accompanying the electronic version of this manuscript. As mentioned, we are interested in the potential future application of using SMx computed solvation free energies to parameterize a binary interaction excess Gibbs free energy model for phase equilibrium calculations. Therefore, in the supporting information the results are additionally categorized following the system classification used in DECHEMA’s “Vapor-Liquid Equilibrium Data Collection” volumes [68].

4. Results and Discussion

The central goal of this work is to assess the ability of electronic structure calculations with the SMx (x = 12, 8, or D) universal solvent models to predict limiting activity coefficients. As seen in Equation (4), this is composed of contributions due to the solvation free energy ($\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$) and self solvation free energy ($\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$). We will therefore first assess each individual contribution. We will quantify the agreement between the reference data and the SMx predictions using the root mean square error (RMSE):

$$\mathit{RMSE}=\sqrt{\frac{1}{N}\sum _{i=1}^N{\left(Y^{\mathrm{pred}}-Y^{\mathrm{ref}}\right)}^2}$$

(9)

the average absolute percent deviation (AAPD):

$$\mathit{AAPD}=\frac{100}{N}\sum _{i=1}^N\left|\frac{Y^{\mathrm{pred}}-Y^{\mathrm{ref}}}{Y^{\mathrm{ref}}}\right|$$

(10)

and the mean unsigned error (MUE):

$$\mathit{MUE}=\frac{1}{N}\sum _{i=1}^N\left|Y^{\mathrm{pred}}-Y^{\mathrm{ref}}\right|$$

(11)

where the superscript “pred” and “ref” correspond to the predicted and reference values, respectively, Y corresponds to the property of interest, and N corresponds to the number of systems for which comparison is being made.

4.1. Self-Solvation Free Energy

We begin with $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$. A parity plot of the predicted versus reference values of $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}/\left(RT\right)$ and a summary of the computed statistics is provided in Figure 1 and

Table 1. A summary of the squared correlation coefficient ($R^2$), root mean square error (RMSE), average absolute percent deviation (AAPD), and mean unsigned error (MUE) for the values of $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}/\left(RT\right)$ predicted using the SM12, SM8, and SMD universal solvent model as compared to the reference data.

	$\mathsf{\Delta }{\mathit{G}}_{\mathit{i},\mathit{i}}^{\mathbf{self}}/\left(\mathit{R}\mathit{T}\right)$
	SM12	SM8	SMD
$R^2$	0.713	0.707	0.566
RMSE	1.526	1.600	1.820
AAPD	10.888	11.230	14.655
MUE	1.077	1.131	1.372

, respectively. The statistics in Table 1 are computed for $Y=\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}/\left(RT\right)$, with the summation carried out over all $N=148$ molecules for which reference data is available. Note that we are interested in the dimensionless $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ as it is the quantity of interest in Equation (4).

Overall, we find that SM12 and SM8 outperform SMD, with SM12 performing slightly better than SM8. The top performance of SM12 is in agreement with the recent SM12 parameterization [29]. Visually, it is clear in Figure 1 that SM12 and SM8 offer a noticeably smaller spread in the data as compared to SMD, corresponding to a smaller overall error.

During the parameterization of SM12, it was pointed out that the largest deviation was observed for the self-solvation free energy of water, which was overly negative by 1 to 2 kcal/mol (or 1.7 to 3.4 $RT$ in dimensionless units) using all methods, and the solvation free energy of the water dimer in water was overly negative by 3 to 4 kcal/mol (or 5.1 to 6.8 $RT$ in dimensionless units) using all methods [29]. We likewise observe large errors for the self-solvation free energy of water, although our predictions are slightly more negative. The reference value of $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}/\left(RT\right)$ for water is −10.689, while we predict −14.473, −15.485, and −15.024 with SM12, SM8, and SMD, respectively. This would correspond to predicting a vapor pressure that is too small, and that we are over estimating the affinity of water for the liquid phase or equivalently the strength of the water–water interactions.

4.2. Solvation Free Energy

Next we consider $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$, making comparison for the 1140 systems, consisting of 107 unique solutes and 73 unique solvents, from the Minnesota Solvation Database. A parity plot of the predicted versus reference values of $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}/\left(RT\right)$ and a summary of the computed statistics is provided in Figure 2 and

Table 2. A summary of the squared correlation coefficient ($R^2$), root mean square error (RMSE), average absolute percent deviation (AAPD), and mean unsigned error (MUE) for the values of $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}/\left(RT\right)$ predicted using the SM12, SM8, and SMD universal solvent model as compared to the reference data from the Minnesota Solvation Database.

	$\mathsf{\Delta }{\mathit{G}}_{\mathit{i},\mathit{j}}^{\mathbf{solv}}/\left(\mathit{R}\mathit{T}\right)$
	SM12	SM8	SMD
$R^2$	0.838	0.826	0.820
RMSE	1.120	1.125	1.181
AAPD	20.985	24.593	26.093
MUE	0.828	0.843	0.882

. The statistics in Table 2 are computed for $Y=\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}/\left(RT\right)$, with the summation carried out over all $N=1140$ systems for which reference data is available.

We again find that SM12 is the top performing model, although the difference in performance between the models is much smaller for $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$ as compared to $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$. Moreover, we find that except for the AAPD, the overall errors are all improved for $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$ relative to $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$. Comparing Figure 1 and Figure 2, there appears to be a smaller spread in the $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$ data, which would correspond to the smaller RMSE and MUE. This at first had us surprised. In theory, $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ corresponds to a single component, homogeneous system. On the other hand, $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$ corresponds to a single solute infinitely dilute in a solvent. In this case we need to account for solute-solvent and solvent-solvent interactions which are not the same. We would therefore assume that the latter case is more difficult to predict and would result in a lower accuracy. However, this was not the case. We do not offer an explanation for this difference, but would like to point out two differences in the reference data sets. (1) Our data set for $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$ consists of 1140 systems, while we have just 148 values of $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$. Of the 148 values of $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$, only 62 were included in the Minnesota Solvation Database. (2) Comparing Figure 1 and Figure 2, we see there is a difference in the range of values in each data set. While a majority of the values of $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}/\left(RT\right)$ are over the range −15 to −5, there is an appreciable number of values of $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}/\left(RT\right)$ greater than −5. We also point out that this last point contributes to the larger AAPD for $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$. The computed errors and top performance of SM12 are consistent with those published in the recent SM12 parameterization [29].

Considering again the original SM12 parameterization, it was found that the errors for the solvation free energy (both aqueous and non-aqueous) for the solute class containing molecular hydrogen, ammonia, water, and the water dimer were relatively large compared to the other solute classes [29]. As already mentioned, the original SM12 parameterization pointed out a particular challenge with the self-solvation free energy of water which was predicted to be overly negative. To explore this further, we broke down $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$ to cases where water was the solvent ($j=$ water), water was the solute ($i=$ water), and when neither the solvent or solute was water. Note that $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ was not include in the set of $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$. This resulted in 104 systems wherein water was the solvent, 14 systems wherein water was the solute, and 1022 systems wherein water was neither the solvent or solute. A parity plot of the predicted versus reference values of $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}/\left(RT\right)$ with SM12 wherein these three cases are differentiated is provided in Figure 3.

We find that in general when water is the solute, $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$ is predicted to be overly negative. We find that $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}/\left(RT\right)$ is overly negative on average by 3.4, 3.5, and 2.4 (dimensionless units) when using SM12, SM8, and SMD, respectively. Compare this to the case of $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}/\left(RT\right)$ for water predicted to be overly negative by 3.8, 4.8, and 4.3 (dimensionless units) for SM12, SM8, and SMD, respectively.

4.3. Limiting Activity Coefficients

When predicting limiting activity coefficients (${\gamma }_{i,j}^{\infty }$) with the SMx (x = 12, 8, or D) universal solvent models, two different strategies were evaluated. (1) We predicted both $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ and $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$, then using reference values of the molar volume [48], ${\gamma }_{i,j}^{\infty }$ was computed via Equation (4). (2) In the second approach, since the pure component vapor pressure ($P_i^{\mathrm{sat}}$) is commonly available [48], we used this reference data to compute $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ via Equation (8). We then combined this reference value of $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ with our SMx predicted value of $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$ to predict ${\gamma }_{i,j}^{\infty }$. When computing statistics, for AAPD we use $Y={\gamma }_{i,j}^{\infty }$. For all other cases ($R^2$, RMSE, and MUE) we use $Y=ln{\gamma }_{i,j}^{\infty }$.

4.3.1. Minnesota Solvation Database

Limiting activity coefficients (${\gamma }_{i,j}^{\infty }$) were computed for the same 1140 systems from the Minnesota Solvation Database for which $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$ was computed. Reference values of ${\gamma }_{i,j}^{\infty }$ were computed via Equation (4) using the reference values of $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$ from the Minnesota Solvation Database, the reference values of $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ already described, and reference values of the molar volume [48]. A parity plot of the predicted versus reference values of $ln{\gamma }_{i,j}^{\infty }$ is provided in Figure 4. In the top pane results are shown when we predict both $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ and $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$, and the bottom pane shows the case when $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ was computed using reference data.

In general, there does not appear to be a significant difference between the predictions wherein $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ was predicted or computed using reference data. The spread in results does appear to decrease when $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ is computed using reference data, which may reflect the spread in the predicted values of $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ (see Figure 1). A comparison of the statistics is provided in

Table 3. A summary of the squared correlation coefficient ($R^2$), root mean square error (RMSE), and mean unsigned error (MUE) for $ln{\gamma }_{i,j}^{\infty }$, and the average absolute percent deviation (AAPD) for ${\gamma }_{i,j}^{\infty }$ predicted using the SM12, SM8, and SMD universal solvent model as compared to the reference data from the Minnesota Solvation Database. For the predictions, we compare the case where (1) $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ is predicted using the SM12, SM8, and SMD universal solvent model, and where (2) $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ is computed using reference data, as indicated.

	$\mathsf{\Delta }{\mathit{G}}_{\mathit{i},\mathit{i}}^{\mathbf{self}}$ Predicted			$\mathsf{\Delta }{\mathit{G}}_{\mathit{i},\mathit{i}}^{\mathbf{self}}$ Reference
	SM12	SM8	SMD	SM12	SM8	SMD
$R^2$	0.796	0.740	0.627	0.806	0.796	0.788
RMSE	1.155	1.288	1.801	1.120	1.125	1.181
AAPD	89.856	107.821	357.706	92.482	134.797	131.626
MUE	0.840	0.883	1.376	0.828	0.842	0.882

.

Again, while the difference in the two sets of predictions does not appear to be large, in all cases we find that $R^2$ increases and the computed overall errors decrease when using reference values of $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$. One might expect this result. Consider Equation (4) and the errors summarized in Table 1, Table 2 and Table 3, where we will take RMSE as a measure of the error in the predictions. When using reference values of $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$, the only predicted quantity in Equation (4) is $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$. As a result, we find for this case that the RMSE for the predicted values of $ln{\gamma }_{i,j}^{\infty }$ is identical to that for $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}/\left(RT\right)$. Next, consider the case of when $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ is predicted. Using propagation of error, we would expect the error (taken to be RMSE) in $ln{\gamma }_{i,j}^{\infty }$ to be equal to $\sqrt{\mathit{RMS}{E}_{\mathrm{solv}}^2+\mathit{RMS}{E}_{\mathrm{self}}^2}$, where $\mathit{RMS}{E}_{\mathrm{solv}}$ and $\mathit{RMS}{E}_{\mathrm{self}}$ are the RMSE for the predicted values of $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}/\left(RT\right)$ and $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}/\left(RT\right)$, respectively. This would result in an RMSE of 1.893, 1.970, and 2.170 for SM12, SM8, and SMD, respectively. However, we find that the actual error is less than the predicted error. The expression used assumes that the error is random [69]. The smaller observed uncertainty suggests the presence of systematic error.

With this in mind, let us again look at the specific case of water. We found earlier that when water is the solute, on average $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}/\left(RT\right)$ is predicted to be overly negative by 3.4, 3.5, and 2.4 (dimensionless units) when using SM12, SM8, and SMD, respectively. Further, $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}/\left(RT\right)$ for water is predicted to be overly negative by 3.8, 4.8, and 4.3 (dimensionless units) for SM12, SM8, and SMD, respectively. From Equation (4), this leads to a convenient cancellation of errors in the calculation of $ln{\gamma }_{i,j}^{\infty }$. For the case of SM12, using the predicted value of $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$, this would suggest that on average the value of $ln{\gamma }_{i,j}^{\infty }$ is overly positive by just 0.4. On the other hand, using the reference value of $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$, the corresponding value would be overly negative by 3.4. This is illustrated for SM12 in Figure 5 where the predictions were again broken down to cases where water was the solvent ($j=$ water), water was the solute ($i=$ water), and when neither the solvent or solute was water.

Again, in general we find that the predicted values of $ln{\gamma }_{i,j}^{\infty }$ do not change substantially when $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ was predicted or computed using reference vapor pressure data. However, when water is the solute, there is a clear difference. This would suggest that when water is the solute of interest, both $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$ and $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ should be computed consistently in the calculation of $ln{\gamma }_{i,j}^{\infty }$.

In Figure 5 we additionally point out that we observe an apparent break in our trend for reference values of $ln{\gamma }_{i,j}^{\infty }>9$. For all cases when the reference values of $ln{\gamma }_{i,j}^{\infty }$ are greater than 9, we under-predict the value of $ln{\gamma }_{i,j}^{\infty }$. These would correspond to highly non-ideal (unfavorable) systems. For all 26 cases, water is the solvent. A possible source of error could be due to the difficulty of experimentally measuring such non-ideal values [70].

4.3.2. DECHEMA

Next, as a secondary test, values of ${\gamma }_{i,j}^{\infty }$ at 298.15 K were computed for 2163 systems, consisting of 117 unique solutes and 77 unique solvents, from Parts 1 to 6 of DECHEMA’s “Activity Coefficients at Infinite Dilution” volumes [50,51,52,53,54,55]. A parity plot of the predicted versus reference values of $ln{\gamma }_{i,j}^{\infty }$ is provided in Figure 6. In the top pane of Figure 6 results are shown when we predict both $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ and $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$, and the bottom pane shows the case when $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ was computed using reference data. A summary of the statistics is provided in

Table 4. A summary of the squared correlation coefficient ($R^2$), root mean square error (RMSE), and mean unsigned error (MUE) for $ln{\gamma }_{i,j}^{\infty }$, and the average absolute percent deviation (AAPD) for ${\gamma }_{i,j}^{\infty }$ predicted using the SM12, SM8, and SMD universal solvent model as compared to the reference data from DECHEMA. For the predictions, we compare the case where (1) $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ is predicted using the SM12, SM8, and SMD universal solvent model, and where (2) $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ is computed using reference data, as indicated.

	$\mathsf{\Delta }{\mathit{G}}_{\mathit{i},\mathit{i}}^{\mathbf{self}}$ Predicted			$\mathsf{\Delta }{\mathit{G}}_{\mathit{i},\mathit{i}}^{\mathbf{self}}$ Reference
	SM12	SM8	SMD	SM12	SM8	SMD
$R^2$	0.911	0.913	0.874	0.865	0.889	0.886
RMSE	1.163	1.065	1.382	1.354	1.196	1.271
AAPD	86.823	120.158	308.173	177.079	166.345	289.780
MUE	0.847	0.768	1.012	1.061	0.875	0.910

.

We again breakdown the predictions to cases where water was the solvent ($j=$ water), water was the solute ($i=$ water), and when neither the solvent or solute was water, and plot the results for SM12 in Figure 7. The breakdown resulted in 437 systems where water is the solvent, 52 systems where water is the solute, and 1668 systems where water was neither the solvent or solute.

In general, our observations are the same as when we compared to the Minnesota Solvation Database. First, there again does not appear to be a significant difference between the predictions wherein $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ was predicted or computed using reference data. However, for the specific case of water as a solute, there is a noticeable difference. When $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ was computed using reference vapor pressure data, $ln{\gamma }_{i,j}^{\infty }$ when water was a solute was under-predicted on average by 3.02, 3.68, and 2.84 with SM12, SM8, and SMD, respectively. For all of theses cases, the value of $ln{\gamma }_{i,j}^{\infty }$ was always under-predicted except for water in tetrachloroethene with SMD. Additionally, we again notice an abrupt deviation in our trend for reference values of $ln{\gamma }_{i,j}^{\infty }$ greater than 9. As before, all of these highly non-ideal (unfavorable) systems, water was the solvent. Interestingly, the value of $R^2$ does increase. However, we are dealing with a larger data set that does contain duplicates.

4.3.3. Yaws’ Aqueous Systems

Having identified challenges in predicting $ln{\gamma }_{i,j}^{\infty }$ for aqueous systems, we sourced an additional data set composed of organics in water and water in organics at 298.15 K from “Yaws’ Handbook of Properties for Aqueous Systems” [56]. Data was available for 121 systems where water is the solvent and 90 systems where water is the solute, with no duplicates. A parity plot of the results for SM12 is provided in Figure 8. In the top pane of Figure 8 results are shown when we predict both $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ and $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$, and the bottom pane shows the case when $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ was computed using reference data. A summary of the statistics is provided in

Table 5. A summary of the squared correlation coefficient ($R^2$), root mean square error (RMSE), and mean unsigned error (MUE) for $ln{\gamma }_{i,j}^{\infty }$, and the average absolute percent deviation (AAPD) for ${\gamma }_{i,j}^{\infty }$ predicted using the SM12, SM8, and SMD universal solvent model as compared to the reference data from “Yaws’ Handbook of Properties for Aqueous Systems” [56]. For the predictions, we compare the case where (1) $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ is predicted using the SM12, SM8, and SMD universal solvent model, and where (2) $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ is computed using reference data, as indicated. The results are additionally separated into case where water was the solvent and where water was the solute, as indicated.

	$\mathsf{\Delta }{\mathit{G}}_{\mathit{i},\mathit{i}}^{\mathbf{self}}$ Predicted			$\mathsf{\Delta }{\mathit{G}}_{\mathit{i},\mathit{i}}^{\mathbf{self}}$ Reference
	SM12	SM8	SMD	SM12	SM8	SMD
	water as solute			water as solute
$R^2$	0.889	0.793	0.645	0.889	0.793	0.645
RMSE	1.339	1.836	2.482	3.721	3.600	3.175
AAPD	168.491	645.592	2598.544	94.691	93.947	97.359
MUE	1.116	1.625	2.168	3.485	3.397	2.733
	water as solvent			water as solvent
$R^2$	0.913	0.928	0.907	0.916	0.946	0.951
RMSE	1.854	1.658	1.652	1.464	1.300	1.090
AAPD	117.058	160.651	327.313	409.655	197.774	176.344
MUE	1.353	1.191	1.317	1.115	1.017	0.792

, broken down into the case where water is the solvent and water is the solute.

The results here are consistent with the other data sets. Considering first $ln{\gamma }_{i,j}^{\infty }$ for the case of water as a solute when $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ was computed using reference vapor pressure data. For all cases with SM12 and SM8 $ln{\gamma }_{i,j}^{\infty }$ is under predicted with an average error of 3.485 and 3.397, respectively. With SMD, $ln{\gamma }_{i,j}^{\infty }$ is under predicted for 85 of the 90 systems, with an average error of 2.834. On the other hand, when $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ was also predicted, with SM12, SM8, and SMD the number of systems for which $ln{\gamma }_{i,j}^{\infty }$ was under predicted was 22, 13, and 17, respectively. Most importantly, with SM12, SM8, and SMD, the overall MUE noticeably decreased.

When water is the solvent, we find that with SM12, SM8, and SMD, the value of $R^2$ increases, and the value of the MUE and RMSE decreases for $ln{\gamma }_{i,j}^{\infty }$. The difference in behavior results from a systematic error when water is the solute, where when $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ for water is also predicted, we obtain a convenient cancellation of errors. We note that we find the values of AAPD to be less informative as the value can be dominated by just a few systems.

4.3.4. Free SMx ${\gamma }_{i,j}^{\infty }$ Calculator

To facilitate the use of the SM12, SM8, and SMD solvation models to predict ${\gamma }_{i,j}^{\infty }$ and ultimately phase equilibrium, we are actively developing HTML (HyperText Markup Language) based software that can run in a web browser. The software is currently housed on GitHub at: https://github.com/SpencerSabatino/SolvationFreeEnergyWebsite. For the 148 molecules for which QChem 5.1.2 [45,46] has SMx parameters (which excludes iodine and carboxylic acid containing molecules), the user can select any of the 148 molecules to form a binary pair and compute $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$, $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$, and ${\gamma }_{i,j}^{\infty }$ for each component. The QChem 5.1.2 calculation results used by the software were generated in the present study. The software will continue to be developed with additional features added.

4.3.5. Reference Calculations

Overall, we found that the SMx (x = 12, 8, or D) predicted values of $ln{\gamma }_{i,j}^{\infty }$ did not change substantially when $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ was predicted compared to when it was computed using reference vapor pressure data, except in the case when water is the solute. Therefore, moving forward the selected method for calculating $ln{\gamma }_{i,j}^{\infty }$ will be to compute $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$ and $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ consistently; put differently, both $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$ and $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ will be predicted.

To provide further insight into the potential utility of electronic structure calculations with SMx universal solvent models to predict $ln{\gamma }_{i,j}^{\infty }$, the top performing SM12 was compared to two well-known methods, the Universal Quasichemical Functional-group Activity Coefficients (UNIFAC) method [15,16,71] and the Modified Separation of Cohesive Energy Density (MOSCED) model [19,20,72]. UNIFAC was selected as a comparative method because it is perhaps the most common predictive excess Gibbs free energy method [73]. Additionally, it is widely used and integrated into most process simulators. Brouwer and Schuur [73] recently conducted a comparative study of eight predictive methods to compute $ln{\gamma }_{i,j}^{\infty }$ at 298.15 K. Of the methods evaluated, the solubility parameter method MOSCED was the top performer. This strong performance of MOSCED is consistent with the findings of the 2005 re-parameterization [20,74] and our recent work [8,75]. MOSCED was therefore selected as an additional comparison alongside UNIFAC. Note that for water, we adopted the revised MOSCED parameters from our recent work [75].

The Minnesota Solvation Database, DECHEMA, and Yaws’ reference databases were constrained to match the compound availability provided by the two models, UNIFAC and MOSCED. In comparison to the Minnesota Solvation Database, UNIFAC was able to make predictions for 1065 (of 1140) systems with 99 unique solutes and 68 unique solvents. MOSCED was able to make predictions for 732 (of 1140) systems from the Minnesota Solvation Database with 58 unique solutes and 54 unique solvents. In regards to DECHEMA, UNIFAC was able to model 1988 (of 2163) systems with 112 unique solutes and 70 unique solvents. MOSCED was able to model 1818 (of 2163) systems with 71 unique solutes and 59 unique solvents. We again note that the DECHEMA data set includes repeated systems. All of the UNIFAC calculations were performed with CHEMCAD 7 [76]. A simple GNU Octave/MATLAB [77,78] M-file to predict limiting activity coefficients using MOSCED, a screen cast tutorial and additional references are available on the YouTube Channel of Andrew Paluch [79].

A parity plot of the results for SM12, UNIFAC and MOSCED is provided in Figure 9. In the top pane of Figure 9 results for the Minnesota Solvation Database are displayed, and in the bottom pane the corresponding DECHEMA results are displayed. A summary of the statistics is included in

Table 6. A summary of the squared correlation coefficient ($R^2$), root mean square error (RMSE), and mean unsigned error (MUE) for $ln{\gamma }_{i,j}^{\infty }$, and the average absolute percent deviation (AAPD) for ${\gamma }_{i,j}^{\infty }$ predicted using UNIFAC and MOSCED as compared to the reference data from the Minnesota Solvation Database version 2012 [47] and DECHEMA [50,51,52,53,54,55].

	Minnesota		DECHEMA
	UNIFAC	MOSCED	UNIFAC	MOSCED
$R^2$	0.820	0.884	0.915	0.976
RMSE	1.077	0.802	1.157	0.524
AAPD	52.595	41.251	41.044	131.758
MUE	0.634	0.422	0.559	0.241

.

For both the Minnesota Solvation Database and DECHEMA, MOSCED achieves the largest $R^2$ and smallest RMSE and MUE for $ln{\gamma }_{i,j}^{\infty }$. For DECHEMA, MOSCED does have a noticeably larger AAPD for ${\gamma }_{i,j}^{\infty }$, but we again note that the value AAPD can be heavily influenced by a few systems. Comparing to SM12 (see Table 3 and Table 4), UNIFAC does slightly outperform SM12. However, considering $R^2$ and RMSE, the UNIFAC and SM12 results are comparable, especially with regards to the DECHEMA data set.

In following the same methodology conducted previously, we sought to analyze the models’ ability to predict limiting activity coefficients when water is the solute and when water is the solvent. The reference data from “Yaws’ Handbook of Properties for Aqueous Systems” again served as the basis for systems that were modeled using the SM12 universal continuum solvent model, UNIFAC and MOSCED. A parity plot is provided in Figure 10 where the top pane contains data where water is the solute and the bottom pane where water is the solvent. The corresponding statistics are recorded in

Table 7. A summary of the squared correlation coefficient ($R^2$), root mean square error (RMSE), and mean unsigned error (MUE) for $ln{\gamma }_{i,j}^{\infty }$, and the average absolute percent deviation (AAPD) for ${\gamma }_{i,j}^{\infty }$ predicted using UNIFAC and MOSCED as compared to the reference data from “Yaws’ Handbook of Properties for Aqueous Systems” [56]. The results are separated for the case of water as the solute (water in organics) and water as the solvent (organics in water).

	Water As Solute		Water As Solvent
	UNIFAC	MOSCED	UNIFAC	MOSCED
$R^2$	0.713	0.923	0.816	0.962
RMSE	1.423	0.688	2.479	1.000
AAPD	51.395	55.749	77.925	138.245
MUE	0.824	0.532	1.647	0.761

.

Within the Yaws’ reference set where water is the solute, the SM12 universal continuum solvent model made predictions for 90 systems, UNIFAC modeled 87 systems, and MOSCED 53 systems. When water was the solvent, the SM12 universal continuum solvent model made predictions for 121 systems, UNIFAC 114 systems, and MOSCED 61 systems. MOSCED in both cases was the highest performer, yet modeled the lowest number of systems due to its limited compound inventory. We additionally note that the MOSCED parameters for water were optimized using limiting activity coefficients from “Yaws’ Handbook of Properties for Aqueous Systems” [75]. Overall, all models seem to have more difficulty predicting $ln{\gamma }_{i,j}^{\infty }$ when water is the solute. UNIFAC is the lowest performer in both cases. A plausible reason for the SM12 universal continuum solvent model’s improved accuracy over UNIFAC could be due to the cancellation of errors that occurs in computing $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ and $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$ consistently. Again, when water is the solvent with reference values of $ln{\gamma }_{i,j}^{\infty }$ > 9, we notice a tapering or an under prediction in $ln{\gamma }_{i,j}^{\infty }$ for both the SM12 and UNIFAC models. As stated previously, these are highly non-ideal systems.

The SMx universal solvent models, MOSCED, and UNIFAC are all limited in the system they can treat. MOSCED can only be used to model binary systems for which the molar volume and five molecular descriptors are available for each component. At present, MOSCED parameters are available for 130 organic solvents and water [20,75]. Efforts have been made to use electronic structure calculations to predict missing MOSCED parameters, including the use of solvation free energy calculations with the SMx universal solvent models [42,80,81,82,83,84]. While we would not expect the predicted MOSCED parameters to be as good as those regressed using experimental data, the results of the present study nonetheless would suggest they would yield acceptable results, especially for early stage process conceptualization and design applications. Interestingly, Phifer et al. [42] suggested that results obtained with MOSCED parameters predicted using solvation free energy calculations with the SMx universal solvent model were superior to predictions made with SMx alone. This may result from implicitly including reference data by using MOSCED parameters that were regressed using reference data.

UNIFAC is a group contribution-based method, and in general is able to cover a wider range of chemical compound space. Nonetheless, some parameters are still missing from the parameter matrix. Recently, it has likewise been demonstrated how electronic structure calculations, specifically COSMO-based (conductor-like screening model) calculations, can be used to predict missing UNIFAC parameters [85,86,87]. While the SMx universal solvent models, MOSCED, and UNIFAC are all limited in the system they can treat, there is opportunity for synergy among the methods.

5. Summary and Conclusions

The SMx (x = 12, 8, or D) universal solvent models are implicit solvent models that can be used to compute the solvation free energy at 298.15 K using electronic structure calculations. Fundamentally, the solvation free energy is an important thermophysical property, and corresponds to the transfer free energy of taking a single solute from a non-interacting ideal gas state to solution (at infinite dilution) at the same molar density. In going from a non-interacting state to a state with a single solute molecule in solution, the solvation free energy is representative of the solute-solvent interactions in solution. If the solvation free energy of a solute is computed in two different solvents, they can readily be used (via a thermodynamic cycle) to compute the transfer free energy between the two solvents. This transfer free energy is directly related to the partition coefficient, which can be used to relate the relative composition of the solute in the two (immiscible) solvents at equilibrium, or more generally can be used to compare the relative affinity of the solute for the two solvents.

In conventional chemical engineering thermodynamic modeling of phase equilibrium, we are interested in the activity coefficient of the components in the mixture, and the excess Gibbs free energy of the mixture. Coincidentally, values of the solvation free energy are not directly measured experimentally. Rather, it is common to measure limiting activity coefficients from which one may obtain the solvation free energy. The limiting activity coefficient is related to the transfer free energy of a single solute from a system of pure liquid solute to a solution phase at infinite dilution. Moreover, limiting activity coefficients may be related directly to the parameters in common binary interaction excess Gibbs free energy models used to model phase equilibrium.

In the present study, we therefore assessed the ability to use the SMx universal solvent models to predict limiting activity coefficients. We considered 148 molecules for which SMx (solvent) parameters were available in QChem 5.1.2. As seen in Equation (4), the calculation of limiting activity coefficients is dependent upon both the self-solvation free energy ($\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$) and solvation free energy terms ($\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$). Therefore we assessed the contribution of each individual term. The Minnesota Solvation Database was used as our reference set, since it was used to parameterize the SM12 model, which was consistently the top performer throughout the analysis.

In assessing the contribution of $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ and $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$ terms, water was found to be a large source of deviation. The $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ of water was predicted to be overly negative, suggesting we over-estimate the strength of the water-water interactions. We gained further insight into this contribution by breaking down the $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$ analysis into cases where water was the solute ($i=$ water), water was the solvent ($j=$ water), and where water was neither the solute or solvent. We found that when water was the solute, $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}/\left(RT\right)$ is predicted to be overly negative on average by 3.4 (dimensionless units) in the SM12 model as compared to the self-solvation free energy term $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}/\left(RT\right)$ which was overly negative by 3.8 (dimensionless units). This deviation ends up being advantageous in discerning how to calculate limiting activity coefficients, as there were two methods in which to do so: (1) prediction of $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ and $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$ and use of reference values for the molar volume, and (2) compute $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ using pure component vapor pressure data, which is readily available. We evaluated these two methods using three reference data sets: the Minnesota Solvation Database, Parts 1 to 6 of DECHEMA’s “Activity Coefficients at Infinite Dilution” volumes [50,51,52,53,54,55], and “Yaws’ Handbook of Properties for Aqueous Systems” [56]. All three methods told the same story, where at first glance there does not appear to be a significant difference in the predictions when $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ is predicted or computed from reference data. However, upon isolating systems where water is the solute, we see a that $ln{\gamma }_{i,j}^{\infty }$ is under-predicted when $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ is computed from reference vapor pressure data. Yet when both $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ and $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$ are predicted there is a convenient cancellation of errors in the $ln{\gamma }_{i,j}^{\infty }$. Therefore this suggests that the best method to calculate limiting activity coefficients moving forward is to compute $\mathsf{\Delta }{G}_{i,i}^{\mathrm{self}}$ and $\mathsf{\Delta }{G}_{i,j}^{\mathrm{solv}}$ consistently; that is, to predict both terms.

In conclusion, the findings were compared to two popular and successful predictive methods, UNIFAC and MOSCED. Comparison to these models allows us to assess the viability of taking the SMx models into further chemical engineering applications. A recent 2019 study conducted by Brouwer and Schuur [73], previously showed MOSCED was the top performing model in an comparative study to predict limiting activity coefficients at 298.15 K. Our comparison further confirms that MOSCED is the top performing model in comparison to UNIFAC and the SM12 universal continuum solvent model, yet it is limited by the number of compounds and systems it can model. UNIFAC and the SM12 universal continuum solvent model are comparable in performance among the three reference sets. This suggests that the SMx models have potential to be used and tested in chemical engineering applications such as the prediction of phase equilibrium. There is also opportunity for synergy between SM12, UNIFAC, and MOSCED.