Application of MOSCED To Predict Limiting Activity Coefficients

Article Cite This: J. Chem. Eng. Data 2018, 63, 352−364

pubs.acs.org/jced

Application of MOSCED To Predict Limiting Activity Coefficients, Hydration Free Energies, Henry’s Constants, Octanol/Water Partition Coefficients, and Isobaric Azeotropic Vapor−Liquid Equilibrium Pratik Dhakal, Sydnee N. Roese, Erin M. Stalcup, and Andrew S. Paluch* Department of Chemical, Paper and Biomedical Engineering, Miami University, Oxford, Ohio 45056, United States S Supporting Information *

ABSTRACT: Modified Separation of Cohesive Energy Density (MOSCED) is a solubility parameter-based method to predict limiting activity coefficients. In addition to making quantitative predictions, MOSCED may additionally be used to understand the underlying molecular-level driving forces for intuitive solvent selection and formulation. A major improvement of MOSCED over similar solubility parameter methods is that it splits the association term. We show by example how this change allows MOSCED to better model the molecular interactions of associating fluids. While parametrized to predict limiting activity coefficients, we demonstrate the ability to predict hydration-free energies, Henry’s constants in water, and octanol/water partition coefficients. We focus on water as MOSCED was previously found to perform substantially worse when water was the solvent. Comparison is made to molecular simulation and mod-UNIFAC, and we find for the studied reference set that predictions with MOSCED were in better agreement with experimental data. Comparison to mod-UNIFAC was additionally made to model binary isobaric azeotropic vapor−liquid equilibrium. Overall, for these nonideal systems, the quantitative agreement with experiment is better with mod-UNIFAC, although the difference with MOSCED is relatively small. An interactive software suite (MOSCED-CAD) has also been developed to calculate the properties described in the present study, in addition to binary interaction parameters for the Wilson and the nonrandom two-liquid (NRTL) equations for use in a process simulator.

■

Solubility parameter-based methods were first introduced in the Scatchard−Hildebrand regular solution theory. Regular solution theory assumes that the excess volume and entropy are zero, and therefore that the excess Gibbs free energy is equal to the excess internal energy.1,2 While the contribution to the excess Gibbs free energy due to the excess volume is typically minute, the contribution of the excess entropy may be significant, especially for cases in which the size and shape of the solvent and solute differ greatly. To remedy this shortcoming, one may add to the regular solution theory equation the athermal Flory−Huggins equation (or one of its various extensions), derived assuming the excess Gibbs free energy is equal to the excess entropy.1 For a binary mixture of species 1 and 2, this leads to the following expression for the limiting (or infinite dilution) activity coefficient of component 2, γ∞ 2 :

INTRODUCTION The ability to model and predict phase-equilibria, in addition to understanding the underlying molecular-level driving forces, is of central importance for the design of industrial separation processes, understanding the toxicity and the removal of environmental pollutants, the design of pharmaceuticals, and many other applications. Conventional methods commonly use group contribution methods to predict phase-equilibrium, which lack a molecular interpretation of the predicted phase behavior. State-of-the-art molecular simulations can provide molecular level insight, but are often computationally demanding. In the present study we investigate the use of solubility parameter-based methods. Solubility parameter-based methods have long been an important design tool for solvent selection and formulation. These methods allow one to (quantitatively) predict the activity coefficient of the components of a mixture using only pure component properties (i.e., solubility parameters). More importantly, the solubility parameter of component i is related to i−i intermolecular interactions. These methods therefore have the additional advantage of allowing one to understand (qualitatively) the behavior of a mixture by comparing solubility parameters (or intermolecular interactions). This allows one to intuitively design solvent mixtures based on intermolecular interactions. © 2018 American Chemical Society

ln γ2∞ =

⎛v ⎞ ⎛v ⎞ v2 (δ1 − δ2)2 + ln⎜ 2 ⎟ + 1 − ⎜ 2 ⎟ RT ⎝ v1 ⎠ ⎝ v1 ⎠

(1)

where R is the ideal gas constant, T is the absolute temperature, δi is the solubility parameter of component i, defined as the Received: August 22, 2017 Accepted: December 26, 2017 Published: January 26, 2018 352

DOI: 10.1021/acs.jced.7b00748 J. Chem. Eng. Data 2018, 63, 352−364

Journal of Chemical & Engineering Data

Article

⎡ q 2q 2(τ1(T ) − τ2(T ))2 v2 ⎢ (λ1 − λ 2)2 + 1 2 ψ1 RT ⎢⎣ (T ) (T ) (T ) (T ) ⎤ ⎛ v ⎞aa2 ⎛ v ⎞aa2 (α1 − α2 )(β1 − β2 ) ⎥ + + ln⎜ 2 ⎟ + 1 − ⎜ 2 ⎟ ⎥⎦ ξ1 ⎝ v1 ⎠ ⎝ v1 ⎠

square root of the component’s cohesive energy density, and vi is the liquid molar volume of component i, where i = {1, 2}. An equivalent expression may be written for component 1 by switching the subscript indices. Following the work of Abrams and Prausnitz,3 eq 1 may be interpreted as the sum of a residual (RES) and combinatorial (COMB) contribution: ln γ2∞ = ln γ2∞,RES + ln γ2∞,COMB

ln γ2∞ =

aa 2 = 0.953 − 0.002314[(τ2(T ))2 + α2(T )β2(T )]

(2)

⎛ 293K ⎞0.8 (T ) ⎛ 293K ⎞0.8 (T ) ⎛ 293K ⎞0.4 ⎟ ⎟ ⎟ αi(T ) = αi⎜ , βi = βi ⎜ , τi = τi⎜ ⎝ T ⎠ ⎝ T ⎠ ⎝ T ⎠ (5)

The first term of eq 1 (i.e., regular solution theory) is the residual contribution and results from intermolecular interactions; it corresponds to the partial molar excess internal energy. For all cases, ln γ∞,RES ≥ 0. The greater is the 2 dissimilarity in intermolecular interactions between two components, the greater is the value of ln γ∞,RES . The latter 2 terms (i.e., Flory−Huggins theory) are the combinatorial contribution, and results from the size dissimilarity of the components, which corresponds to the partial molar excess entropy. For all cases, ln γ∞,COMB ≤ 0. The greater the size 2 dissimilarity between the two components, the more negative ∞ the value of ln γ∞,COMB . An ideal solution (ln γ∞ 2 2 = ln γ1 = 0) results when the components have similar intermolecular interactions (δ1 = δ2) and are comparable in size (v1 = v2). While eq 1 has been successfully used for a wide range of applications, its major limitation is that the residual term only accounts for molecules interacting via dispersion interactions. Higher order interactions may be accounted for by rewriting the cohesive energy density as a sum of contributions (or partial solubility parameters).2,4−6 Perhaps the most wellknown expansion is that of the Hansen Solubility Parameters (HSP):5,7 δi2 = δi2,D + δi2,P + δi2,H = λi2 + τi2 + (αiβi )2

where i = {1 or 2} ψ1 = POL + 0.002629α1(T )β1(T ) ξ1 = 0.68(POL − 1) 2 + [3.4 − 2.4 exp(− 0.002687(α1β1)1.5 )](293K/ T ) POL = q14[1.15 − 1.15 exp(− 0.002337(τ1(T ))3 )] + 1

where an equivalent expression can be written for component 1 by switching the subscript indices. The induction parameter, qi, reflects the ability of the nonpolar part of a molecule to interact with a polar part. The terms ψ1 and ξ1 are empirical (solvent dependent) asymmetry terms to modify the residual contribution for polar and hydrogen-bonding interactions, and aa2 is an empirical (solute dependent) term to modify the size disimilarity in the combinatorial contribution for polar and hydrogen bonding interactions. These additional empirical terms are not adjustable but are functions of the other parameters (τi, αi, βi, and qi). For all cases aa2 ≤ 0.953, effectively reducing the size dissimilarity and magnitude of ln γ∞,COMB , with the value smaller for polar and associating 2 compounds. The superscript (T) is used to indicate temperature dependent parameters, where the temperature dependence is computed using the empirical correlations provided in eq 5. As suggested by the equations, MOSCED adopts a reference temperature of 293 K (20 °C). An important difference between MOSCED and HSP is that by splitting the association (or hydrogen-bonding) term, the contribution to ln γ∞,RES due to association (or hydrogen2 bonding interactions) may be either positive or negative. The term resulting from intermolecular interactions, ln γ∞,RES , may 2 now contribute toward both positive and negative deviations from Raoult’s law. This subtle difference allows MOSCED to capture the underlying physics of solvation of associating compounds. This will be explored in greater detail in the subsection “Association/Hydrogen Bonding Term,” including a discussion of the system acetone/chloroform. Moreover, in the present study we relate γ∞ 2 to the solute solvation free energy. This relation between solvation free energy and γ∞ 2 is valuable for two reasons. First, the solute solvation free energy is a direct measure of solute−solvent intermolecular interactions, which will allow us to focus solely on MOSCED’s ability to capture the challenging solute− solvent intermolecular interactions. The (Lewis/Randall normalized) limiting activity coefficient is a measure of solute−solvent intermolecular interactions, however, these interactions are normalized with respect to the interaction of the solute with itself (i.e., solute−solute interaction). Second, the solvation free energy is a quantity from statistical mechanics that may readily be computed using molecular simulation free energy calculations or electronic structure calculations in a

(3)

where δi,D, δi,P, and δi,H correspond to the solubility parameter resulting from dispersion, polar, and association (or hydrogen bonding) interactions, respectively. We equivalently write our dispersion and polar solubility parameters as λi and τi, respectively, and write the association solubility parameter as the product of the contribution due to hydrogen-bond acidity (αi) and basicity (βi). This results in ln γ2∞,RES =

v2 [(λ1 − λ 2)2 + (τ1 − τ2)2 RT + (α1β1 − α2β2)2 ]

(4)

where we still have for all cases ln γ∞,RES ≥ 0. It follows that ln 2 γ∞,RES , the term resulting from intermolecular interactions, can 2 only contribute toward positive deviations from Raoult’s law, while ln γ2∞,COMB can only contribute toward negative deviations from Raoult’s law. As an example of where this is not correct, consider the binary system acetone/chloroform which is well-known to exhibit negative deviations from Raoult’s law and a maximum boiling azeotrope.8 Recently, using neutron diffraction, this was confirmed to be caused by the association of acetone and chloroform molecules.9 When HSP is used, ln γ∞,RES and the association term may only 2 contribute toward positive deviations from Raoult’s law and is therefore unable to capture the underlying physics of this system. A major difference between MOSCED (Modified Separation of Cohesive Energy Density) and HSP is that MOSCED separates the association term, resulting in the following system 10,11 of equations for γ∞ 2 in a binary mixture: 353



Article

v1(T , p)p2sat (T ) 1 ΔG2solv (T , p) = ln γ2∞(T , p) + ln RT RT

continuum solvent. This therefore provides a route by which molecular simulation and electronic structure calculations may be used to parametrize MOSCED, creating a truly predictive method,12−16 or to use MOSCED as a tool to guide molecular simulation and to help interpret results. Moreover, reference solvation free energy values are valuable in the optimization of molecular simulation force fields; the importance of such has resulted in the recent development of the FreeSolv and CompSol Databanks for this purpose.17−19 In the present study, we use hydration free energies (i.e., solvation free energies in water) and Henry’s constants for solutes in water to demonstrate the ability of MOSCED to accurately model solute−water interactions. Additionally, using octanol/water partition coefficients, we demonstrate the ability of MOSCED to properly balance solute hydrophobic and hydrophilic interactions. The hydration free energies are compared to state-of-the-art molecular simulations,17,18,20 and all calculations are compared to mod-UNIFAC.21−26 In all cases, MOSCED is found to agree best with available experimental reference data. We then further use MOSCED to predict binary isobaric (Txy) azeotropic vapor−liquid equilibrium (VLE). For these nonideal systems, predictions made using MOSCED are in good quantitative agreement with experiment. Although modUNIFAC performs slightly better, use of MOSCED is advantageous as it directly compares the intermolecular interactions of the components in the mixture, shedding insight into the underlying molecular driving forces responsible for the observed phase behavior. To facilitate the adoption and use of MOSCED, an interactive software suite (MOSCED-CAD) has been developed. MOSCED-CAD can be used to calculate all of the properties described in the present study, in addition to binary interaction parameters for the Wilson and nonrandom twoliquid (NRTL) equations for use in a process simulator.1 A copy of the current user manual is provided in the Supporting Information accompanying the electronic version of this manuscript, which includes instructions for downloading from our personal Web site.27

(8)

Note that we define the solvation free energy as the change in free energy of taking a solute from an ideal gas to solution at the same molecular density (or concentration). Henry’s Constant. The activity coefficient of the solute can be related to the Henry’s constant of the solute (/2,1) by equating the solution phase fugacity computed using a Lewis/ Randall (Raoult’s law) and Henry’s law standard state:1 f2L (T , p , x 2) = x 2γ2(T , p , x 2)f 20 (T , p) = x 2γ2*(T , p , x 2)/2,1(T , p)

where γ*2 is the (asymmetical) activity coefficient of the solute normalized with respect to a Henry’s law standard state; in the infinite dilution limit, γ2* = 1. At low pressures we may again assume that the vapor phase in equilibrium with the liquid phase at T is an ideal gas (such that the fugacity coefficient is unity) and the Poynting correction is negligible such that f L2 ≈ psat 2 . Solving for / 2,1 in the infinite dilution limit we have /2,1(T , p) = γ2∞(T , p)p2sat (T )

ln /2,1(T , p) =

x 2Iγ2I(x 2I , T , p)f 20 (T , p) = x 2IIγ2II(x 2II , T , p)f 20 (T , p)

K 2I/II(T ,

(12)

p) =

x 2I x 2II

=

γ2II, ∞(T , p) γ2I, ∞(T , p)

(13)

The partition coefficient may equivalently be seen as the relative solubility in solvent I relative to solvent II. While KI/II 2 is a dimensionless quantity, it is dependent on the units of concentration used. It is common to use molar (mol/volume) or mass concentrations (mass/volume) instead of mole fractions. Having assumed the solute is infinitely dilute, we can readily convert from mole fractions to molar or mass concentrations as cI2 = xI2/vI and cII2 = xII2 /vII. This results in the final expression for the partition coefficient using molar or mass concentrations:

v (T , p ) 1 ΔG2solv (T , p) − ln 1 RT RT

v (T , p) 1 ΔG2self (T , p) − ln 2 RT RT

(11)

Solving for the relative concentration in phase I relative to phase II, and assuming the solute concentration in each solvent is sufficiently small so as to be considered infinitely dilute, we obtain the solute partition coefficient (KI/II 2 ):

(6)

where f 02 is the pure liquid fugacity of the solute. The pure liquid fugacity may be related to the solute “self”-solvation free 12 energy (ΔGself 2 ) as ln f 20 (T , p) =

1 RT ΔG2solv (T , p) + ln v1(T , p) RT

Partition Coefficient. The distribution of a solute between two immiscible solvents may be described by the solute partition coefficient.1 Equating the solute fugacity in solvent I and II at equilibrium, we have

COMPUTATIONAL METHODS Solvation Free Energy. In the context of classical molecular simulation free energy calculations, we have shown previously that the solvation free energy of a solute (component 2, ΔGsolv 2 ), or equivalently the residual chemical potential at infinite dilution, may be related to the limiting (or infinite dilution) activity coefficient of the solute as12,28−30

− ln f 20 (T , p)

(10)

Following the previous section, the Henry’s constant can be related to the solute solvation free energy as31

■

ln γ2∞(T , p) =

(9)

P2I/II(T , p) = (7)

c 2I c 2II

=

γ2II, ∞(T , p) vII(T , p) γ2I∞(T , p) vI(T , p)

(14)

The partition coefficient can equivalently be computed using solvation free energies as32,33

At low pressures we may assume the vapor phase in equilibrium with the liquid phase at T is an ideal gas (such that the fugacity coefficient is unity) and the Poynting correction is negligible sat such that f L2 ≈ psat 2 , where p2 is the (pure component) liquid 1 saturation pressure. Solving eq 6 for ΔGsolv 2 , we have

ln P2I/II(T , p) =

1 1 ΔG2II,solv (T , p) − ΔG2I,solv (T , p) RT RT (15)

354



Article

Written in this form, we see that ln PI/II 2 is a direct measure of solute−solvent II intermolecular interactions relative to solute− solvent I intermolecular interactions. Binary Vapor−Liquid Equilibrium. At low pressures we may again assume the vapor phase in equilibrium with the liquid phase at T is an ideal gas (such that the fugacity coefficient is unity) and the Poynting correction is negligible, resulting in the “modified”-Raoult’s law (γ − ϕ approach) to model VLE. Equating the fugacity of each species i in the liquid and vapor phase, we have1 xiγi(xi , T , p)pisat = yp i

For all VLE calculations, the pure component vapor pressure of each component was computed using the Antoine equation parameters from ref 37. When predicting isothermal VLE (pxy), MOSCED calculations were performed at the temper∞ ature of the reference data (ln γ∞ 1 and ln γ2 ), which were then used to obtain Λ12 and Λ21 for Wilson’s equation. When predicting isobaric VLE (Txy), MOSCED calculations were performed at 10 equally spaced temperatures between the saturation temperature of component 1 and 2. At each ∞ temperature, ln γ∞ 1 and ln γ2 were used to solve for the BIPs for Wilson’s, TK−Wilson, and NRTL (α12 = 0.3) equations. While all three equations contain a temperature dependence, we performed a linear correlation of the temperature dependent BIPs (i.e., a12 = c0+c1T). With the Wilson’s and TK−Wilson equations, we assume the ratio of molar volumes is insensitive to temperature, and use the molar volumes used by MOSCED at 20 °C.11,38 VLE was then computed by specifying x1 over the range 0 to 1 in increments of 0.05, and then solving for y1 and p or T, for the isothermal and isobaric case, respectively. While we used an equal spacing of x1 for our calculations, this was not necessarily the case in the reference data. To quantitatively compare the predicted values of y1 and p or T to the reference data, we used (1-D) spline interpolation to interpolate the predicted values of y1 and p or T with respect to x1 using GNU Octave.39,40 Spline interpolation was used to make predictions at the values of x1 where reference data was available, allowing a direct comparison of the predicted values of y1 and p or T to the reference data. For all of the calculations (of all properties), reference calculations were additionally performed with mod-UNIFAC.21−26 All of these calculations were performed using CHEMCAD 7.1.0.9402.36 For VLE calculations, calculations were performed in increments of x1 of 0.05 mole fractions. To calculate solvation free energies, Henry’s constants, and partition coefficients, mod-UNIFAC was used to predict limiting activity coefficients. For all solvation free energy and Henry’s constant calculations (using MOSCED and modUNIFAC), vapor pressure data were taken from ref 37. The necessary molar volumes were assumed constant and taken directly from MOSCED.11,38 All of the calculations (of all properties) using MOSCED and mod-UNIFAC are tabulated in the Supporting Information of this manuscript along with the reference data to which comparison is made. We additionally tabulate the MOSCED11 and Antoine parameters37 for all 130 organic solvents and water for which MOSCED is parametrized.

(16)

where yi is the vapor-phase mole fraction of component i and i = {1, 2}. While MOSCED predicts limiting activity coefficients, composition-dependent activity coefficients are required for eq 16. As suggested by Schreiber and Eckert,34 the infinite dilution activity coefficients of a binary pair can be used to parametrize a binary interaction excess Gibbs free energy model. In the present study we will investigate the use of Wilson’s, TK− Wilson, and NRTL (α12 = 0.3) equations.1,35 For Wilson’s equation, we have1 ⎛ ⎞ Λ 21 Λ12 − ln γ1 = −ln(x1 + Λ12x 2) + x 2⎜ ⎟ Λ 21x1 + x 2 ⎠ ⎝ x1 + Λ12x 2 ⎛ ⎞ Λ12 Λ 21 − ln γ2 = −ln(x 2 + Λ 21x1) − x1⎜ ⎟ Λ 21x1 + x 2 ⎠ ⎝ x1 + Λ12x 2 (17)

where Λ12 and Λ21 are adjustable parameters which may be related to the binary (intermolecular) interaction parameters (BIPs) of the system (a12 and a21).

Λ12 =

v2 ⎡ a12 ⎤ exp⎢ − ⎥ v1 ⎣ RT ⎦

Λ 21 =

v1 ⎡ a 21 ⎤ exp⎢ − ⎥ v2 ⎣ RT ⎦

(18)

At infinite dilution, eq 17 reduces to ln γ1∞ = −ln(Λ12) + 1 − Λ 21 ln γ2∞ = −ln(Λ 21) + 1 − Λ12

γ∞ 1

(19)

γ∞ 2

When ln and ln are calculated using MOSCED, eq 19 provides a system of two equations with two unknowns (Λ12 and Λ21 or a12 and a21). Equation 17 may be used to calculate composition-dependent activity coefficients when solving for Λ12 and Λ21. Similar expressions for the TK−Wilson and NRTL equations are provided in the Supporting Information accompanying this manuscript. Note that we have not evaluated the universal quasi-chemical model (UNIQUAC) as we seek to use methods that may be parametrized solely using information from MOSCED.3 While the present study is limited to binary VLE, this approach may readily be generalized to model multicomponent VLE. The Wilson’s, TK−Wilson, and NRTL equations all have two binary parameters for each pair in the multicomponent system, in which MOSCED may be used to parametrize each binary pair.1 This is the approach used by the process simulator CHEMCAD,36 for which binary interaction parameters optimized for a binary system are combined to model multicomponent systems.

■

RESULTS AND DISCUSSION Association/Hydrogen Bonding Term. As already mentioned, solubility parameter-based methods are powerful in that they can both make quantitative predictions and give qualitative, molecular-level insight which may be used for intuitive solvent selection and formulation. Comparing MOSCED and HSP, a major difference is MOSCED separates the association (or hydrogen bonding) term in ln γ∞,RES . We 2 will illustrate how this is able to better describe the underlying physics of associating fluids using three examples: acetone/ chloroform, methanol/acetone, and methanol/N-methylpyrrolidone. In this discussion remember that ln γ∞,RES results from 2 intermolecular interactions, and compares 2−1 (solute− solvent) interactions relative to solute−solute (2−2) interactions. MOSCED αi and βi values may be found in Table 1. 355



Article

Table 1. MOSCED Association Parameters11 for Chloroform, Acetone, Methanol, and N-Methylpyrrolidone, Where αi and βi Have Units of MPa1/2 or (J/cm3)1/2 species (i)

αi

βi

chloroform acetone methanol N-methylpyrrolidone

5.80 0.00 17.43 0.00

0.12 11.14 14.49 24.22

αβself

(21)

In this light, we observe positive deviation from Raoult’s law for binary mixtures in which the components prefer to associate with themselves, where the mean self-association term is greater than the mean cross-association term (αβself > αβcross). On the other hand, we observe negative deviation from Raoult’s law for binary mixtures in which the components prefer to associate with each other rather than themselves, where the mean crossassociation term is greater than the mean self-association term (αβcross > αβself ). This can be compared directly to HSP where the squared difference of the self-association term of each species is computed. In the top pane of Figure 1 we plot the pxy phase-diagram for acetone(1)/chloroform(2). The system exhibits negative

Here our focus is on the ability of MOSCED to capture positive and negative deviations from Raoult’s law, and the interpretation of the results. While it is possible to observe deviations from Raoult’s law using isobaric VLE, we will consider isothermal VLE for which it is easier to visually observe the deviations. First, consider the binary system of acetone(1)/chloroform(2). Chloroform has three weak C−H hydrogen-bond donating sites (proton donor, α2), and a weaker C−Cl halogen bonding site (proton donor and acceptor, α2 and β2). On the other hand, acetone has a moderately strong CO hydrogen bond accepting site (β1).41 We find that α2 > α1 while β2 < β1, indicating that chloroform would prefer to associate with acetone and acetone would prefer to associate with chloroform, than either component with itself. Favorable cross-interactions are expected to lead to negative deviations from Raoult’s law (ln γ1 < 0 and ln γ2 < 0). Looking at the association term in ln , we find that (α1 − α2)(β1 − β2) < 0. It follows that the γ∞,RES 2 association term is negative; the favorable cross-interactions result in a negative contribution to ln γ1∞ and ln γ2∞, contributing toward negative deviations from Raoult’s law. The self-association of each species is given by the product α1β1 and α2β2. For acetone, this product is 0, indicating that it does not associate with itself. For chloroform, this product is 0.696 J/cm3, suggesting weak association (C−H···Cl−C and C−Cl···Cl−C) interactions. It is exactly this difference in selfassociation terms that is compared in HSP: (α1β1 − α2β2)2. ∞ This results only in contributions to ln γ∞ 1 and ln γ2 that are ∞,RES and ln greater than 0. HSP can only predict values of lnγ1 greater than or equal to zero. Therefore, the only term γ∞,RES 2 that can contribute to negative deviations from Raoult’s law is the combinatorial term (COMB) which results from size and shape dissimilarities, and is not due to intermolecular interactions. We see that as compared to HSP, the splitting of the association term in MOSCED allows us to better capture the underlying physics of the process, which we ultimately expect to lead to better quantitative predictions for associating fluids. We may additionally obtain an alternative (but equivalent) interpretation of the MOSCED association term by expanding and rearranging it:

Figure 1. Plots of pxy diagrams for acetone(1)/chloroform(2) at 298.15 K, acetone(1)/methanol(2) at 323.15 K, and methanol(1)/Nmethylpyrrolidone(2) at 393.15 K. Circles are experimental reference data,8,42,43 and solid lines correspond to predictions made using MOSCED11 to parametrize Wilson’s equation. The dotted red line corresponds to the bubble line predicted using Raoult’s law and is drawn for convenience.

deviations from Raoult’s law and a minimum boiling azeotrope, with the predictions made using MOSCED in close agreement with experiment.8 Observed recently using neutron scattering, the dominant contribution to the minimum boiling azeotrope is the strong association (hydrogen bonding) between chloroform and acetone.9 This is in exact agreement with MOSCED; α2β1 = 64.612 J/cm3, which is 2 orders of magnitude greater than either self-association term. Consider next the binary system of methanol(1)/acetone(2). The O−H group of methanol is a moderately strong hydrogenbond donor and acceptor, and we find that α1 > α2 and β1 > β2. We therefore expect methanol to prefer to associate with itself rather than with acetone, which would result in positive deviations from Raoult’s law (ln γ1 > 0 and ln γ2 > 0). This is

(α1 − α2)(β1 − β2) = −(α1β2 + α2β1) + (α1β1 + α2β2) = −2(αβcross − αβself )

1 (α1β2 + α2β1) 2 1 = (α1β1 + α2β2) 2

αβcross =

(20)

where we have defined αβcross and αβself to be the mean “cross” and “self” association term (or energy), respectively: 356



Article

exactly captured by MOSCED where (α1 − α2) (β1 − β2) > 0. In the middle pane of Figure 1 we plot the pxy phase-diagram for methanol(1)/acetone(2), and find the predictions made using MOSCED result in positive deviations from Raoult’s law in agreement with experiment.42 It follows that if we seek a species (component 2) that when mixed with methanol (component 1) results in negative deviations from Raoult’s law, we need to find a species such that either α2 > α1 or β2 > β1, such that (α1 − α2) (β1 − β2) < 0. One candidate that fits this requirement is N-methylpyrrolidone, for which β2 > β1. In the bottom pane of Figure 1 we plot the bubble line (px) for methanol(1)/N-methylpyrrolidone(2), and find the predictions made using MOSCED result in negative deviations from Raoult’s law in agreement with experiment.43 Limiting Activity Coefficient, Solvation Free Energy, and Henry’s Constant. In the 2005 MOSCED parametrization, parameters were regressed for 130 organic solvents using 6441 reference limiting activity coefficients. The root mean squared error (RMSE) for ln γ∞ 2 was 0.148, and the 11,44 average absolute error (AAE) for γ∞ With the 2 was 10.6%. MOSCED parameters for the 130 organic solvents set, the authors then regressed MOSCED parameters for water using reference limiting activity coefficients for organic solutes in 11,44 water, and obtained a much larger AAE for γ∞ 2 of 41.1%. Given the inferior performance for water, we will consider the specific case of solvation in water. We first compare the ability of MOSCED to predict ln γ∞ 2 for organics in water to reference data from ref 45, mostly at 25 ◦C. Of the systems in the MOSCED database, predictions could be made for 88 solutes. The results of the predictions are summarized in Table 2 and

Figure 2. Parity Plot of ln γ∞ 2 for organics in water predicted using MOSCED11 and mod-UNIFAC21−26,36 versus reference values from ref 45. The dashed lines indicate an error of ±2 log units.

Table 3. Summary of the Number of Solutes (N) for which Hydration Free Energies Were Calculated, Average Absolute Error (AAE) and the Corresponding Standard Deviation (SD), and Pearson’s Correlation Coefficient (R2) of the Calculations versus the Reference Data,17,18,20 Where All Hydration Free Energy Calculations Were in Units of kJ/ mol

Table 2. Summary of the Number of Solutes (N) for Which Log Limiting Activity Coefficients (ln γ∞ 2 ) Were Calculated in Water, the Average Absolute Error (AAE) and the Corresponding Standard Deviation (SD), and Pearson’s Correlation Coefficient (R2) of the Calculations versus the Reference Data45 method

N

AAE

SD

R2

MOSCED mod-UNIFAC

88 58

0.83 2.31

0.75 2.15

0.96 0.78

method

N

AAE

SD

R2

MOSCED molecular simulation mod-UNIFAC

73 73 60

2.04 2.85 7.95

2.73 3.69 9.31

0.94 0.94 0.70

much more computationally expensive and more detailed molecular simulation results, with MOSCED resulting in a smaller overall AAE. The FreeSolv database additionally includes the uncertainty of the experimental measurements.20 For the set of 73 solutes studied here, the mean uncertainty is 2.33 kJ/mol with an SD of 0.51 kJ/mol. We find that the mean uncertainty in the experimental data is greater than the AAE for the MOSCED predictions. With mod-UNIFAC, of the 73 solutes examined using MOSCED, groups were available to model 60 of the solutes. For this reduced set we obtain an AAE of 7.95 kJ/mol with an SD of 9.31 kJ/mol, and R2 = 0.70. While in the 2005 MOSCED parametrization MOSCED was found to perform the worst when water was the solvent, we find here that the predicted hydration free energies are superior to mod-UNIFAC and comparable to results obtained using molecular simulation. This is additionally illustrated in the parity plot of the predictions versus the reference data provided in Figure 3. Likewise, we next compute log Henry’s constants (ln /2,1) for organic solutes (2) in water (1). From eq 11, we see that ln /2,1 is equal to the dimensionless solvation free energy (here hydration free energy) plus the log of the solvent molar volume. ln /2,1 may therefore additionally be used as a direct

Figure 2. We obtain an AAE for ln γ∞ 2 of 0.83 with a standard deviation (SD) of 0.75, and a Pearson’s correlation coefficient of R2 = 0.96. Predictions could be made for 58 of the systems using mod-UNIFAC. With this reduced set we obtain an AAE of 2.31 with an of 2.15, and R2 = 0.78. Next we focus the calculation of solvation free energies, which are a direct measure of solute−solvent intermolecular interactions. For the specific case of solvation in water, this corresponds to hydration free energies. We compare predictions using MOSCED to reference values at 25 °C from the FreeSolv database,17,18,20 with the results summarized in Table 3. Of the systems in the database, MOSCED predictions could be made for 73 solutes. We obtain an AAE of 2.04 kJ/mol with an SD of 2.73 kJ/mol, and R2 = 0.94. The FreeSolv database also contains predictions using state-of-theart molecular simulation methods and has been used to assess the accuracy of conventional molecular models. For the same set of 73 solutes, the simulation results obtain an AAE of 2.85 kJ/mol with an SD of 3.69 kJ/mol, and R2 = 0.94. We find that the predictions made using MOSCED are comparable to the 357



Article

Figure 3. Parity plot of the solute solvation free energy in water 11 mod(hydration free energy, ΔGsolv 2 ), predicted using MOSCED, 21−26,36 17,18,20 and molecular simulation versus reference UNIFAC, data.17,18,20 The dashed lines indicate an error of ±5 kJ/mol or approximately 2RT.

Figure 4. Parity plot of log Henry’s constants (ln /2,1) for organics in water, predicted using MOSCED11 and mod-UNIFAC21−26,36 versus reference data from ref 45. The dashed lines indicate an error of ±2 log units.

measure of solute−solvent intermolecular interactions. This affords us the opportunity to both demonstrate the calculation of an additional property and to compare to an independent database from ref 45. As compared to our hydration free energy calculations, with both MOSCED and mod-UNIFAC, we are able to make predictions for an additional eight solutes. When ln /2,1, was calculated, /2,1 had units of kPa. The results of predicting ln /2,1 using both MOSCED and mod-UNIFAC are summarized in Table 4 and are shown

volumes of reference data.46 More importantly, in the present study, octanol/water parition coefficients allow us to probe the ability of MOSCED to model the intermolecular interactions between the solute and water relative to the solute and 1octanol (a lipophilic reference solvent). We will use eq 14 and ignore the mutual solubility of water and 1-octanol, consistent with refs 46 and 47. We compare predictions of lnPI/II using MOSCED to 2 reference values at 25 °C from ref 45, with the results summarized in Table 5. Of the systems in the database,

Table 4. Summary of the Number of Solutes (N) for which log Henry’s Constants (ln /2,1) Were Calculated in Water, Average Absolute Error (AAE), and the Corresponding Standard Deviation (SD), and Pearson’s Correlation Coefficient (R2) of the Calculations versus the Reference Data,45 Where /21 Was Calculated in Units of kPa

Table 5. A Summary of the Number of Solutes (N) for which log Octanol(I)/Water(II) Partition Coefficients (ln PI/II 2 ) Were Calculated, The Average Absolute Error (AAE) and the Corresponding Standard Deviation (SD), and Pearson’s Correlation Coefficient (R2) of the Calculations versus the Reference Data,45 where Molar Concentrations (mol/ Volume) Were Used in the Calculation of ln PI/II

method

N

AAE

SD

R2

MOSCED mod-UNIFAC

81 68

0.74 1.44

1.07 1.93

0.97 0.71

graphically in the parity plot of the predictions versus the reference data in Figure 4. With MOSCED, we obtain an AAE of 0.74 with an SD of 1.07, and R2 = 0.97. With mod-UNIFAC, we obtain an AAE of 1.44 with an SD of 1.93, and R2 = 0.71. The results for ln /2,1 are consistent with our hydration free energy and ln γ2∞ calculations; predictions made using MOSCED have a smaller AAE and SD as compared to modUNIFAC, and have an R2 value closest to 1. Partition Coefficients. Next, we consider the ability of MOSCED to predict octanol/water (I/II) partition coefficients (lnPI/II 2 ). The study of octanol/water partition coefficients are of particular interest because of their utility to correlate biological activity, which has resulted in compilations of

method

N

AAE

SD

R2

MOSCED mod-UNIFAC

117 96

0.92 1.95

1.43 3.04

0.93 0.60

MOSCED predictions could be made for 117 solutes. We obtain an AAE of 0.92 with an SD of 1.43, and R2 = 0.93. Of the 117 solutes examined using MOSCED, mod-UNIFAC groups were available to model 96 of the solutes. For this reduced set we obtain an AAE of 1.95 with a SD of 3.04, and R2 = 0.60. The results are also summarized graphically in Figure 5. Notice in Figure 5 that the predictions break from the trend I/II for reference values of lnPI/II 2 < 0. A value of lnP2 < 0 would indicate a hydrophilic solute; a solute that prefers water. For these cases the MOSCED predictions appear to overpredict the solute affinity for water (hydrophilicity). However, we can understand this deviation if we consider as an example ternary liquid−liquid equilibrium with water/methanol/1-octanol at 25 358



Article

searching through parts 1, 2a, 2b, 2h, and 7 of the DECHEMA “Vapor-Liquid Equilibrium Data Collection”51−55 and the DECHEMA “Recommended Data of Selected Compounds and Binary Mixtures”56 for all possible binary pairs that could be modeled using the existing MOSCED parameters for 130 organic solvents plus water.11,44,57 The agreement between the predictions and reference data was quantified by comparing the computed volatility of component 1, K1, when component 1 is dilute, where here we take x1 < 0.3 to mean dilute.58,59 We also compare the computed volatility of component 2, K2, when component 2 is dilute (x2 < 0.3), where Ki =

yi xi

= γi

pisat p

(22)

We see that the volatility is directly proportional to the activity coefficient, where p is constant and psat i is computed using reference Antoine equation parameters from ref 37. A plot of the relative volatility (α1,2 = K1/K2), K1, and K2 versus x1 is provided in Figure 6 for the representative system cyclohexane(1)/1,2-dichloroethane(2). This system was chosen because it had the largest number of reference data. We find in the plots for K1 and K2 that the predictive methods exhibit the largest deviation from the reference data in the dilute region, with all of the predictions in close agreement with the reference data outside this region. We additionally find in the plot for α1,2 that all of the predictive methods are able to predict the composition of the azeotrope in good agreement with the reference data. Txy and McCabe−Thiele (yx) diagrams are also provided in the Supporting Information of this manuscript for four representative systems. A summary of the quantitative agreement between the predictions and reference data is provided in Table 6. The percent average absolute relative error (AARE) is computed as

Figure 5. Parity plot of log octanol(I)/water(II) parition coefficients 11 and mod-UNIFAC21−26,36 (ln PI/II 2 ), predicted using MOSCED versus reference data from ref 45. The dashed lines indicate an error of ±2.8 log units. This is indicative of an error of approximately 2RT in and ΔGII,solv . ΔGI,solv 2 2

°C. Starting with the binary water/1-octanol, the aqueous-phase is pure water while the organic phase is 0.732 mole fractions of 1-octanol. As the amount of methanol added to this system increases, the amount of water in the organic phase also increases. When we have a mole fraction of methanol of 0.3995 in the organic phase, the solute-free mole fraction of octanol has decreased to 0.2390 mole fractions. At the same time, the solute-free mole fraction of water in the aqueous phase is 0.9815, very close to 1.48 We expect similar behavior for the other hydrophilic solutes. For these cases, the assumptions of immiscible solvents and an infinitely dilute solute used to derive eq 14 are not appropriate. The actual value of the reported partition coefficient will also be sensitive to the conditions used to make the measurement.48 Likewise, since the assumption of a pure 1-octanol phase is inadequate, ref 49 has cautioned against using octanol/water partition coefficients as a means to compute limiting activity coefficients. Nonetheless, for these systems MOSCED still predicts the correct qualitative behavior that the solutes prefer water over 1-octanol. Overall, we find that MOSCED is able to model a wider range of solutes and offers superior predictions. While the MOSCED parameters for the 117 solutes for which predictions were made were regressed using experimental data, this result is encouraging as it suggests that MOSCED is able to accurately balance the hydrophilic and hydrophobic interactions of a molecule. Azeotropic Vapor/Liquid Equilibrium. Last, we consider the use of MOSCED to predict binary isobaric (Txy) azeotropic VLE. Azeotropic VLE results from deviations from ideality, with either favorable of disfavorable cross interactions. Additionally, the prediction of isobaric VLE offers the additional challenge of accurately calculating temperaturedependent activity coefficients. We looked at a total of 17 binary systems,50 which are summarized in Table 6. These systems were selected by

N

AARE =

∑ i=1

|K iref − K ipred| K iref

(23)

where the superscript “ref” and “pred” indicated reference and predicted data, respectively, where the summation is over all N compositions for which reference data is available for comparison. A summary of agreement between the predictions and reference data in terms of error in y1 and T is additionally provided in the Supporting Information. For the MOSCED predictions, we find that the best overall predictions are made using Wilson’s equation. This is consistent with ref 60 which found for the entire DECHEMA “Vapor-Liquid Equilibrium Data Collection”, that Wilson’s equation was on average able to best correlate VLE data. Comparing our MOSCED predictions using Wilson’s equation to mod-UNIFAC, we find that modUNIFAC on average performs better. However, the difference is small. The average differance in K1 and K2 in the dilute region is 5.5 and 2.7%, respectively. The average difference in AAE for y1 and T is 0.009 mol fractions and 0.568 K, respectively. Comparing to mod-UNIFAC, we therefore conclude that MOSCED with Wilson’s equation performs well, especially given the nonideality of the systems studied. Earlier, we found that MOSCED predictions of limiting activity coefficients in water, hydration free energies, Henry’s constants in water, and octanol/water partition coefficients were in better agreement with reference data than modUNIFAC. All of these calculated properties were directly 359



Article

Table 6. A Summary of the Percent Average Absolute Relative Error (AARE) in the Volatility of Components 1 (K1) and 2 (K2) in the Dilute Region, x1 < 0.3 and x2 < 0.3, Respectively, of the Predictions versus Reference Data51−56a Ki

system cyclohexane/1,2-dichloroethane

ethanol/1,2-dichloroethane

ethanol/water

1-propanol/water

2-propanol/water

2-butanone/water

acetonitrile/water

benzene/2-propanol

chloroform/methanol

cylohexanone/phenol

ethyl-acetate/water

methanol/benzene

methanol/cyclohexane

water/pyridine

ethanol/acetonitrile

benzene/acetonitrile

methanol/tetrahydrofuran

average

1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1

and 2

and 2

and 2

and 2

and 2

and 2

and 2

and 2

and 2

and 2

and 2

and 2

and 2

and 2

and 2

and 2

and 2

and 2

mod-UNIFAC

NRTL

Wilson

TK−wilson

N

8.644 3.943 6.294 5.307 4.480 4.831 2.666 1.033 2.258 2.341 2.293 2.323 25.266 3.590 15.632 30.151 6.870 22.391 30.924 15.612 23.620 9.219 26.027 17.623 3.103 2.788 2.958 52.045 46.290 50.123 45.123 64.983 53.635 5.944 11.538 8.741 7.025 29.209 19.125 18.761 45.916 43.201 1.089 2.385 1.459 2.036 10.168 4.538 −50 −50 −50 15.603 17.327 17.422

3.446 6.643 5.044 8.055 13.025 10.541 15.578 7.407 13.670 18.509 18.414 18.473 55.217 11.224 33.321 37.359 4.586 25.442 42.422 25.29 34.637 6.835 25.039 14.799 28.122 8.702 19.159 168.822 279.096 210.175 50.134 93.242 68.609 13.187 17.200 15.194 9.305 20.714 15.528 20.978 66.722 62.148 12.554 16.616 13.714 12.432 4.078 9.647 3.484 7.029 6.016 29.791 36.766 32.775

6.953 4.322 5.638 4.547 6.163 5.355 22.203 9.414 19.006 15.265 12.953 14.398 52.771 2.904 25.067 32.213 8.143 24.189 44.970 33.2566 39.646 14.737 26.277 19.786 14.637 5.730 10.526 39.047 40.025 39.423 49.849 65.229 56.440 11.800 20.103 15.952 7.260 14.702 11.319 24.308 70.709 66.609 1.570 4.812 5.293 3.330 9.519 5.580 10.643 6.683 7.814 20.947 20.056 21.884

6.924 4.322 5.623 4.764 5.993 5.378 29.462 26.632 28.755 17.333 12.825 15.642 52.220 38.269 46.020 37.316 28.219 34.284 45.259 34.144 40.207 14.717 26.277 19.775 14.479 5.605 10.392 39.070 40.025 39.547 44.375 110.067 72.529 5.944 11.538 8.741 4.328 5.558 4.998 28.547 73.532 69.034 1.550 4.801 5.293 2.990 8.251 5.014 15.523 10.046 11.601 21.459 25.653 24.873

14 14 28 6 6 12 6 2 8 5 3 8 5 4 9 8 4 12 6 5 11 9 7 16 7 6 13 5 3 8 4 3 7 7 7 14 5 6 11 1 9 10 4 3 7 8 5 13 2 5 7

a

For each system, N corresponds to the number of experimental data points available for comparison, and predictions were made using modUNIFAC and using MOSCED to parameterize NRTL, Wilson’s, and TK−Wilson equations as described in the text. For all cases, in the system name the first component corresponds to component 1. In the last line of the table, we calculate the average of AARE for all of the systems studied. 360



Article

overall AAE of y1 decreased to 0.0079 mol fractions. However, both sets of parameters performed poorly when predicting excess enthalpies. When excess enthalpy data was included in the parametrization, the predictions improved considerably. It follows that to improve the prediction of a particular property, experimental values of that property should be included in the parametrization. This comes at the slight expense of the other properties used in the parametrization. This is consistent with the results of the present study. We expect MOSCED to perform best when predicting limiting activity coefficients because MOSCED was parametrized using only limiting activity coefficient data, with only about 15% of the limiting activity coefficients from VLE.11 Nonetheless, the value of the limiting activity coefficients from VLE are sensitive to the excess Gibbs free energy model used to correlate the data. On the other hand, mod-UNIFAC21−26 is parametrized using limiting activity coefficients, VLE, liquid−liquid equilibrium, solid−liquid equilibrium, azeotropic data, excess enthalpies, and excess heat capacities. While MOSCED is in better agreement with experimental data when predicting limiting activity coefficients, mod-UNIFAC performs better when predicting azeotropic VLE.

Figure 6. A plot of the relative volatility (α1,2 = K1/K2) and the volatility of each component (K1 and K2) for the binary system cyclohexane(1)/1,2-dichloroethane(2) at 1.013 bar. The solid black line is experimental reference data51−56 and the dotted black lines are drawn for convenience as ±20% error in the reference data. The dashed lines correspond to predictions made using mod-UNIFAC21−26,36 and using MOSCED to parametrize NRTL (α = 0.3), TK−Wilson, and Wilson’s equation, as indicated. The horizontal gray line in the top plot corresponds to α1,2 = 1 and is useful for locating the azeotrope. The vertical gray line in the bottom two plots corresponds to x1 = 0.3 and x2 = 0.3, respectively, and is used to identify the dilute region.

■

CONCLUSION In the present study, we assessed the ability of MOSCED to predict solvation free energies (specifically hydration free energies) and Henry’s constants for solutes in water, and octanol/water partition coefficients. It was shown how Henry’s constants and partition coefficients may be directly related to solvation free energies, allowing us to probe directly MOSCED’s ability to accurately capture solute−solvent intermolecular interactions. Comparison of the computed hydration free energies to state-of-the-art molecular simulations demonstrated that MOSCED was able to predict hydration free energies in greater agreement with experiment. We propose that MOSCED and molecular simulation may be used as two complementary tools. MOSCED may be used to save computational time by performing solvation free energy calculations and computing phase-equilibria, computationally intensive tasks to perform using molecular simulation. Additionally, by relating solvation free energies to infinite dilution activity coefficients, we present the opportunity to use molecular simulation free energy calculations or electronic structure calculations in a continuum solvent to obtain missing MOSCED parameters, making MOSCED a truly predictive method. The computed hydration free energies, Henry’s constants, and partition coefficients were additionally compared to modUNIFAC. For all of these cases, the predictions made with MOSCED were found to be superior to that of mod-UNIFAC. Of the reference sets examined, MOSCED was able to model more compounds than mod-UNIFAC. However, at present, MOSCED is limited to 130 organic solvents and water while mod-UNIFAC is a group contribution method, so modUNIFAC would be able to model additional compounds outside this set. Be that as it may, the ability to use molecular simulation or electronic structure calculations in a continuum solvent to parametrize MOSCED presents an opportunity to expand the MOSCED parameter database. Moreover, MOSCED is parametrized for a wide range of chemical classes, making the development of a group contribution MOSCED an intriguing possibility.

related to limiting activity coefficients, which could be directly computed using MOSCED. This suggests that MOSCED was better able to predict limiting activity coefficients, while modUNIFAC is better able to predict VLE. When predicting VLE with MOSCED, we are dependent on an excess Gibbs free energy model to predict compositiondependent activity coefficients. Schreiber and Eckert 34 previously looked at the ability to use Wilson’s equation parametrized using infinite dilution activity coefficients to predict binary isothermal (pxy) VLE for 31 miscible pairs. Overall, they found the predictions were in good agreement with experiment, and that one may use relatively imprecise ∞ values of γ∞ 1 and γ2 and still obtain good VLE predictions. Correlating isothermal VLE data, they obtained an overall AAE in y1 of 0.0075 mol fractions using Wilson’s equation. When ∞ they instead used experimental values of γ∞ 1 and γ2 to parametrize Wilson’s equation, they obtained an overall AAE in y1 of 0.0086 mol fractions, an increase of 15%. Here, for isobaric (Txy) VLE, we expect the deviation to be larger as a ∞ result of requiring temperature-dependent values of γ∞ 1 and γ2 . This result is consistent with the findings of the first modUNIFAC paper.21 In that work, the authors first attempted to use only limiting activity coefficient data to parametrize modUNIFAC. The authors were then able to predict this set of limiting activity coefficients with an overall AAE of 5.3%. Using this set of parameters, they were able to to predict VLE with an overall AAE of y1 of 0.0096 mol fractions. If they additionally included VLE data in the parametrization, the overall AAE for limiting activity coefficients increased slightly to 5.6%, and the 361



Article

The excellent performance of MOSCED as compared to molecular simulation and mod-UNIFAC is additionally encouraging because during the 2005 parametrization of MOSCED, errors in the predictions in which water was a solvent were substantially larger than for any other case. Here we find that while MOSCED performs worst when water is the solvent, this is still improved as compared to molecular simulation and mod-UNIFAC. Additionally, MOSCED’s ability to predict octanol/water partition coefficients in good agreement with experiment suggests that MOSCED is able to accurately balance a compound’s hydrophobic and hydrophilic interactions. Comparison to mod-UNIFAC was additionally made to compute binary isobaric azeotropic vapor−liquid equilibrium. Azeotropic vapor−liquid equilibrium corresponds to significant deviations from ideality, and isobaric vapor−liquid equilibrium presents the additional challenge of predicting temperaturedependent activity coefficients. Overall for the examined reference set, we find that predictions made using modUNIFAC are in better agreement with the reference data than MOSCED, although the difference is not large. This performance is expected. While MOSCED is parametrized solely using infinite dilution activity coefficients, mod-UNIFAC additionally includes vapor−liquid equilibrium and azeotropic data in its parametrization. In addition to making quantitative predictions, the use of solubility parameter-based methods are advantageous as they provide qualitative insight into the underlying molecular-level driving forces for the behavior of a system. This allows us to both understand the observed behavior from a molecular level, in addition to allowing for intuitive solvent selection and design by comparing solubility parameters. As compared to regular solution theory and Hansen Solubility Parameters, the use of MOSCED is advantageous as it splits the association term into contributions due to hydrogen-bond acidity and basicity. We demonstrated how this subtle difference allows MOSCED to better model the molecular interactions of associating fluids. Additionally, to facilitate the adoption and use of MOSCED, an interactive software suite (MOSCED-CAD) has been developed. MOSCED-CAD can be used to calculate all of the properties described in the present study, in addition to binary interaction parameters for the Wilson and NRTL equations for use in a process simulator. A copy of the current user manual is provided in the Supporting Information which includes instructions for downloading from our personal Web site.27

■

■

from experiment and predicted using mod-UNIFAC and MOSCED; tabulation of MOSCED and Antoine parameters for all compounds for which MOSCED is parametrized (ZIP)

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel.: (513) 529-0784. Fax: (513) 529-0761. ORCID

Andrew S. Paluch: 0000-0002-2748-0783 Funding

Acknowledgment is made to the donors of the American Chemical Society Petroleum Research Fund (56896-UNI6) for support of this research. S.N.R. is additionally thankful for financial support from the Office of Research for Undergraduates at Miami University through the Undergraduate Summer Scholars (USS) program. E.M.S. is thankful for support from the College of Engineering and Computing at Miami University. Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Prausnitz, J. M.; Lichtenthaler, R. N.; de Azevedo, E. G. Molecular Thermodynamics of Fluid-Phase Equilibria, 2nd ed.; Prentice-Hall, Inc.: Englewood Cliffs, NJ, 1986. (2) Hildebrand, J. H.; Prausnitz, J. M.; Scott, R. L. Regular and Related Solutions; Van Nostrand Reinhold Company: New York, NY, 1970. (3) Abrams, D. S.; Prausnitz, J. M. Statistical thermodynamics of liquid mixtures: A new expression for the excess Gibbs energy of partly or completely miscible systems. AIChE J. 1975, 21, 116−128. (4) Blanks, R. F.; Prausnitz, J. M. Thermodynamics of Polymer Solubility in Polar and Nonpolar Systems. Ind. Eng. Chem. Fundam. 1964, 3, 1−8. (5) Hansen, C. M. The Universality of the Solubility Parameter. Ind. Eng. Chem. Prod. Res. Dev. 1969, 8, 2−11. (6) Tijssen, R.; Billiet, H. A. H.; Schoenmakers, P. J. Use of the solubility parameter for predicting selectivity and retention in chromatography. J. Chromatogr. A 1976, 122, 185−203. (7) Hansen Solubility Parameters. http://hansen-solubility.com/ (accessed July 13, 2017). (8) Tamir, A.; Apelblat, A.; Wagner, M. A New Device for Measuring Isothermal Vapor-Liquid Equilibria. Fluid Phase Equilib. 1981, 6, 237− 259. (9) Shephard, J. J.; Callear, S. K.; Imberti, S.; Evans, J. S. O.; Salzmann, C. G. Microstructures of negative and positive azeotropes. Phys. Chem. Chem. Phys. 2016, 18, 19227−19235. (10) Thomas, E. R.; Eckert, C. A. Prediction of limiting activity coefficients by a modified separation of cohesive energy density model and UNIFAC. Ind. Eng. Chem. Process Des. Dev. 1984, 23, 194−209. (11) Lazzaroni, M. J.; Bush, D.; Eckert, C. A.; Frank, T. C.; Gupta, S.; Olson, J. D. Revision of MOSCED Parameters and Extension to Solid Solubility Calculations. Ind. Eng. Chem. Res. 2005, 44, 4075−4083. (12) Ley, R. T.; Fuerst, G. B.; Redeker, B. N.; Paluch, A. S. Developing a Predictive Form of MOSCED for Nonelectrolyte Solids Using Molecular Simulation: Application to Acetanilide, Acetaminophen, and Phenacetin. Ind. Eng. Chem. Res. 2016, 55, 5415−5430. (13) Phifer, J. R.; Solomon, K. J.; Young, K. L.; Paluch, A. S. Computing MOSCED parameters of nonelectrolyte solids with electronic structure methods in SMD and SM8 continuum solvents. AIChE J. 2017, 63, 781−791. (14) Cox, C. E.; Phifer, J. R.; da Silva, L. F.; Nogueira, G. G.; Ley, R. T.; O’Loughlin, E. J.; Barbosa, A. K. P.; Rygelski, B. T.; Paluch, A. S.

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jced.7b00748. Additional discussion of the relationship of the limiting activity coefficient and solvation free energy; summary of error in predicted y1 and T; Txy and McCabe-Thiele (yx) diagrams (PDF) MOSCED-CAD user manual (PDF) Tabulation of hydration free energies from experiment and predicted using molecular simulation, mod-UNIFAC, and MOSCED; tabulation of Henry’s constants in water and octanol/water partition coefficients from experiment and predicted using mod-UNIFAC and MOSCED; tabulation of vapor−liquid equilibrium data 362



Article

Combining MOSCED with molecular simulation free energy calculations or electronic structure calculations to develop an efficient tool for solvent formulation and selection. J. Comput.-Aided Mol. Des. 2017, 31, 183−199. (15) Phifer, J. R.; Cox, C. E.; da Silva, L. F.; Nogueira, G. G.; Barbosa, A. K. P.; Ley, R. T.; Bozada, S. M.; O’Loughlin, E. J.; Paluch, A. S. Predicting the equilibrium solubility of solid polycyclic aromatic hydrocarbonds and dibenzothiophene using a combination of MOSCED plus molecular simulation or electronic structure calculations. Mol. Phys. 2017, 115, 1286−1300. (16) Diaz-Rodriguez, S.; Bozada, S. M.; Phifer, J. R.; Paluch, A. S. Predicting cyclohexane/water distribution coefficients for the SAMPL5 challenge with MOSCED and the SMD solvation model. J. Comput.-Aided Mol. Des. 2016, 30, 1007−1017. (17) Mobley, D. L.; Guthrie, J. P. FreeSolv: a database of experimental and calculated hydration free energies, with input files. J. Comput.-Aided Mol. Des. 2014, 28, 711−720. (18) Duarte Ramos Matos, G.; Kyu, D. Y.; Loeffler, H. H.; Chodera, J. D.; Shirts, M. R.; Mobley, D. L. Approaches for Calculating Solvation Free Energies and Enthalpies Demonstrated with an Update of the FreeSolv Database. J. Chem. Eng. Data 2017, 62, 1559−1569. (19) Moine, E.; Privat, R.; Sirjean, B.; Jaubert, J.-N. Estimation of Solvation Quantities from Experimental Thermodynamic Data: Development of the Comprehensive CompSol Databank for Pure and Mixed Solutes. J. Phys. Chem. Ref. Data 2017, 46, 033102. (20) FreeSolv: Experimental and Calculated Small Molecule Hydration Free Energies, version 0.51. https://github.com/MobleyLab/FreeSolv (accessed June 11, 2017). (21) Weidlich, U.; Gmehling, J. A Modified UNIFAC Model. 1. Prediction of VLE, hE, and γ∞. Ind. Eng. Chem. Res. 1987, 26, 1372− 1381. (22) Gmehling, J.; Li, J.; Schiller, M. A Modified UNIFAC Model. 2. Present Parameter Matrix and Results for Different Thermodynamic Properties. Ind. Eng. Chem. Res. 1993, 32, 178−193. (23) Gmehling, J.; Lohmann, J.; Jakob, A.; Li, J.; Joh, R. A Modified UNIFAC (Dortmund) Model. 3. Revision and Extension. Ind. Eng. Chem. Res. 1998, 37, 4876−4882. (24) Gmehling, J.; Wittig, R.; Lohmann, J.; Joh, R. A Modified UNIFAC (Dortmund) Model. 4. Revision and Extension. Ind. Eng. Chem. Res. 2002, 41, 1678−1688. (25) Jakob, A.; Grensemann, H.; Lohmann, J.; Gmehling, J. Further Development of Modified UNIFAC (Dortmund): Revision and Extension 5. Ind. Eng. Chem. Res. 2006, 45, 7924−7933. (26) Constantinescu, D.; Gmehling, J. Further Development of Modified UNIFAC (Dortmund): Revision and Extension 6. J. Chem. Eng. Data 2016, 61, 2738−2748. (27) A Better Solution Ahead: Welcome to MOSCED. https://sites. google.com/view/mosced, (accessed December 1, 2017). (28) Paluch, A. S.; Maginn, E. J. Predicting the Solubility of Solid Phenanthrene: A Combined Molecular Simulation and Group Contribution Approach. AIChE J. 2013, 59, 2647−2661. (29) Fuerst, G. B.; Ley, R. T.; Paluch, A. S. Calculating the Fugacity of Pure, Low Volatile Liquids via Molecular Simulation with Application to Acetanilide, Acetaminophen, and Phenacetin. Ind. Eng. Chem. Res. 2015, 54, 9027−9037. (30) Noroozi, J.; Paluch, A. S. Microscopic Structure and Solubility Predictions of Multifunctional Solids in Supercritical Carbon Dioxide: A Molecular Simulation Study. J. Phys. Chem. B 2017, 121, 1660− 1674. (31) Paluch, A. S.; Maginn, E. J. Efficient Estimation of the Equilibrium Solution-Phase Fugacity of Soluble Nonelectrolyte Solids in Binary Solvents by Molecular Simulation. Ind. Eng. Chem. Res. 2013, 52, 13743−13760. (32) Paluch, A. S.; Vitter, C. A.; Shah, J. K.; Maginn, E. J. A comparison of the solvation thermodynamics of amino acid analogues in water, 1-octanol and 1-n-alkyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide ionic liquids by molecular simulation. J. Chem. Phys. 2012, 137, 184505.

(33) Paluch, A. S.; Parameswaran, S.; Liu, S.; Kolavennu, A.; Mobley, D. L. Predicting the excess solubility of acetanilide, acetaminophen, phenacetin, benzocaine, and caffeine in binary water/ethanol mixtures via molecular simulation. J. Chem. Phys. 2015, 142, 044508. (34) Schreiber, L. B.; Eckert, C. A. Use of Infinite Dilution Activity Coefficients with Wilson’s Equation. Ind. Eng. Chem. Process Des. Dev. 1971, 10, 572−576. (35) Tsuboka, T.; Katayama, T. Modified Wilson Equation for Vapor-Liquid and Liquid-Liquid Equilibria. J. Chem. Eng. Jpn. 1975, 8, 181−187. (36) CHEMCAD version 7.1.0.9402; Chemstations, http://www. chemstations.com. (37) Yaws, C. L.; Narasimhan, P. K.; Gabbula, C. Yaws’ Handbook of Antoine Coefficients for Vapor Pressure (2nd Electronic ed.); Knovel, 2009. (38) For water, in MOSCED calculations we use a molar volume of v = 36 cm3/mol. For VLE calculations with Wilson’s and TK−Wilson equations, and octanol/water partition coefficient calculations, we use a molar volume of v = 18 cm3/mol. (39) Eaton, J. W.; Bateman, D.; Hauberg, S. GNU Octave version 3.0.1 manual: a high-level interactive language for numerical computations; CreateSpace Independent Publishing Platform, 2009; ISBN 1441413006. (40) When quantitatively comparing the predicted T and y1 to reference data, data at x1 = 0 and x1 = 1 were excluded (41) Jeffrey, G. A. An Introduction to Hydrogen Bonding; Oxford University Press, Inc.: New York, NY, 1997. (42) Góral, M.; Kolasińska, G.; Warycha, P. O. S. Vapor-Liquid Equilibria. II. The Ternary System Methanol-Chloroform-Acetone at 313.15 and 323.15 K. Fluid Phase Equilib. 1985, 23, 89−116. (43) Horstmann, S.; Fischer, K.; Gmehling, J. Isothermal VaporLiquid Equilibrium and Excess Enthalpy Data for the Binary System Water+Sulfolane and Methanol+N-Methyl-2-pyrrolidone. J. Chem. Eng. Data 2004, 49, 1449−1503. (44) Lazzaroni, M. J. Optimizing Solvent Selection for Separation and Reaction. Ph.D. Dissertation, Georgia Institute of Technology, 2004. (45) Yaws, C. L. Yaws’ Handbook of Properties for Aqueous Systems; Knovel, 2012. (46) Sangster, J. Octanol-Water Partition Coefficients of Simple Organic Compounds. J. Phys. Chem. Ref. Data 1989, 18, 1111−1227. (47) Campbell, J. R.; Luthy, R. G. Prediction of Aromatic Solute Partition Coefficients Using the UNIFAC Group Contribution Model. Environ. Sci. Technol. 1985, 19, 980−985. (48) Arce, A.; Blanco, A.; Souza, P.; Vidal, I. Liquid-Liquid Equilibria for Water + Methanol + 1-Octanol and Water + Ethanol + 1-Octanol at Various Temperatures. J. Chem. Eng. Data 1994, 39, 378−380. (49) Sherman, S. R.; Trampe, D. B.; Bush, D. M.; Schiller, M.; Eckert, C. A.; Dallas, A. J.; Li, J.; Carr, P. W. Compilation and Correlation of Limiting Activity Coefficients of Nonelectrolytes in Water. Ind. Eng. Chem. Res. 1996, 35, 1044−1058. (50) Confirmed using the “Group Assignment” tool online provided by DDBST at http://www.ddbst.com/unifacga.html, mod-UNIFAC groups exist for tetrahydrofuran (THF). However, these groups were not available in CHEMCAD. VLE calculations with mod-UNIFAC were therefore not performed for methanol/THF. (51) Gmehling, J., Onken, U., Eds. Vapor-Liquid Equilibrium Data Collection, Part 1: Aqueous-Organic Systems; DECHEMA: Frankfurt am Main, Germany, 1977. (52) Gmehling, J., Onken, U., Eds. Vapor-Liquid Equilibrium Data Collection, Part 2a: Organic Hydroxy Compounds: Alcohols; DECHEMA: Frankfurt am Main, Germany, 1977. (53) Gmehling, J., Onken, U., Arlt, W., Eds. Vapor-Liquid Equilibrium Data Collection, Part 2b: Organic Hydroxy Compounds: Alcohols and Phenol; DECHEMA: Frankfurt am Main, Germany, 1978. (54) Gmehling, J., Onken, U., Eds. Vapor-Liquid Equilibrium Data Collection, Part 2h: Alcohols: Ethanol and 1,2-Ethanediol, Supplement 6; DECHEMA: Frankfurt am Main, Germany, 2006. 363



Article

(55) Gmehling, J., Onken, U., Arlt, W., Eds. Vapor-Liquid Equilibrium Data Collection, Part 7: Aromatic Hydrocarbons; DECHEMA: Frankfurt am Main, Germany, 1980. (56) Stephen, K., Hildwein, H., Eds. Recommended Data of Selected Compounds and Binary Mixtures, Part 1 + 2; DECHEMA: Frankfurt am Main, Germany, 1987. (57) When more than one set of data was available for a binary pair, we chose the set with a pressure closest to 1 bar. If more than one set was available at this pressure, we chose the set for which DECHEMA was able to best correlate the data. For each binary, we compare to only a single set of reference data. The reference data used is tabulated along with the page number from DECHEMA in the Supporting Information accompanying this manuscript. (58) Mathias, P. M. Guidlines for the Analysis of Vapor-Liquid Equilibrium Data. J. Chem. Eng. Data 2017, 62, 2231−2233. (59) Mathias, P. M. Effect of VLE uncertainties on the design of separation sequences by distillation − Study of the bezene-chloroformacetone system. Fluid Phase Equilib. 2016, 408, 265−272. (60) Walas, S. M. Phase Equilibria in Chemical Engineering; Butterworth Publishers: Stoneham, MA, 1985.

364


Application of MOSCED To Predict Limiting Activity Coefficients

Recommend Documents