Extensive Accuracy Test of the Force-Field-Based Quasichemical

Article pubs.acs.org/jced

Extensive Accuracy Test of the Force-Field-Based Quasichemical Method PAC−MAC Augustinus J. M. Sweere,*,† Ruben Serral Gracia,‡ and Johannes G. E. M. Fraaije† †

Soft Matter Chemistry, Leiden University, Einsteinweg 55, 2333CC Leiden, The Netherlands Culgi B.V., P.O. Box 252, 2300AG Leiden, The Netherlands

‡

S Supporting Information *

ABSTRACT: We present a comprehensive evaluation of the recently developed pair configuration to molecular activity coefficient (PAC−MAC) method. PAC−MAC is a force-fieldbased quasichemical method for rapid calculation of binary phase diagrams. The accuracy of the method is tested by comparing the calculated excess mixing free energy with experimental data for 1092 binary mixtures. The root mean squared error (RMSE) is shown to be 0.15 kBT. Furthermore, a comparison with UNIFAC and molecular simulations is performed. UNIFAC shows a significantly higher accuracy (RMSE: 0.07 kBT), whereas molecular simulations lead to comparable results. The accuracy of both molecular simulations as well as PAC−MAC depends highly on the used force field. The binary parameters of UNIFAC are optimized using experimental miscibility data, whereas the force field parameters are not. Therefore, a better performance of UNIFAC is expected. A concise data set shows the capacity of PAC− MAC in predicting results obtained using Monte Carlo simulations.

■

105) adjacent molecular pair configurations, representing the possible observable configurations of two molecules within the first coordination shell of a mixture in the condensed phase. For each sampled pair, the intermolecular energy is calculated using classical force fields, by default the OPLS force field,9 and the portion of the surface occupied by the neighboring molecule is tracked. In the final step, an expression for the free energy, based on the quasichemical approximation of Guggenheim,8 is minimized with respect to the pair probabilities, subject to a constraint stating the molecular surface to be fully occupied by neighboring molecules. The constraint is added in order to obtain a proper packing of the molecules. Because it is assumed that a molecule in the condensed phase is completely surrounded, even if the molecule contains associative sites. From the expression of the free energy it is possible to obtain activity coefficients and corresponding miscibility properties. Sampling of the required 4 × 105 pair configurations in PAC−MAC takes about 5 min on a single core (Intel Xeon E52620) and the additional calculation time to obtain activity coefficients is about 20 s for each chosen mole fraction.7 Monte Carlo (MC) simulations often require about 109 MC steps to obtain a precise result for a single mole fraction,10 implying a great reduction of calculation time with the PAC−MAC model

INTRODUCTION A predictive model for the calculation of thermodynamic miscibility properties is the key to interpret, predict, or replace chemical experiments. Throughout the years several methods have been developed, with pros and cons for each of them. UNIFAC,1 COSMO-RS2,3 and SAFT4,5 are revered for their accuracy and calculation time. However, a common drawback of these methods is the lack of a three-dimensional molecular structure. Molecular dynamics,6 on the other hand, explicitly comprises the configurational information, but at the cost of extensive computations. A fast model, in which also the molecular conformation is taken into account, is the primary objective of the recently developed pair configurations to molecular activity coefficients (PAC−MAC) method.7 PAC−MAC is a force-field-based quasichemical method for the calculation of activity coefficients of multicomponent mixtures. A force field refers to simplified potential energy functions and corresponding parameters for the calculation of the potential energy of a molecular system. The quasichemical approximation of Guggenheim8 refers to an improvement of the regular solution model. The regular solution model assumes a random distribution of pairwise interactions between sites or molecules, whereas the quasichemical approximation contains a better description of the pair distribution by taking intermolecular pair energies into account. The operational procedure of PAC−MAC consists of three steps. The first step is defining the surface of an optimized molecular conformation for all components within the mixture. The second step involves the sampling of many 4m (with m = © XXXX American Chemical Society

Special Issue: Proceedings of PPEPPD 2016 Received: June 9, 2016 Accepted: September 12, 2016

A

DOI: 10.1021/acs.jced.6b00474 J. Chem. Eng. Data XXXX, XXX, XXX−XXX

Journal of Chemical & Engineering Data

Article

Fibonacci mapping technique.21 We generate a spiral over the convex hull of the molecular surface with a separation distance of 0.8 Å between the turns. A convex hull is equal to the solvent accessible surface using a spherical probe with an infinite radius.22 Subsequently, the centers of the segments are placed on the spiral separated by a distance of 0.8 Å from each other. Then, by using a repulsive force between the centers of the segments a more uniform distribution is obtained. Finally, the polygonal surface segments are obtained using a Voronoi diagram23 and the vertices of the segments are translated from the convex hull to the molecular surface. Figure 1 shows the segmentation of the molecular surface of acetonitrile (total area: 100.9 Å2) into 191 surface segments.

for the prediction of miscibility properties. A proof of concept of the PAC−MAC model was provided by comparing the obtained results with experimental data from the NIST Thermodynamic Data Engine database11 for a concise set of binary mixtures. In this paper, we extend the accuracy test of PAC−MAC by comparing the obtained excess Gibbs free energy of mixing with experimental data for a diverse data set containing 1092 binary mixtures. In the extension, we have found it expedient to introduce two modifications to the original PAC−MAC method. The first modification involves the introduction of vacuum surface elements to soften the constraint on molecular surface coverage. The other modification involves the Staverman−Guggenheim expression12 rewritten in PAC−MAC variables, to have a description of the combinatorial entropy of athermal mixtures. The obtained results are compared with UNIFAC13−15 and, for a few binary mixtures, also with molecular simulations.10,16,17 The derivation of the Staverman− Guggenheim expression, the derivation of the expression for the activity coefficients, and calculated data are included in the Supporting Information.

■

THEORETICAL BASIS The PAC−MAC method has been extensively discussed in a previous paper.7 The general outline of PAC−MAC is formulated in this section; two modifications to the original model are explained in detail. The first adjustment involves the introduction of a vacuum surface fraction to allow partial occupation around the molecule. The other adjustment involves a reformulation of the Staverman−Guggenheim12 expression. The method consists of three parts, explained in three subsections. The first subsection (“Segmentation of the Molecular Surface”) presents our definition of the molecular surface and the ensuing segmentation of the surface. The second section (“Sampling of Pair Configurations”) explains the sampling procedure of the pair configurations. In the last section (“Solving System of Equations”), the expression for the free energy is explained. Segmentation of the Molecular Surface. First, the molecular surfaces of all components in the mixture are required to perform a PAC−MAC calculation for obtaining miscibility properties. The molecular surface is defined by the exterior of spheres centered on the atomic nuclei.18 The radius of the sphere around atom k is given by R k = 0.66 × σkk

Figure 1. Acetonitrile molecule. Left: Ball-and-stick representation. Right: Molecular surface in PAC−MAC, divided into 191 surface segments of approximately 0.5 Å2.

Sampling of Pair Configurations. In the PAC−MAC method m pair configurations are sampled, with m = 105 by default, between all possible combinations of adjacent molecules. The sampling procedure consists of 6 steps: (1) The vertices of the polygonal surface segments are translated from the molecular surface to the convex hull of the surface to avoid indentation. (2) A random point on the convex hull surface is chosen on both molecules. (3) The molecules are rotated around their geometrical center in order to obtain the normal vector of both points on the convex hull surface pointing in opposite direction. (4) The molecules are translated in order to have both points on the convex hull surface on the same location. (5) Both molecules obtain a random rotation around the normal vector. (6) Both molecules are randomly translated along the normal vector according to a uniform distribution in an interval controlled by atom type-independent distance scale parameters. See Figure 2 for an illustration of the pair sampling procedure. For each sampled pair configuration two calculations are performed. The first calculation is the intermolecular energy wi using the OPLS all-atom force field.9 Within the OPLS force field the total intermolecular energy consists of Coulomb and Lennard-Jones interactions. Hydrogen bonding effects are included within the point charges of the atoms. The second calculation involves tracking of the surface segments that are blocked (we also refer to such segments as “occupied”) by the neighboring molecule. A surface segment is blocked if its normal vector, starting from the center of the segment and perpendicular to the convex hull surface, pierces the molecular surface of the neighboring molecule. Figure 3 shows the

(1)

where σkk is the Lennard-Jones distance parameter and the factor 0.66 (referred to as sseg in the previous paper) is an atomtype independent scaling coefficient obtained by fitting the coordination numbers of 12 selected solvent molecules calculated using molecular dynamics (MD) simulations.7 The molecular surface is divided into polygonal surface segments (mostly pentagon- and hexagon-shaped) with an area AIj that is approximately equal to 0.5 Å2. The chosen surface segment size is a trade-off between a fast PAC−MAC calculation (large size) and a high resolution of the molecule (small size). In a previous article, we used 200 surface segments independent of the size of the molecule.7 A segment size of 0.5 Å2 results in a comparable resolution for small molecules like acetonitrile (191 surface segments) and ethanol (206 surface segments). Various algorithms are developed to create a uniform distribution of segments over the molecular surface.3,19,20 Within PAC−MAC the segmentation is based on the spherical B



Article

where z is the coordination number. We also introduce in this paper the term vacuum fraction: the fraction of the molecular surface that is not in contact with a neighboring molecule. Because of steric effects the molecular surface will never be fully occupied by surrounding molecules in the first coordination shell. Even in the fcc close-packed structure of equal spheres 20% of the area is not covered by neighboring spheres, see the illustration in Figure 4.

Figure 2. Pair sampling procedure. (1) Convex hull around the molecular surface. (2) Choosing a random point on the convex hull surface. (3) Rotation of molecules to obtain two normal vectors in opposite direction. (4) Translation of molecules to have both points on the same location. (5) Rotation of molecules around their normal vector. (6) Random translation in interval. Colors of surface segments refer to the closest atomic centers: hydrogen (white), carbon (gray), nitrogen (blue), or oxygen (red).

intermolecular energy and occupied surface segments of the pair configuration obtained in Figure 2.

Figure 4. Fcc close-packing of equal spheres. The green spheres cover 80.4% of the surface of the central red sphere.

We assume that no intermolecular energy is involved in the vacuum-contribution. This leads to a free energy contribution that is purely entropic F vac = NkBT

F int 1 = z NkBT 2

∑

xI ln(xI ) (3)

⎛φ ⎞ ∑ xI ln⎜ I ⎟ + ⎝ xI ⎠ I∈A ,B

x −φ ⎛y⎞ ∑ yI I I ln⎜⎜ I ⎟⎟ yI − φI ⎝ φI ⎠ I∈A ,B

1 z 2

⎛ IJ ⎞ x ·m ⎟ xiIJ ln⎜⎜ i ⎟ ⎝ yI yJ ⎠ I∈A ,B J∈A ,B i=1

I∈A ,B J∈A ,B i=1

wiIJ kBT

m

∑ ∑ xiIJoijIJ + J∈A ,B i=1

1 z 2

(7)

m

∑ ∑ xiJIoijJI + xIxjvac , I = xI ∀ J∈A ,B i=1

j ∈ [1..LI ], ∀ I ∈ A , B (4)

(8)

So the introduction of a vacuum fraction weakened the former strict constraint stating the molecule to be completely surrounded. The free energy, eq 2, is minimized by optimizing variables xIJi , xvac j,I and z subject to the constraint, eq 8. Finally, the expression for the activity coefficient of molecule I ∈ A, B is obtained using the “intercept rule”

with φI and yI being respectively the volume fraction and coordination fraction of component I in the mixture. Equation 4 is an analytically derived expression to describe the total number of possible ways to arrange molecules on a lattice, taking size and shape into account.24 The quasichemical combinatorial term is8 F comb 1 = z NkBT 2

m

∑ ∑ ∑ xiIJ

In which is the intermolecular energy of pair configuration i between molecules I and J calculated using the OPLS force field. In the original PAC−MAC model the vacuum terms are absent. Because of the vacuum inclusion the surface coverage constraint is now that each surface segment is covered by neighboring molecules or by the segmental vacuum fraction

(2)

where xI is the mole fraction of component I in the mixture of molecules A and B. FSG is the Staverman−Guggenheim correction (see Supporting Information) F SG = NkBT

(6)

w iIJ

Fid is the entropy of mixing of an ideal solution:

I∈A ,B

I∈A ,B j=1

⎞ ⎛ xjvac ,I ⎟ ⎜ ln vac ⟨Ai ⟩I ⎝ xI ⎠ AjI

with being the vacuum surface fraction of surface segment j, I xvac I being the average vacuum surface fraction of molecule I, Aj being the surface area of surface segment j and ⟨Ai⟩I being the average occupied surface area in an interaction on molecule I. Fint is the total intermolecular energy of the binary mixture

Solving System of Equations. The probability xIJi to obtain the sampled pair configuration i between molecules I and J is calculated by minimizing an expression for the free energy. The expression for the free energy is given by

F id = NkBT

∑ ∑

xIxjvac ,I

xvac j,I

Figure 3. Calculation of the intermolecular energy and blocked surface segments (black colored) for a sampled pair configuration.

F mix = F id + F SG + F comb + F vac + F int

LI

m

F ex ln(γI ) = + (1 − xI ) NkBT

∑ ∑ ∑

(5) C

F ex NkBT

( )

∂

∂xI

(9)



Article

Figure 5. Distribution of the 1092 binary mixtures in the experimental data set. Molecules are clustered in 9 categories.

We refer to the Supporting Information for more details concerning the derivation of the activity coefficients, eq 9, in PAC−MAC.

■

RESULTS AND DISCUSSION A quantitative comparison of the original PAC−MAC model is performed by comparing the results for the excess free energy with experimental data from the NIST database.25 The experimental data set consists of 1092 vapor−liquid equilibrium diagrams of binary mixtures at constant temperature (Pxy diagrams) and contains 166 different molecules. The distribution of the combinations of molecules is depicted in Figure 5. The 2184 components in the experimental data set are clustered in 9 categories. The 12 ketones, 4 aldehydes, 16 ethers, and 11 esters in the set are grouped in the category “molecules containing oxygen (excluding alcohols)”. The 10 amines, 3 nitriles, and 3 nitrogen-containing aromatic compounds are grouped in the category “molecules containing nitrogen”. Three examples of the 22 molecules with different or multiple functional groups are DMSO, 2-aminoethanol, and nitrobenzene. The temperatures in the experimental data set are on average 336.3 K and 95% of the binary mixtures have a temperature between 280 and 400 K. The temperature distribution of the 1092 Pxy diagrams is presented by a histogram in Figure 6. Using the extended Raoult’s Law, assuming gaseous ideality, the activity coefficient of component I is calculated from the experimental Pxy diagram26 γI =

Figure 6. Temperature distribution of the experimental data set of 1092 Pxy diagrams. Average temperature: 336.3 K. Lowest temperature: 253.2 K (acetone−hexane). Highest temperature: 523.1 K (quinoline−meta-cresol).

The PAC−MAC model is tested by comparing the calculated ex liq excess free energy, Fex PAC−MAC, with Fexp at a mole fraction of xA liq = xB = 0.5 for all 1092 data points. If the activity coefficients γA liq and γB are not available at xliq A = xB = 0.5, then an estimation is performed by linear interpolation between the nearest neighboring activity coefficients. The accuracy is represented by the correlation coefficient and the root mean squared error (RMSE) given by

xIgasP(xIliq) xIliqP(xIliq = 1)

xIgas

(10)

RMSE =

xIliq

With and representing the mole fraction of component I in respectively the vapor and liquid phase, P(xliq I ) is the vapor pressure of the binary mixture at mole fraction xliq I . Subsequently, the experimental excess free energy of a binary mixture containing molecules A and B is calculated using ex Fexp

kBT

= xAliq ln(γA) + xBliq ln(γB)

1 n

n ex ex 2 ∑ (FPAC −MAC, k − Fexp, k) k=1

(12)

With n being the number of data points and k being the index of a data point. The experimental values for the excess free energy of mixing are included in the Supporting Information. ex A scatter plot of Fex PAC−MAC versus Fexp for 1092 binary mixtures using the original vacuum-free PAC model is shown in Figure 7. The data set is divided into 5 different categories to

(11) D



Article

Figure 7. Scatter plot of excess free energies of mixing, original PAC− MAC versus experimental.

Figure 8. Root mean squared error of the modified PAC−MAC method using various combinations of sseg and xvac I , (neglecting 10% of the most positive and negative deviating outliers).

obtain more insight in the plot. The 77 mixtures containing water are grouped together in the category “Water−X”. All 474 remaining mixtures containing one or more hydrogen bonding donors (alcohols or amines) are represented by the category “HB donor−X”. Mixtures without any hydrogen bonding involved are divided in mixtures in which both molecules contain heteroatoms (“Heteroatom−Heteroatom”, 130 data points), mixtures in which only one molecule contains heteroatoms (“Heteroatom−Hydrocarbon”, 285 data points) and mixtures of hydrocarbons (“Hydrocarbon−Hydrocarbon”, 126 data points). Figure 7 shows reasonable correlation between the calculated and experimental excess free energy. However, the plot contains many outliers significantly reducing the correlation coefficient to 0.558 and increasing the RMSE to 0.273 kBT. The outliers are primarily attributable to mixtures containing water. Omission of the 77 binary mixtures containing water results in a reduction of the RMSE to 0.177 kBT. As predicted by the previous paper7 the PAC−MAC method poorly describes mixtures containing water, with the used s-parameters. One of the reasons for the obtained outliers is the dominant constraint that each surface segment is fully occupied. With the introduction of a vacuum surface fraction the constraint is weakened: partial occupation of the surface segments is made possible. A modified PAC−MAC model also requires modified general parameters. The values for smin and smax remain unaffected to respectively 0.36 and 0.72. The value for size scaling coefficient sseg and the value of the parameter for the vacuum fraction xvac I are optimized by minimizing the RMSE. The 10% most deviating PAC−MAC results, in positive and negative direction, are omitted in the calculation of the RMSE to avoid a misleading contribution of the outliers in the optimization of sseg and xvac I ; a calculation technique often used in jury sports.27 Figure 8 shows the obtained RMSE of PAC− MAC for 133 combinations of sseg and xvac I . Note that the observed minimum RMSE of 0.069 kBT does not represent the actual accuracy of PAC−MAC because the most deviating results are omitted. The original PAC−MAC method uses sseg = 0.66 and xvac I = 0.0, so Figure 8 shows a potential reduction of the RMSE by using a vacuum surface fraction within the PAC−MAC method. Figure 8 also shows that the value of xvac I , minimizing the RMSE, increases for a decreasing value of sseg. Decreasing sseg

reduces the occupied surface area and therefore implies an increasing xvac I to preserve an equal coordination number. The global minimum RMSE is obtained for sseg = 0.62 and xvac I = 0.6. A scatter plot presenting the accuracy of the modified PAC− MAC model using the optimized parameters is given in Figure 9.

Figure 9. Scatter plot of excess free energies of mixing, vacuum modified PAC−MAC versus experimental.

As shown in Figure 9 we find that the introduction of a vacuum fraction leads to a significant reduction of the outliers, with a reduction of the root mean squared error to 0.153 kBT and an increase of the correlation coefficient to 0.695. On the basis of the PAC−MAC calculation, not all binary mixtures are thermodynamically stable at equimolar mole fraction. According to Flory−Huggins lattice theory, a calculated Fex above 0.5 kBT always results in phase separation for a binary mixture of two equal-sized components. Within the quasichemical PAC−MAC model the critical value of Fex depends on the components in the mixture. A homogeneous mixture is unconditionally stable over the whole fraction range if E



Article

Figure 10. Miscibility prediction of the PAC−MAC model. (Left) Original PAC−MAC model without an additional vacuum surface fraction. (Right) Modified PAC−MAC model with an additional vacuum surface fraction. The dashed line represents the critical value for Fex in Flory− Huggins theory for a binary mixture of equal-sized molecules.

∂ 2F mix ≥ 0 ∀ xA ∈ [0.0, 1.0] ∂xA2

(13)

Therefore, we define a binary mixture to be miscible if eq 13 is fulfilled and immiscible otherwise. Figure 10 shows whether or not phase separation is predicted by PAC−MAC: It is assumed that all 1092 binary mixtures of the used experimental data set are thermodynamically stable. According to Figure 10, the amount of incorrectly predicted immiscible mixtures is significantly reduced by the addition of a vacuum surface fraction in the PAC−MAC model: the original model wrongly predicts 17.6% to be immiscible, whereas the modified model predicts 5.5% to be immiscible. Using the group-contribution method UNIFAC, it is also possible to obtain excess free energies of mixing.13−15 The UNIFAC method combines a thermodynamic model based on the quasi chemical theory (UNIQUAC28) with groupinteraction parameters optimized by fitting experimental data of a large number of multicomponent mixtures. Various parameter sets have been developed for the UNIFAC model. For comparison with experimental data, we used a public list of 56 groups and corresponding parameters extracted from the AIM Aerosol Thermodynamics Model Web site.29 A scatter plot of the results obtained with UNIFAC is given in Figure 11. The root mean squared error and correlation coefficient obtained with the UNIFAC model for the same data set of binary mixtures is given by 0.07 kBT and 0.932, respectively. So, despite the increasing correlation coefficient with the introduction of a vacuum surface fraction in PAC−MAC, the accuracy is still much less than the accuracy of UNIFAC. The calculation with UNIFAC is very fast: a few seconds on a single core for all points. However, UNIFAC is a group contribution method containing over one hundred parameters optimized using experimental miscibility data. PAC−MAC contains only two parameters optimized using experimental free energies of mixing: the nonbonded OPLS force field parameters are derived from thermodynamic properties of pure liquids.9 More parameters optimized based on miscibility data will result in increasing accuracy of PAC−MAC. UNIFAC cannot be used if a molecule contains undefined groups. Therefore, the scatter plot in Figure 11 contains 1042 data points whereas our data set contains 1092 binary mixtures. We realize that because PAC−MAC is a force-field-based method, its accuracy is limited to the accuracy of the force

Figure 11. Scatter plot of excess free energies of mixing, UNIFAC versus experimental.

fields. PAC−MAC in principle predicts results of molecular dynamics (MD) or Monte Carlo (MC) simulations and therefore PAC−MAC can never achieve a higher accuracy than MD or MC simulations. Hence, further improvement in PAC−MAC could be possible if we consider reparameterization of the underlying force fields directly. A quantitative comparison of MD with experimental data is required to test the accuracy of the used force fields; a quantitative comparison of MD with PAC−MAC proves the plausibility of the assumptions made within the model. But here, we encounter a considerable challenge, which in fact has been known for a long time in the molecular dynamics community. A quantitative comparison with molecular simulations is extremely difficult simply because of a lack of tabulated data, which immediately goes back to the great computational expense of simulations. Fortunately, for our purpose, Jedlovszky et al. did perform various MC simulations for the calculation of mixing Helmholtz free energies of binary mixtures through thermodynamic integration.10,16,17,30 Albeit, the data set is still small, the results are very useful for comparison with PAC−MAC. The comparison for a binary mixture of acetone and methanol is shown in Figure 12. Mixing free energies are obtained at 298 K for three different force fields. Corresponding experimental data is provided by Wilsak, Campbell, and Thodos.31 F



Article

Figure 12. Free energy of mixing for a binary mixture of acetone−methanol at 298 K. Calculated data obtained by MC (left graph) and PAC−MAC (right graph). Three different force fields are used: OPLS acetone−OPLS methanol (blue), TraPPE acetone−TraPPE methanol (green) and TraPPE-mod acetone−TraPPE methanol (red). MC results are drawn from tabulated data kindly provided by Jedlovszky.16 Original experimental data is obtained by Wilsak et al.31

Figure 13. Free energy of mixing for a binary mixture of DMSO−water at 298 K. Calculated data obtained by MC (left graph) and PAC−MAC (right graph). Water is described by the SPC force field. DMSO is described by eight different force fields: P1 (red), P2 (green), RS (blue), VLL (light blue), LMPvG (pink), VB (yellow), GROMOS (olive), and OPLS (dark blue). MC data provided by Jedlovszky et al.10 Experimental data is obtained by Lam and Benoit.32

As shown in Figure 12, neither MC nor PAC−MAC fit the experimental data perfectly. However, a comparison of MC with PAC−MAC reflects the similarities of both methods. Sorting the force fields by the minimum of the mixing free energy results in the same sequence. Moreover, the maximum difference between both methods is only 0.32 kJ/mol (= 0.13 kBT). A similar comparison between MC and PAC−MAC is performed for a binary mixture of DMSO−water at 298 K. The results are shown in Figure 13. Experimental data is provided by Lam and Benoit.32 PAC−MAC is inaccurate for mixtures containing water. Figure 13 shows that, depending on the force field, the mixing free energies obtained by MC and PAC−MAC can deviate a lot from experimental data: up to 2.3 kJ/mol (= 0.91 kBT). For DMSO−water, the difference in the obtained mixing free energy between both methods is not as close as for acetone− methanol. Still, sorting the force fields by the minimum of the mixing free energy results in a comparable sequence. Plotting Fmix calculated by PAC−MAC versus Fmix obtained with molecular simulations results in the scatter plot given in Figure 14. The scatter plot compares the mixing free energy of MC and PAC−MAC for three different equimolar (xA = xB = 0.5) binary mixtures with various combinations of force fields. The three molecular simulations methanol−acetone16 (shown

Figure 14. Scatter plot of free energies of mixing, PAC−MAC versus MC. All molecular simulations are performed by Jedlovszky et al.10,16,17

in Figure 12), DMSO−water10 (partly shown in Figure 13) and acetone−water17 (not shown) are all performed by Jedlovszky et al. G



■

Figure 14 shows a high correlation between PAC−MAC and MC in the calculation of mixing free energies. Also, the difficulties of PAC−MAC in handling water are revealed; the data obtained by PAC−MAC for the two mixtures containing water appear to be biased, molecule specific s-parameters are suggested as a solution to overcome this problem.7 A quantitative comparison with more binary mixtures is planned for future work. On the basis of Figures 12, 13, and 14, it can be concluded that the force field has a big influence on the accuracy of PAC− MAC. When the fast calculation time of PAC−MAC is taken into account, it may open doors in the optimization of force fields. These force fields will not only improve the accuracy of PAC−MAC but, most probably, will also improve the accuracy of the underlying MD or MC simulations. Then, the use of the optimized force fields will greatly exceed the calculation of only mixing free energies.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Funding

We are very grateful to SABIC and, in particular, K. Remerie for providing financial support for our research. Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS We would like to thank J. v. Male for useful discussions. REFERENCES

(1) Fredenslund, A.; Jones, R. L.; Prausnitz, J. M. GroupContribution Estimation of Activity Coefficients in Nonideal Liquid Mixtures. AIChE J. 1975, 21, 1086−1099. (2) Klamt, A. Conductor-like Screening Model for Real Solvents: A New Approach to the Quantitative Calculation of Solvation Phenomena. J. Phys. Chem. 1995, 99, 2224−2235. (3) Klamt, A. COSMO-RS: From Quantum Chemistry to Fluid Phase Thermodynamics and Drug Design, 1st ed.; Elsevier: Amsterdam, 2005. (4) Chapman, W. G.; Gubbins, K. E.; Jackson, G.; Radosz, M. SAFT: Equation-of-State Solution Model for Associating Fluids. Fluid Phase Equilib. 1989, 52, 31−38. (5) Gross, J.; Sadowski, G. Perturbed-Chain SAFT: An Equation of State Based on a Perturbation Theory for Chain Molecules. Ind. Eng. Chem. Res. 2001, 40, 1244−1260. (6) Alder, B. J.; Wainwright, T. E. Studies in Molecular Dynamics. I. General Method. J. Chem. Phys. 1959, 31, 459−466. (7) Sweere, A. J. M.; Fraaije, J. G. E. M. Force-Field Based Quasichemical Method for Rapid Evaluation of Binary Phase Diagrams. J. Phys. Chem. B 2015, 119, 14200−14209. (8) Guggenheim, E. A. Statistical Thermodynamics of Mixtures with Non-Zero Energies of Mixing. Proc. R. Soc. London, Ser. A 1944, 183, 213−227. (9) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225−11236. (10) Idrissi, A.; Marekha, B.; Barj, M.; Jedlovszky, P. Thermodynamics of Mixing Water with Dimethyl Sulfoxide, as Seen from Computer Simulations. J. Phys. Chem. B 2014, 118, 8724−8733. (11) Frenkel, M.; Chirico, R. D.; Diky, V.; Muzny, C. D.; Kazakov, A. F.; Magee, J. W.; Abdulagatov, I. M.; Kroenlein, K.; Diaz-Tovar, C. A.; Kang, J. W.; Gani, R. ThermoData Engine, NIST Standard Reference Database #103b, Version 7.0 (Pure Compounds, Binary Mixtures, Ternary Mixtures, and Chemical Reactions); National Institute of Standards and Technology: Gaithersburg, MD, 2011. (12) Staverman, A. J. The Entropy of High Polymer Solutions. Generalization of Formulae. Rec. Trav. Chim. Pays-Bas 1950, 69, 163− 174. (13) Hansen, H. K.; Rasmussen, P.; Fredenslund, A.; Schiller, M.; Gmehling, J. Vapor-Liquid Equilibria by UNIFAC Group Contribution. 5. Revision and Extension. Ind. Eng. Chem. Res. 1991, 30, 2352− 2355. (14) Wittig, R.; Lohmann, J.; Gmehling, J. Vapor-Liquid Equilibria by UNIFAC Group Contribution. 6. Revision and Extension. Ind. Eng. Chem. Res. 2003, 42, 183−188. (15) Balslev, K.; Abildskov, J. UNIFAC Parameters for Four New Groups. Ind. Eng. Chem. Res. 2002, 41, 2047−2057. (16) Idrissi, A.; Polok, K.; Barj, M.; Marekha, B.; Kiselev, M.; Jedlovszky, P. Free Energy of Mixing of Acetone and Methanol: A Computer Simulation Investigation. J. Phys. Chem. B 2013, 117, 16157−16164. (17) Pinke, A.; Jedlovszky, P. Modeling of Mixing Acetone and Water: How Can Their Full Miscibility Be Reproduced in Computer Simulations? J. Phys. Chem. B 2012, 116, 5977−5984.

■

CONCLUSION In this paper an extensive test of the accuracy of the PAC− MAC model is presented. A diverse set of 1092 binary mixtures is used to compare the calculated excess Gibbs free energy of mixing with experimental data. The predictive capacity of PAC−MAC is proven; however, many outliers also are obtained, significantly increasing the root-mean-square error (RMSE) to 0.27 kBT. The outliers are mainly attributable to mixtures containing water. Two adjustments to the model are introduced. The first adjustment involves partial occupation of the molecular surface and takes into account holes between the molecules within the first coordination shell. The second adjustment involves rewriting of the Staverman−Guggenheim expression to PAC−MAC variables and inclusion of this expression into the formula for the free energy. A significant reduction of the outliers is obtained with these adjustments, resulting in a correlation coefficient of 0.70 and a RMSE of 0.15 kBT. The PAC−MAC method shows a lower accuracy than UNIFAC. PAC−MAC in principle predicts molecular dynamics results. As a consequence, the accuracy of PAC−MAC highly depends on the accuracy of the force field. Mixing free energies calculated using PAC−MAC are shown to be comparable with the mixing free energies obtained by MC. The computational time for a PAC−MAC calculation is in the order of minutes. Therefore, PAC−MAC might be used as a quick method to optimize force field parameters for MD or MC simulations based on experimental miscibility data. Then, the calculation of mixing free energies is just one example of the whole range of possibilities with these optimized force fields.

■

Article

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jced.6b00474. Appendix A: Staverman−Guggenheim expression in PAC−MAC. Analytic derivation of the Staverman−Guggenheim correction term, eq 4, used in PAC−MAC. Appendix B: Quasichemical derivation of PAC−MAC. Analytic derivation of the expression for the activity coefficients, eq 9, in PAC−MAC. Appendix C: Experimental and calculated excess free energies of mixing. Table containing the experimental and calculated excess free energies of mixing. (PDF) H



Article

(18) Bondi, A. Van der Waals Volumes and Radii. J. Phys. Chem. 1964, 68, 441−451. (19) Silla, E.; Tuñoń , I.; Pascual-Ahuir, J. L. GEPOL: An Improved Description of Molecular Surfaces II. Computing the Molecular Area and Volume. J. Comput. Chem. 1991, 12, 1077−1088. (20) Zauhar, R. J. SMART: A Solvent-Accessible Triangulated Surface Generator for Molecular Graphics and Boundary Element Applications. J. Comput.-Aided Mol. Des. 1995, 9, 149−159. (21) Keinert, B.; Innmann, M.; Sänger, M.; Stamminger, M. Spherical Fibonacci Mapping. ACM Trans. Graph. 2015, 34, 193. (22) Lee, B.; Richards, F. M. The Interpretation of Protein Structures: Estimation of Static Accessibility. J. Mol. Biol. 1971, 55, 379−IN4. (23) Voronoi, G. Nouvelles Applications des Paramètres Continus à la Théorie des Formes Quadratiques. Premier Mémoire. Sur Quelques Propriétés des Formes Quadratiques Positives Parfaites. J. Reine Angew. Math. 1908, 133, 97−178. (24) Bronneberg, R.; Pfennig, A. MOQUAC, a New Expression for the Excess Gibbs Free Energy Based on Molecular Orientations; Hochschulbibliothek der Rheinisch-Westfälischen Technischen Hochschule Aachen: Aachen, Germany, 2012. (25) Flory, P. J. Thermodynamics of High Polymer Solutions. J. Chem. Phys. 1942, 10, 51−61. (26) Gmehling, J.; Kolbe, B.; Kleiber, M.; Rarey, J. Chemical Themodynamics for Process Simulation; Wiley-VCH: Weinheim, Germany, 2012. (27) Grandi, B.; Gueisbuhler, A. F. Regulations for the Judges’ Evaluation Programme and its’ Application; Fédération Internationale de Gymnastique: Lausanne, Switzerland, 2015; pp 1−5. (28) 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. (29) AIM Aerosol Thermodynamics Model. UNIFAC Structural Groups and Parameters. http://www.aim.env.uea.ac.uk/aim/info/ UNIFACgroups.html (accessed September 1, 2016). (30) Darvas, M.; Jedlovszky, P.; Jancsó, G. Free Energy of Mixing of Pyridine and Its Methyl-Substituted Derivatives with Water, As Seen from Computer Simulations. J. Phys. Chem. B 2009, 113, 7615−7620. (31) Wilsak, R. A.; Campbell, S. W.; Thodos, G. Vapor-Liquid Equilibrium Measurements for the Methanol-Acetone System at 372.8, 397.7 and 422.6 K. Fluid Phase Equilib. 1986, 28, 13−37. (32) Lam, S. Y.; Benoit, R. L. Some Thermodynamic Properties of the Dimethylsulfoxide-Water and Propylene Carbonate-Water Systems at 25 °C. Can. J. Chem. 1974, 52, 718−722.

I


Extensive Accuracy Test of the Force-Field-Based Quasichemical

Recommend Documents