Force-Field Based Quasi-Chemical Method for Rapid Evaluation of

Article pubs.acs.org/JPCB

Force-Field Based Quasi-Chemical Method for Rapid Evaluation of Binary Phase Diagrams Augustinus J. M. Sweere and Johannes G. E. M. Fraaije* Soft Matter Chemistry, Leiden University, Einsteinweg 55, 2333CC Leiden, The Netherlands ABSTRACT: We present the Pair Configurations to Molecular Activity Coefficients (PAC-MAC) method. The method is based on the pair sampling technique of Blanco (Fan, C. F.; Olafson, B. D.; Blanco, M.; Hsu, S. L. Application of Molecular Simulation to Derive Phase Diagrams of Binary Mixtures. Macromolecules 1992, 25, 3667−3676) with an extension that takes the packing of the molecules into account by a free energy model. The intermolecular energy is calculated using classical force fields. PAC-MAC is able to predict activity coefficients and corresponding vapor−liquid equilibrium diagrams at least 4 orders of magnitude faster than molecular simulations. The accuracy of the PAC-MAC method is tested by comparing the results with experimental data and with the results of the COSMO-SAC model (Lin, S.-T.; Sandler, S. I. A Priori Phase Equilibrium Prediction from a Segment Contribution Solvation Model. Ind. Eng. Chem. Res. 2002, 41, 899−913). PAC-MAC (using the OPLS-aa force field) is shown to be comparable in accuracy to COSMO-SAC, at the considerable advantage that PAC-MAC in principle does not require quantum calculation, provided proper force fields to be available.

■

INTRODUCTION Activity coefficients of multicomponent fluid mixtures provide a useful source of information for the calculation of various chemical properties, such as vapor−liquid equilibrium data and Flory−Huggins interaction parameters. Experimental databases are one way to obtain these properties, for example, the NIST database1 (the source of the data we will use for validation) which contains data for more than 42 500 binary mixtures. However, for cases where experimental data are not available, one has to find other ways to predict the thermodynamic behavior of fluid mixtures. The statistical associating fluid theory (SAFT)2,3 is a semiempirical equation-of-state model predicting phase equilibria accurately in a broad range of conditions.4 However, to obtain accurate results, SAFT requires molecular-dependent parameters, typically fitted to vapor pressure and density data of the pure components. In that case, experimental data is still necessary. An older contribution to the prediction of miscibility properties is made by Fowler and Guggenheim; their quasichemical theory forms the basis for various thermodynamic models, 5 as presented in the textbook of Prausnitz, Lichtenthaler, and de Azevedo.6 We address three variants: UNIFAC,7 COSMO (COSMO-RS8,9 and COSMO-SAC10), and the novel F-SAC.11 In the UNIFAC model, a molecule is completely split in independent functional groups. As with all group contribution methods, UNIFAC relies on two factors: First, on the definition of the functional groups, which is not unique.12,13 Second, on the data set of experimental phase diagrams used for the parametrization of the functional groups. In the COSMO variants, a molecule is generally not described as a collection of independent functional groups, but as a © 2015 American Chemical Society

collection of independent surface panels. In a mixture, the surface panels will form pairwise interactions with each other. A single quantum mechanical calculation is required to obtain the charge density of each surface panel of one molecule in a perfect conductor. With the distribution of the surface area as a function of the charge density, the so-called sigma profile, the COSMO-based methods are able to predict activity coefficients with a fraction of the amount of parameters used in UNIFAC.14 The third model, F-SAC, is basically a combination of COSMO-SAC and UNIFAC. Identical to COSMO-SAC, the F-SAC model uses sigma profiles to calculate the energetic interaction between two molecules, although the sigma profile is not obtained with a quantum mechanical calculation. Instead, three empirically calibrated electrostatic parameters for each functional group are used to describe the sigma profile of the entire molecule. Note that UNIFAC, COSMO, and F-SAC completely neglect the 3D structure of the molecule. Meaning that steric effects are not included and the pair formation at neighboring surface panels or neighboring functional groups is completely independent. For COSMO, this inadequacy can be partially tackled by using local sigma profiles.15,16 One can also obtain thermodynamic properties using forcefield-based molecular simulations. Simulations are more ab initio than any of the engineering models mentioned before. Several studies showed very accurate results;17−21 however, other studies showed that these results highly depend on the used force field.22−24 Moreover, with simulations requiring up to one billion Monte Carlo cycles for each composition, it is Received: June 25, 2015 Revised: September 26, 2015 Published: September 29, 2015 14200

DOI: 10.1021/acs.jpcb.5b06100 J. Phys. Chem. B 2015, 119, 14200−14209

Article

The Journal of Physical Chemistry B

distance D. This sampling procedure is repeated m times for molecules A and A, m times for molecules A and B, and m times for molecules B and B. The chosen value for m has to be large enough in order to obtain a diverse set of pair configurations, but not too large since the calculation cost scales linearly with m (in the calculations, we use m = 105). The distance of closest proximity is defined as the minimum intermolecular distance between two atoms in a pair configuration. Assume the position of atom k of the stationary molecule is given by a vector rk and the position of atom l of the displacing molecule is given by a vector rl, then the required distance of translation r along the x-axis to place atom l at a distance D from atom k can be calculated by solving the equation:25

also a time-consuming process. To reduce the computational cost, Blanco et al.25,26 introduced a pair sampling method for a quick estimation of the Flory−Huggins interaction parameter. In their technique, the intermolecular energy is calculated just between two molecules in direct contact with each other, in vacuum, as opposed to a full calculation in the condensed phase. However, because only pair configurations are sampled, it is difficult to take in a proper way the packing of molecules in the liquid phase into account. A possible solution was found by Akkermans27 who proposed to sample complete clusters of molecules instead of pairs. The method of Akkermans then uses an expression for the cluster free energy to calculate the Flory− Huggins interaction parameter of two-component mixtures. Here, we combine the power of the quasi-chemical derived activity coefficients directly with the use of a force field in the Pair Configurations to Molecular Activity Coefficients (PACMAC) method. The method is based on Blanco’s pair sampling technique, but it is extended in two ways. First, a better packing model is obtained by tracking the occupied surface in each pair configuration. Second, the derived activity coefficients are calculated from a quasi-chemical free energy model. The method is fast, because only pair configurations are sampled, and in principle with a relatively small parameter set, because it does not contain ad hoc functional groups. The paper is organized as follows. First, we discuss the sampling of pair energies, the calculation of the molecular surface, and the free energy model. Then we compare the calculated vapor−liquid equilibrium diagrams with experimental data from the NIST database and with the results of the COSMO-SAC model using sigma profiles from the VT-2005 Sigma Profile Database.28 This paper contains a proof of concept of the PAC-MAC model with limited validation. A comparison with a wide range of experimental data will be presented elsewhere.

|r × i + rl − rk|2 = D2

(1)

with i being the unit vector codirectional with the x-axis. The solution of eq 1 is given by r± = rk ·i − rl·i ±

(rl·i − rk ·i)2 − |rl − rk|2 + D2

(2)

The two solutions of eq 2 represent a translation in positive direction and a translation in negative direction, r+ and r− respectively. Equation 2 has to be calculated for every set of atoms k and l. The translation in positive direction required to make a configuration of two molecules, in which the minimum intermolecular distance between two atoms is equal to D, is then given by the maximum obtained value for r+; see Figure 1 for a 2D illustration.

■

THEORETICAL BASIS The PAC-MAC method predicts the free energy and related properties of a binary liquid mixture using the assumption that only interactions between neighboring molecules are relevant in a condensed phase.26,27 The possible configurations of two neighboring molecules in the mixture are given by a diverse set of sampled molecular pair configurations. Using classical force fields, the intermolecular energy between the two molecules is calculated and, by introducing the shielding surface of a molecule, PAC-MAC predicts which part on both molecules is occupied by the neighboring molecule. The calculation of the intermolecular energy and occupied shielding surface is performed for each sampled pair configuration. Then, the probability distribution of the pair configurations within a mixture in thermodynamic equilibrium is calculated by minimizing an expression for the free energy, given a constraint that each molecule has to be fully surrounded by other molecules. Sampling of Pair Configurations. The PAC-MAC model requires a varied set of all possible mutual orientations of two molecules within the first coordination shell of each other. First, we build and optimize the structure of the two different molecules (A and B) in the binary mixture using the CULGI software package.29 Then, a pair configuration is formed by giving these two molecules a random rotation, placing both geometric centers at the origin of a Cartesian coordinate system and translating one of the molecules along the x-axis until the distance of closest proximity is equal to a certain predefined

Figure 1. Illustration of the translation r in positive direction along the x-axis to make a pair configuration of a blue and red molecule with a distance of closest proximity being equal to D.

If the value of D is set equal to the sum of the van der Waals radii of atoms k and l, then the sampling procedure is exactly the same as the method of Blanco. Here, we use a different formulation for D: since we use a force-field-based method to calculate the intermolecular pair energy, we relate the value of D to the Lennard−Jones distance parameter, σ, multiplied with a dimensionless scaling factor s: D = s × (σkk + σll)

(3)

The minimum of a Lennard−Jones potential between a pair of atoms k and l, using Lorentz−Berthelot combining rules, is obtained at: 1 |rk − rl| = 21/6 × × (σkk + σll) (4) 2 Combining eq 3 and 4 implies that the distance between two charge-neutral atoms is energetically optimal with a scaling factor s equal to 21/6 × (1/2) = 0.56. However, in fluids, the distance of closest proximity between two neighboring molecules is not fixed. In order to sample pairs of neighboring 14201


Article


Figure 2. Distribution of the amount of pairs n(s), at a specific value of s, with a distance of closest proximity between two molecules being exactly equal to D, obtained with a MD simulation at 298 K and 1 atm. The value of s of a pair of molecules is calculated backward using eq 5. The limits of the domain defining the first coordination shell are for the tested molecules, excluding water, on average given by smin = 0.36 and smax = 0.72.

molecules at several distances D, it is necessary to define a range of possible s-values. Molecular dynamics simulations in a condensed phase simulation (NPT) were used to calculate the range wherein s should vary. We perform a backward calculation of s for which the corresponding value of D, calculated with eq 3, is exactly equal to the distance of closest proximity between two selected molecules. We can then express s by the following formula: s=

⎛ |r − rl| ⎞ min ⎜ k ⎟ k , l ∈ {atoms}⎝ σkk + σll ⎠

Table 1. Upper and Lower Limit of the First Coordination Shell and Coordination Number of the Molecules Given in Figure 2 Obtained with a MD Simulation at 298 K and 1 atm range of the first coordination shell

(5)

with the minimization performed over all combinations of atom k of one molecule and atom l of the other molecule. By calculating the scaling factor s for multiple pairs of molecules in the condensed phase, a distribution of the amount of pairs n(s) at a specific value of s is obtained. The distribution of a concise set of molecules is shown in Figure 2. Defining neighboring molecules as being the molecules within the first coordination shell implies that the limits of the range wherein s should vary, for the formation of a pair configuration, are given by the boundaries of the first peak in n(s). The lower limit of the first coordination shell is defined by the maximum value of s at which n(s) = 0 and the upper limit is defined by the first local minimum of n(s). The upper and lower limit of the first coordination shell of the 12 distributions presented in Figure 2 are given in Table 1 as well as the coordination number z defined by

∫lower‐limit

n(s) ds

(6)

Figure 2 and Table 1 show that, besides water, the entire selection of molecules have an approximately equal range of s defining the first coordination shell; the average values of the upper and lower limit, neglecting water, are respectively given by smax = 0.72 and smin = 0.36. The sampling of pair configurations is performed with a scaling factor s being a uniformly generated random number between smin and smax: si = U (smin , smax ) ∀ i ∈ {1, 2, .., m}

minimum s

maximum s

coordination number z

0.35 0.35 0.36 0.35 0.37 0.37 0.36 0.36 0.36 0.37 0.36 0.38

0.72 0.68 0.73 0.74 0.74 0.73 0.70 0.79 0.73 0.65 0.70 0.51

10.97 9.62 12.00 11.87 12.19 11.09 10.72 12.26 12.67 8.95 10.93 4.43

between neighboring molecules. Both phenomena are explained by the ability of a water molecule in forming up to four hydrogen bonds with surrounding water molecules. Intermolecular Energy and Occupied Shielding Surface. For each sampled pair configuration i, i ∈ {1,2,..,m}, two calculations are performed by PAC-MAC: calculation of the pairwise intermolecular energy and calculation of the part of both molecules occupied by the neighboring molecule. Tracking of the occupied parts of the molecules is necessary because, to describe a condensed system accurately, we force each molecule to be completely surrounded by a single shell of nearest neighbor molecules. The pairwise intermolecular interaction energy, wi, is calculated using classical force fields. The occupied part of a pair of molecules is defined by the area of the molecular shielding surface covered by the other molecule. The shielding surface of a molecule is given by the external surface of spheres around the atomic nuclei. The shielding radius is given by R shield, k = ssegσkk (8)

upper‐limit

z=

molecule acetone acetonitrile 1-butanol n-butylamine 1,2-dichloroethane diethyl ether ethanol n-hexane 1-hexanol methanol 1-propanol water

(7)

The low coordination number of water indicates a liquid phase preserving much of the ice-like tetrahedral structure. The small value for smax in water suggests a stronger attractive force

The molecular shielding surface is divided in L surface segments of approximately equal size. In principle, the chosen 14202


Article

The Journal of Physical Chemistry B value for L could be different for both molecules. In general, the value is a trade-off between a high resolution of the molecule (high L) and a fast PAC-MAC calculation (low L). For each sampled configuration of a pair of molecules it is recorded which of the L surface segments are occupied and which not. L needs to be sufficiently large in order to have in each pair configuration a decent amount of surface segments, around 20, occupied at both molecules. In the calculations, we use L = 200 independent of molecular type. Given a coordination number of around 10 (see Table 1), implies a resolution per molecular contact of about 20 surface panels, which is of sufficient accuracy. See Figure 3 for an illustration of the shielding surface of acetone:

Figure 5. Three sampled pair configurations of acetone and dimethyl ether. The occupied surface segments are black colored, and the corresponding intermolecular energies are indicated at the top.

number within PAC-MAC is given in the subsection Distribution of Pair Configurations. Because the sampled pairs of neighboring molecules are defined to be the molecules within the first coordination shell implies that the coordination number obtained by PAC-MAC should be equal to the coordination number obtained in a MD simulation. It is shown in Figure 6 that the optimal values of sseg, for which the coordination number obtained by PAC-MAC is exactly equal to the coordination number given in Table 1, are close to each other for all 12 selected molecules except water. The arithmetic mean, neglecting water, is given by sseg = 0.66 which is therefore the value we use within the PAC-MAC model. Note that eqs 3 and 4 show an energetically optimal distance of closest proximity D between two charge-neutral atoms equal to the summation of both atomic Lennard−Jones σ-parameters multiplied with a scaling factor s = 0.56. The shielding surface of a molecule is defined using a comparable scaling function for the atomic shielding radius: Rshield,k is given by the Lennard− Jones σ-parameter of atom k multiplied with a scaling factor sseg (eq 8). This also indicates a comparable scaling factor, sseg ≈ 0.56, in the representation of the shape of a molecule. However, in reality, the molecular surface in a condensed system is not fully covered; holes are present because steric effects between neighboring molecules avoid a complete occupation of the surface. In order to correct for this effect, the molecules need to be made artificially bigger in a PAC-MAC simulation to obtain a proper coordination number. This proves the observed average value for the atomic radius scaling factor: sseg = 0.66. Free Energy in PAC-MAC. We take a liquid binary solution with a total number of N molecules consisting of NA molecules A and NB molecules B:

Figure 3. Acetone molecule. Left: Ball-and-stick representation. Right: Shielding surface, using sseg = 0.66, divided in L = 200 surface segments.

A surface segment is occupied if there is no space for another molecule to be at that position. This is the case for the surface segments of which a line through their center, coming from the geometric center of the molecule, will pierce at least one surface segment of the neighboring molecule. See Figure 4 for

N = NA + NB

(9)

The N molecules in the mixture will form pairs between each other. The set of sampled pair configurations represent all the possible pairs of neighboring molecules that can be obtained. The total amount of pair configuration i observed in the mixture between molecules I and molecules J is denoted by Ni,IJ. The expression for the mixing free energy is given by a summation of the ideal free energy, excess combinatorial free energy and interaction free energy:

Figure 4. 2D illustration of the shielding surface area occupied (red) in a pair configuration of two simplified molecules. The left molecule consists of two atoms, and the right molecule of a single atom. The dashed lines start in the geometric center of the molecules and touch the surface of the neighboring molecule. The intercept is the boundary of the occupied surface.

F mix = F id + F comb + F int

clarification using a 2D illustration. A visualization of the occupied surface segments for various randomly sampled pair configurations of acetone and dimethyl ether is shown in Figure 5. By increasing the value of sseg in eq 8, the occupied surface of a molecule also increases and therefore the total amount of pairs a molecule can form (the coordination number) decreases. Figure 6 shows the relationship between sseg and the coordination number for the set of molecules of Figure 2. More information about the calculation of the coordination

(10)

id

The ideal free energy F is given by the ideal entropy of mixing: F id = NA ln(xA) + NB ln(xB) kBT

(11)

with xA and xB being the mole fractions of, respectively, molecule A and molecule B in the binary mixture. A random probability distribution of pairs between the molecules implies a total entropy of the solution equal to Fid. The quasi-chemical approximation is used to correct the entropy for mixtures with 14203


Article


Figure 6. Lines: Coordination number obtained by PAC-MAC as a function of the parameter sseg. Diamonds: Point whereby the coordination number obtained by PAC-MAC is equal to the coordination number obtained in a MD simulation at 298 K and 1 atm. Vertical dashed line: Value for parameter sseg in PAC-MAC, sseg = 0.66; arithmetic mean of sseg of the diamonds neglecting water.

their value is obtained by minimization of the free energy function (eq 10). The calculated values of xi,IJ and Ni,IJ therefore represent the pair fraction and total pair quantity of configuration i observed in a mixture at thermodynamic equilibrium, provided a correct expression for the free energy. Packing in PAC-MAC. The PAC-MAC model enforces every surface segment to be fully occupied by neighboring molecules:

a nonrandom distribution of pair configurations. The excess combinatorial free energy is given by30 F comb = kBT

m

∑ Ni ,AA i=1

⎛x m⎞ i , AA ⎟+ ln⎜⎜ 2 ⎟ ⎝ yA ⎠ ⎛x

m

+

∑ Ni ,BB ln⎜⎜ ⎝

i=1

⎛x

⎞ ⎟⎟ ⎝ 2yA yB ⎠

m

∑ Ni ,AB ln⎜⎜ i=1

i , ABm

⎞

i , BBm ⎟ yB2 ⎟⎠

(12)

m

NA =

in which m is the total number of sampled pair configurations between two molecules, yA and yB are the “coordinationequivalent” fractions of, respectively, molecule A and molecule B: m

yA =

∑ xi ,AA + i=1

1 2

i=1

m

∑ xi ,AB

NB =

(13)

i=1

m

(15)

i=1

wi , AA kBT

m

+

∑ Ni ,AB i=1

wi , AB kBT

m

+

∑ Ni ,BB i=1

i=1

The term oj,I,i,IJ describes whether surface segment j of molecule I is occupied (oj,I,i,IJ = 1) or not (oj,I,i,IJ = 0) in pair configuration i between molecule I and J. The packing constraints are an adaption to the method of Blanco et al. to force a complete surrounding of the molecules. Distribution of Pair Configurations. To fulfill the packing constraints we add eqs 17 and 18 to eq 10 together with a set of Lagrange multipliers ξj,A and ξj,B (j ∈ {1,2,..,L}):

with z being the average coordination number of the N molecules. The equation is comparable to the formulation used in the modified quasi-chemical model of Pelton and Blander,30,31 with the difference that the model of Pelton and Blander contains only a single possible pair configuration between a molecule I and a molecule J whereas PAC-MAC contains m possible pairs. The interaction free energy Fint is given by

∑ Ni ,AA

m

(18)

(14)

1 zNxi , IJ 2

(17)

∑ Ni ,BBoj ,B,i ,BB + ∑ Ni ,ABoj ,B,i ,AB ∀ j ∈ {1, 2, .., L} i=1

Ni,IJ is the amount of pair configurations i, i ∈ {1,2,..,m}, between molecules I and J in the mixture, and xi,IJ is the corresponding pair fraction. The relation between Ni,IJ and xi,IJ is given by

F int = kBT

i=1

j ∈ {1, 2, .., L}

m

yB = 1 − yA

Ni , IJ =

m

∑ Ni ,AAoj ,A ,i ,AA + ∑ Ni ,ABoj ,A ,i ,AB ∀

3=

F mix − kBT

m

⎞

m

⎝ i=1

j=1

⎛

L

−

⎛

L

∑ ξj , A⎜⎜∑ Ni , AAoj , A , i , AA + ∑ Ni , ABoj , A , i , AB − NA ⎟⎟ m

⎠

i=1

⎞

m

∑ ξj , B⎜⎜∑ Ni , BBoj , B, i , BB + ∑ Ni , ABoj , B, i , AB − NB⎟⎟ j=1

⎝ i=1

⎠

i=1

(19)

The constrained free energy is minimized with respect to Ni,IJ and the surface segment packing correction terms ξj,I. This results in a set of equations:

wi , BB kBT (16)

xi , AA =

The interaction free energy is obtained by summing the intermolecular energy of all pair configurations in the mixture. Long-range interactions between non-neighboring molecules are not taken into account by the PAC-MAC model. Note that the set of pair quantities Ni,IJ within eqs 12 and 16 are variables. Ni,IJ is related to the pair fraction xi,IJ via eq 15 and

xi , AB =

⎞ ⎛ w i , AA exp⎜⎜ − + ∑j ξj , Aoj , A , i , AA ⎟⎟ m ⎠ ⎝ kBT

yA 2

2yA yB m

(20)

⎛ wi , AB ⎞ exp⎜ − + ∑j ξj , Aoj , A , i , AB + ∑j ξj , Boj , B, i , AB⎟ ⎝ kBT ⎠ (21)

14204


Article


Figure 7. Distribution of the amount of pairs n(s), at a specific value of s, with a distance of closest proximity between two molecules being exactly equal to D of eq 3, obtained with MD (blue curve) and with PAC-MAC (red curve) at 298 K for diethyl ether (left), acetonitrile (middle), and water (right) in a pure condensed phase. The limits of the domain of the PAC-MAC distribution function are given by smin = 0.36 and smax = 0.72.

⎞ ⎛ wi , BB exp⎜⎜ − + ∑j ξj , Boj , B, i , BB⎟⎟ m kBT ⎠ ⎝

m 1 ⎛ xA = z⎜⎜∑ xi , AAoj , A , i , AA + 2 ⎝ i=1

xB =

1=

L

yB2

xi , BB =

ln(γIres) =

⎞ ∑ xi ,ABoj ,A ,i ,AB⎟⎟ ⎠ i=1

∑ xi ,BBoj ,B,i ,BB⎟⎟

m

⎞

i=1

⎠

m

m

∑ xi ,AA + ∑ xi ,AB + ∑ xi ,BB i=1

i=1

i=1

is given by the difference of the summation of the packing correction terms belonging to the surface segments of molecule I between the mixed phase and the pure phase. This suggests that in PAC-MAC the Lagrange multipliers ξj,I represent the segmental free energies of molecule I.

(23)

■

RESULTS AND DISCUSSION The Lennard−Jones σ-parameters, required for the formation of pair configurations and for defining the shielding surface, are taken from the OPLS-aa force field.33 The same force field is also used for the calculation of the intermolecular energy of each pair configuration. The calculation time for the sampling of 3m pair configurations (m = 105) is about 5 min on a single core (Intel Xeon E5-2620). Solving the systems of equations takes about 20 s for each chosen mole fraction. First Coordination Shell. As an internal consistency check, we calculated the molecular pair distribution function of the first coordination shell from PAC-MAC. The molecular pair distribution is obtained by multiplying the calculated pair fractions, xi,II, with the coordination number, z, from a PACMAC calculation of a pure condensed molecule I and plotting the obtained values as a function of the scaling factor s, calculated using eq 5. The distribution functions of PAC-MAC closely fit the distribution functions of the first coordination shell of MD, given in Figure 2, for the molecules for which the chosen values of parameters smin (smin = 0.36), smax (smax = 0.72), and sseg (sseg = 0.66) are in good agreement with the optimal values (Table 1 and Figure 6) for these molecules; see Figure 7. It is also shown that PAC-MAC yields a less sharp density peak of the first coordination shell than the peak obtained by MD. This is because in PAC-MAC only nearest neighbor pair interactions are involved, so there is no pressure of second or higher coordination shells for compressing the first one. Comparison with Experimental Data and COSMOSAC. The bubble point equation for a binary mixture under the extended Raoult’s law, assuming gaseous ideality, is given by

(24)

(25)

The solution for xi,AA, xi,AB, xi,BB, ξj,A, ξj,B, and z is found by first substituting eqs 20−22 into eqs 23−25 and then solving iteratively using the Newton−Raphson method. Calculating Activity Coefficients with PAC-MAC. The solutions of eqs 20−25 can subsequently be used to calculate the activity coefficients of the molecules in the mixture. The activity coefficient of molecule I is split into two contributions:7 γI = γIcombγIres

(26)

The combinatorial contribution to the activity coefficient is given by the Staverman−Guggenheim equation:32 ln(γIcomb,SG) = 1 − −

⎛ϕ ⎞ + ln⎜ I ⎟ xI ⎝ xI ⎠

ϕI

⎛ ϕ ⎞⎞ ϕ zI ⎛ ⎜⎜1 − I + ln⎜ I ⎟⎟⎟ 2⎝ θI ⎝ θI ⎠⎠

(27)

θI and φI are, respectively, the surface fraction and the volume fraction of molecule I in the mixture; calculated from the area of the shielding surface and the volume inside. The variable zI is the coordination number of molecule I, defined by zI = z

yI xI

P(xAliq) = xAliqγAPA* + (1 − xAliq)γBPB*

(28)

1 ⎛ ∂F mix ∂F pure ∂F id ⎞ − − ⎟ ⎜ kBT ⎝ ∂NI ∂NI ∂NI ⎠

(31)

xliq A

with denoting the mole fraction of molecule A in the liquid phase and PI* being the vapor pressure of pure component I. The corresponding mole fraction of molecule A in the gaseous phase at the dew point curve can be calculated using

The value of the residual term, γres I , is given by ln(γIres) =

(30)

j=1

ln(γres I )

m

m 1 ⎛ z⎜⎜∑ xi , ABoj , B , i , AB + 2 ⎝ i=1 m

j=1

(22)

L

pure ∑ ξjmix , I − ∑ ξj , I

(29)

xAgas =

Using eq 10 and taking the derivatives of eq 29 results in a remarkably simple expression for γres I : 14205

xAliqγAPA* P(xAliq)

(32) DOI: 10.1021/acs.jpcb.5b06100 J. Phys. Chem. B 2015, 119, 14200−14209

Article


Figure 8. Concise set of vapor−liquid equilibrium (VLE) diagrams obtained with COSMO-SAC and PAC-MAC compared with experimental data for 1-butanol−n-butylamine (313, 14 K), 1-propanol−n-hexane (323, 137 K), ethanol−acetonitrile (333, 15 K), methanol−water (323, 15 K), diethyl ether−acetone (302, 984 K), and 1-butanol−1,2-dichloroethane (323,137 K).

The experimental data is provided by the NIST database;1 COSMO-SAC results are obtained using sigma profiles from the freely available VT-2005 Sigma Profile Database.28 As a proof of concept, we present in Figure 8 the results of PACMAC for a concise, but varied, set of six binary solvent mixtures. Figure 8 illustrates a comparable accuracy for the prediction of VLE diagrams using the PAC-MAC method and the COSMO-SAC method. The upper three plots are examples in which PAC-MAC results in a better fit of the experimental data, the lower three plots are examples in which COSMO-SAC excels. An explanation of the observed results is given for two of these six VLE’s. PAC-MAC was expected to outperform COSMO-SAC for the case of 1-butanol−n-butylamine because it is known that the unpaired electron in the cavity causes poor results for amines in COSMO-SAC.14 On the other hand, the PAC-MAC method performs rather poorly for the case of methanol−water. Also other VLEs of binary mixtures involving water are poorly described by PAC-MAC. The flaw can be explained by the sparameters: the preset values (smin = 0.36, smax = 0.72, sseg = 0.66) are correct for most of the molecules; however, for water, they are invalid. The right plot of Figure 7 predicts a coordination number of 8.8 for water using PAC-MAC, whereas the MD simulation of Figure 2 yields a coordination number of 4.4. Consequently, the obtained activity coefficients are overestimated. Molecular Specific s-Parameters. The inconsistency of the parameters for water is overcome by using molecular specific s-parameters. These parameters are obtained using the distribution function n(s) calculated by a MD simulation of the pure condensed phase. The limits of the first coordination shell (see Table 1) define the values of smin and smax; the value of sseg is optimized in order to obtain a coordination number, calculated by the PAC-MAC model, being exactly equal to the coordination number obtained by MD. The molecular specific s-parameters are given in Table 2.

Table 2. Molecular Specific Parameters smin, smax, and sseg Calculated Using a MD Simulationa molecule

smin

smax

sseg

acetone acetonitrile 1-butanol n-butylamine 1,2-dichloroethane diethyl ether ethanol n-hexane 1-hexanol methanol 1-propanol water

0.35 0.35 0.36 0.35 0.37 0.37 0.36 0.36 0.36 0.37 0.36 0.38

0.72 0.68 0.73 0.74 0.74 0.73 0.70 0.79 0.73 0.65 0.70 0.51

0.67 0.70 0.64 0.65 0.64 0.68 0.67 0.65 0.61 0.71 0.67 0.77

a

smin and smax are, respectively, given by the lower limit and upper limit of the first coordination shell. The value for sseg is optimized in order to have a coordination number, calculated by PAC-MAC, being equal to the coordination number calculated by MD.

As expected the molecular optimized values of Table 2 are close to the general chosen s-parameters (smin = 0.36, smax = 0.72, sseg = 0.66) except in the case of water. The six VLE diagrams obtained with PAC-MAC of Figure 8 are recalculated using the values of Table 2; see Figure 9. The VLE of methanol−water is strongly improved by the molecular specific s-parameters, because of the large deviation from the general parameters used in Figure 8. As expected, the influence of molecular specific s-parameters on the other VLEs is rather small. In general, the accuracy of the PAC-MAC method increases by using a set of molecularly defined sparameters. However, the VLE of ethanol−acetonitrile shows that improvement is not always achieved. Besides the sparameters, there are other factors having influence on the results; one of these factors is the accuracy of the force field. To 14206


Article


Figure 9. Vapor−liquid equilibrium (VLE) diagrams obtained with COSMO-SAC and PAC-MAC compared with experimental data for six binary solvent mixtures. The used s-parameters are taken from Table 2.

Figure 10. VLE diagrams obtained with PAC-MAC and compared with experimental data and Monte Carlo calculations. The general chosen sparameters are used: smin = 0.36, smax = 0.72, sseg = 0.66. Left: argon−methane (122.89 K), GEMC simulation by Nguyen et al.18 Right: methanol− ́ et al.19 MTBE (338.15 K), RGEMC simulation by Lisal

avoid influence of the force field, we test the PAC-MAC model with molecular simulations using the same force field. Comparison with Molecular Simulations. In Figure 10, we compare the results of PAC-MAC with the results of Gibbs ensemble Monte Carlo simulations for a mixture of argon− methane18 and a mixture of methanol−MTBE19 using the same force field. The PAC-MAC method only requires the sampling of 3m pair configurations (m = 105) between two molecules to obtain a full VLE diagram, whereas the Gibbs ensemble Monte Carlo simulations were performed by sampling 1.5−3 × 105 cycles consisting of 5514 steps in a box containing 512 particles for just a single liquid fraction at the VLE. Assuming the calculation time of a MC step being equal to the sampling of a single pair configuration implies that PAC-MAC is at least 4 orders of magnitude faster than a corresponding MC simulation for obtaining VLE data at ten different liquid fractions. The comparison shows a very good agreement with the noble gas mixture, but with the compound molecules much less so. This is because the PAC-MAC model contains a few assumptions, in addition to condensed phase molecular simulations, affecting the results: (i) Only direct pair interactions are assumed to be relevant. Indirect effects, due to multiple body interactions, are

only taken into account by the excluded surface of each interaction. (ii) Only nearest neighbor pair configurations are sampled. Long-range interactions between non-neighboring molecules are neglected. (iii) A molecule is completely surrounded by other molecules. A possibly inaccurate assumption, especially if a molecule contains strong specific interactions, such as water. Allowing holes around the molecular shielding surface might result in a better representation of the packing around a molecule. (iv) The conformation of the molecules is fixed. It is more realistic to extend the pair sampling with multiple molecular conformations. (v) The expressions for the excess combinatorial free energy (eq 12) and the Staverman−Guggenheim correction (eq 27) are not exact.30,34

■

CONCLUSION In this paper, a force-field-based quasi-chemical method for the prediction of vapor liquid equilibrium diagrams is presented: PAC-MAC. The method is based on Blanco’s pair sampling technique with two adjustments. First, the occupied surface of every sampled pair configuration is measured in order to obtain 14207


Article


(14) Oldland, R. J. Predicting Phase Equilibria Using COSMO-Based Thermodynamic Models and the VT-2004 Sigma-Profile Database. Masters Thesis, Virginia Polytechnic Institute and State University, Blacksburg, VA, 2004. (15) Thormann, M.; Klamt, A.; Wichmann, K. COSMOsim3D: 3DSimilarity and Alignment Based on COSMO Polarization Charge Densities. J. Chem. Inf. Model. 2012, 52, 2149−2156. (16) Klamt, A.; Thormann, M.; Wichmann, K.; Tosco, P. COSMOsar3D: Molecular Field Analysis Based on Local COSMO σ-Profiles. J. Chem. Inf. Model. 2012, 52, 2157−2164. (17) Caleman, C.; van Maaren, P. J.; Hong, M.; Hub, J. S.; Costa, L. T.; van der Spoel, D. Force Field Benchmark of Organic Liquids: Density, Enthalpy of Vaporization, Heat Capacities, Surface Tension, Isothermal Compressibility, Volumetric Expansion Coefficient, and Dielectric Constant. J. Chem. Theory Comput. 2012, 8, 61−74. (18) Nguyen, V. T.; Fan, C.; Razak, M. A.; Do, D. D.; Nicholson, D.; Ustinov, E. Development of Kinetic Monte Carlo and Bin-Monte Carlo Schemes for Simulation of Mixtures - Vapor-Liquid Equilibria & Adsorption. Chem. Eng. Sci. 2013, 102, 220−226. (19) Lísal, M.; Smith, W. R.; Nezbeda, I. Accurate Computer Simulation of Phase Equilibrium for Complex Fluid Mixtures. Application to Binaries Involving Isobutene, Methanol, Methyl tertButyl Ether, and n-Butane. J. Phys. Chem. B 1999, 103, 10496−10505. (20) Agrawal, R.; Wallis, E. P. Molecular Simulation of the VaporLiquid Equilibrium of Methanethiol plus Propane Mixtures. Fluid Phase Equilib. 1997, 131, 51−65. (21) Johansson, E.; Bolton, K.; Theodorou, D. N.; Ahlström, P. Monte Carlo Simulations of Equilibrium Solubilities and Structure of Water in n-Alkanes and Polyethylene. J. Chem. Phys. 2007, 126, 224902. (22) 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. (23) 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. (24) 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. (25) Blanco, M. Molecular Silverware. I. General Solutions to Excluded Volume Constrained Problems. J. Comput. Chem. 1991, 12, 237−247. (26) Fan, C. F.; Olafson, B. D.; Blanco, M.; Hsu, S. L. Application of Molecular Simulation to Derive Phase Diagrams of Binary Mixtures. Macromolecules 1992, 25, 3667−3676. (27) Akkermans, R. L. C. Mesoscale Model Parameters from Molecular Cluster Calculations. J. Chem. Phys. 2008, 128, 244904. (28) Mullins, E.; Oldland, R.; Liu, Y. A.; Wang, S.; Sandler, S. I.; Chen, C.-C.; Zwolak, M.; Seavey, K. C. Sigma-Profile Database for Using COSMO-Based Thermodynamic Methods. Ind. Eng. Chem. Res. 2006, 45, 4389−4415. (29) Fraaije, J. G. E. M.; Nath, S. K.; Male, J. v.; Becherer, P.; Wolterink, J. K.; Handgraaf, J.-W.; Case, F.; Tanase, C.; Gracià, R. S. Culgi, Version 9.0.1; Culgi B.V.: Leiden, the Netherlands, 2014. (30) Pelton, A. D.; Degterov, S. A.; Eriksson, G.; Robelin, C.; Dessureault, Y. The Modified Quasichemical Model I - Binary Solutions. Metall. Mater. Trans. B 2000, 31, 651−659. (31) Blander, M.; Pelton, A. D. Thermodynamic Analysis of Binary Liquid Silicates and Prediction of Ternary Solution Properties by Modified Quasichemical Equations. Geochim. Cosmochim. Acta 1987, 51, 85−95. (32) Staverman, A. J. The Entropy of High Polymer Solutions. Generalization of Formulae. Recl. Trav. Chim. Pays-Bas 1950, 69, 163− 174. (33) 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.

a realistic packing of the molecules. Second, the method is quasi-chemical correct. The obtained VLE diagrams by PAC-MAC are compared with experimental data and with VLEs obtained by COSMOSAC. It is shown that, even with the relatively simple classical OPLS-aa force field and only three additional distance scaling parameters, the PAC-MAC method can produce VLE diagrams fast, within 15 min on a single computational core, and with an accuracy comparable to COSMO-SAC. Binary mixtures involving water are an exception for which the PAC-MAC method is highly inaccurate. By using molecular specific s-parameters, obtained by short MD simulations on pure liquids, the PAC-MAC method also obtains VLE diagrams with high accuracy for water containing systems.

■ ■

AUTHOR INFORMATION

Notes

The authors declare no competing financial interest.

ACKNOWLEDGMENTS We are very grateful to SABIC and, in particular, K. Remerie for providing financial support for our research, and we would like to thank J. v. Male for useful discussions.

■

REFERENCES

(1) 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. (2) 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. (3) 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. (4) Economou, I. G. Statistical Associating Fluid Theory: A Successful Model for the Calculation of Thermodynamic and Phase Equilibrium Properties of Complex Fluid Mixtures. Ind. Eng. Chem. Res. 2002, 41, 953−962. (5) Fowler, R. H.; Guggenheim, E. A. Statistical Thermodynamics of Super-Lattices. Proc. R. Soc. London, Ser. A 1940, 174, 189−206. (6) Prausnitz, J. M.; Lichtenthaler, R. N.; de Azevedo, E. G. Molecular Thermodynamics of Fluid-Phase Equilibria, 3rd ed.; Prentice-Hall, Inc.: Upper Saddle River, NJ, 1999. (7) Fredenslund, A.; Jones, R. L.; Prausnitz, J. M. GroupContribution Estimation of Activity Coefficients in Nonideal Liquid Mixtures. AIChE J. 1975, 21, 1086−1099. (8) 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. (9) Klamt, A. COSMO-RS: From Quantum Chemistry to Fluid Phase Thermodynamics and Drug Design; Elsevier: Amsterdam, 2005. (10) Lin, S.-T.; Sandler, S. I. A Priori Phase Equilibrium Prediction from a Segment Contribution Solvation Model. Ind. Eng. Chem. Res. 2002, 41, 899−913. (11) Soares, R. d. P.; Gerber, R. P. Functional-Segment Activity Coefficient Model. 1. Model Formulation. Ind. Eng. Chem. Res. 2013, 52, 11159−11171. (12) Wu, H. S.; Sandler, S. I. Use of Ab Initio Quantum Mechanics Calculations in Group Contribution Methods. 1. Theory and the Basis for Group Identifications. Ind. Eng. Chem. Res. 1991, 30, 881−889. (13) Wu, H. S.; Sandler, S. I. Use of Ab Inito Quantum Mechanics Calculations in Group Contribution Methods. 2. Test of New Groups in UNIFAC. Ind. Eng. Chem. Res. 1991, 30, 889−897. 14208


Article

The Journal of Physical Chemistry B (34) 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.

14209


Force-Field Based Quasi-Chemical Method for Rapid Evaluation of

Recommend Documents