- Email: [email protected]
Activity coefficient prediction by osmotic molecular dynamics
Fluid Phase Equilibria 193 (2002) 53–73
Activity coefficient prediction by osmotic molecular dynamics Paul S. Crozier∗ , Richard L. Rowley Department of Chemical Engineering, Brigham Young University, Provo, UT 84602, USA Received 6 March 2001; accepted 24 September 2001
Abstract The osmotic molecular dynamics method (OMD) is used to calculate activity coefficients for vapour–liquid equilibria (VLE) and liquid–liquid equilibria (LLE) predictions. The previously reported OMD methodology is refined and applied to mixtures of polar, structured molecular fluids. Other computer simulation approaches to phase equilibria prediction are discussed briefly, and comparison to recent Gibbs Ensemble Monte Carlo (GEMC) results is made. OMD-predicted activity coefficients are compared to experimentally-measured activity coefficients for six industrially-significant binary mixtures (methanol/n-hexane, n-hexane/n-pentane, chloroform/acetone, chloroform/methanol, methanol/water, chloroform/n-hexane). Molecular model inadequacies, especially cross-parameters between unlike molecules, are shown. A single cross-parameter for the acetone/chloroform binary is modified to produce good agreement with experimentally-measured activity coefficients. Also, OMD-derived LLE predictions are produced for the methanol/n-hexane system and compared with experimentally-measured LLE data. © 2002 Elsevier Science B.V. All rights reserved. Keywords: Computer simulation; Methods of calculation; Molecular interactions; Vapour–liquid equilibria; Liquid–liquid equilibria; Chemical potential
1. Introduction For decades, scientists and engineers have sought to predict activity coefficients and excess thermodynamic properties based upon molecular-level chemical characteristics. As computer speed and availability continue to increase, predicting thermodynamic properties by simulating substances on the molecular-level has gained acceptance. Thermodynamic integration [1–6], particle insertion [7,8], particle growing [9–22], distribution-histogram methods [23–27], and osmotic methods [28–32] are some of the more common methods currently used for the calculation of chemical potential by computer simulation. Kofke and Cummings [33,34] offer excellent reviews of the various methods and explain the inherent ∗
Corresponding author. Present address: Sandia National Laboratiories, P.O. Box 5800, MS 0316 Albuquerque, NM 87185-0316, USA. Tel.: +1-505-845-9714; fax: +1-505-284-5451. E-mail address: [email protected] (P.S. Crozier). 0378-3812/02/$ – see front matter © 2002 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 8 - 3 8 1 2 ( 0 1 ) 0 0 7 3 4 - 8
54
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
error in particle deletion schemes [35–40]. Insertion methods are widely used, the Gibbs Ensemble Monte Carlo (GEMC) method being the most widely accepted for VLE prediction. In this study, we refine the osmotic molecular dynamics method (OMD) [29–31] that can be used to calculate chemical potential and phase equilibrium. Unlike GEMC methods, OMD does not require costly particle insertions, and is not limited to the calculation of chemical potentials at phase boundaries. Using the same molecular models, we compare OMD results with recent GEMC phase equilibrium predictions. In addition, we compare OMD activity coefficient predictions to experimentally-measured activity coefficients for VLE and LLE systems. 2. OMD methodology 2.1. OMD thermodynamics As discussed in [29–31], activity coefficients can be determined by osmosis experiments in which a pure liquid solvent is separated from a liquid mixture of solvent and solute by a semi-permeable membrane, permeable only to solvent molecules. The activity coefficient is calculated according to the thermodynamic identity: 1 P0 +Posm ˜ ln γi = −ln xi − Vi dP (1) RT P0 where γ i is the activity coefficient of the permeable component (solvent) at a mole fraction xi in the mixture, P0 the pressure on the pure solvent side of the membrane, and P0 + P osm is the pressure on the liquid mixture side of the membrane. The molar volume of solvent, V˜i , is left inside the integral and considered to be a function of pressure. This molar volume versus pressure relationship can be determined by simulation as described in [30]. 2.2. Membrane interactions We used an artificial fixed-in-space membrane that is permeable only to solvent molecules and of the same functionality as that given by Rowley et al. [29]. Rowley et al. used the OMD method to calculate activity coefficients for Lennard–Jones (LJ) fluid mixtures, and modelled the interaction between the membrane and solute particles with one-dimensional, LJ potentials truncated at the potential-well minimum and shifted up in energy by ε, the LJ energy interaction parameter. This provided a continuous, purely repulsive potential based on the distance of the particle from the membrane. Their method was modified in this study to accommodate heterogeneous interaction sites on larger structured molecules. Each site i on solute (impermeable) molecules was set to interact with the membrane according to the following force formula: σi,i 6 σi,i 6 Fi (z) = 24εi,i (2) − 1 , for z < 21/6 σi,i 2 z z and Fi (z) = 0,
for z ≥ 21/6 σi,i
(3)
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
55
where z is given by z = 21/6 σi,i − d
(4)
and d is the perpendicular distance between the molecule’s centre of mass and the membrane surface. The value of d is defined to be positive only after the centre of mass of the solvent molecule has passed through the membrane, and is negative otherwise (no interaction). Although the value of the activity coefficient obtained from OMD simulations is independent of the form of the purely repulsive potential chosen for the solute-membrane interactions (the membrane must simply be semi-permeable), using a soft repulsive force like Eq. (2) avoids impulses that can occur with harder potentials during the finite time step associated with the numerical integration of the equations of motion. The form of Eq. (4) is chosen so that a molecule uniformly feels the repulsion when it crosses the membrane surface. The repulsion is exerted on all sites simultaneously, regardless of the molecule’s orientation. 2.3. Slot definition The simulation cell must be of adequate length in the z-direction (perpendicular to the membrane surface) to allow the formation of a fluid region that is isolated from, and independent of, the artificial membrane. Even though the chosen membrane described in the earlier sections is relatively undisruptive of the natural liquid structure, a volume of fluid within a fixed distance from the membrane is excluded from property calculations to ensure that only fluid whose structure is not perturbed by the membrane surface is used to calculate the activity coefficient. But it is desirable to choose small values for the z-direction box dimensions in order to minimise the time required for diffusion and equilibration. For the purposes of this work, the central “slot” region of unperturbed fluid was defined as all molecules with centre of mass at least 10 Å away from the membrane surface. Simulation studies [41,42] of net-neutral systems that include long-range Coulumbic interactions, similar to the present study, suggest that a 10 Å interfacial layer adequately shields even long-range interactions. Benchmark tests, discussed in the following sections, confirm that this distance is sufficient for the molecules of this study. Smaller distances were found to introduce error into the calculation. 2.4. Thermostatting Precise temperature control of the system is crucial to obtaining accurate activity coefficient predictions. Thermostatting of the entire simulation cell to the set-point temperature does not guarantee that the individual slot regions will be held at precisely the desired temperature during mass-transfer [43], especially since the membrane artificially adds energy to the system in the excluded volume regions. We independently thermostat the slot and adjacent exclusion volumes of each half-cell to ensure uniform temperature at the set-point throughout. The population of molecules pertaining to any given slot or exclusion volume is updated every 100 time steps, at which point the net momentum of that population of molecules in each direction is set to zero by velocity re-scaling, followed by an additional velocity re-scaling that fixes the population at the set-point temperature. A Gaussian thermostat [44,45] is used to hold each population’s temperature constant between updates. The zeroing of net momenta ensures that a net drift of the population that can accumulate in OMD simulations does not produce an erroneously high temperature measurement [46].
56
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
Temperature is calculated from the translational motion using 2 i mi vi T = k(3N − NRDF)
(5)
where T is the temperature, mi and v i are the mass and velocity of molecule i, k the Boltzmann constant, N the number of molecules in the population for a molecular thermostat or the number of atoms for an atomistic thermostat, and NRDF is the number of removed degrees of freedom (d.f.). One d.f. is removed due to the thermostat. In a regular MD simulation, momentum is conserved in all three directions for the thermostatted population; so three additional d.f. must be removed. However, in an OMD simulation, the thermostatted population is continually changing; new populations have not been previously thermostatted to the set-point temperature. Also, momentum is not conserved for the thermostatted populations since the molecules of each slot and exclusion volume interact with the fixed, rigid membrane and/or molecules outside that region. No d.f. were removed due to the bond and angle constraints because molecular temperatures, rather than atomic temperatures, were computed. Hence, for the simulations of this work, NRDF = 1 in Eq. (5). 2.5. Osmotic pressure calculation The system temperature has a strong effect on the system pressure and on the osmotic pressure difference. Prediction of activity coefficients by Eq. (1) is extremely sensitive to the precise determination of the equilibrium osmotic pressure difference, Posm . Although this pressure difference can be found by taking the difference between the pressure found on the mixture side of the membrane (Side A) and the pressure found on the solvent side of the membrane (Side B), each being calculated using the virial equation, this procedure is computationally expensive and yields large instantaneous fluctuations in the calculation of Posm . We have found that this method of finding Posm is especially prone to error if discontinuous molecular models (e.g. cutoff fluids) are used. A better way to determine Posm is by simply finding the average force that the membrane exerts to hold impermeable molecules on the mixture side and then dividing that force by the membrane surface area: Fmem Posm = (6) Amem We use this method of finding Posm throughout this work. Theoretical justification for the use of Eq. (6) is based upon a simple equilibrium force balance argument: the force per unit area that artificially maintains the pressure gradient must clearly be equivalent to the osmotic pressure difference. In practice, the osmotic pressure difference is computed as the time-average of the instantaneous force per unit area exerted by the membrane. 2.6. Force calculations The molecular models used in this work include LJ and Coulumbic interactions. The LJ interactions and real-space contributions to the Coulumbic interactions were truncated at r cut = 10 Å. Long-range Coulumbic interactions were computed in reciprocal space using the particle–particle–particle–mesh (P3 M) methodology that is described in detail elsewhere [47–49]. Reciprocal space forces were determined
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
57
by ik differentiation (differentiation in Fourier space). The value of α, which determines the partition of the calculation between real and reciprocal space, was recalculated at each volume change according to the prescription given by Deserno and Holm [49] in order to optimise calculation accuracy. The P3 M optimal influence function was also updated at each volume change. 2.7. Fast equilibration Previously, the primary difficulty associated with the OMD method has been the difficulty of establishing chemical potential equilibrium across the membrane due to slow mass-transfer rates. In this work, it was found that the problem of slow mass-transfer could be greatly mitigated by pre-emptively changing the volume (and, thus, the pressure) on the mixture side (Side B) to null out the net migration of solvent. In essence, we apply the appropriate osmotic pressure to null the molecular flux instead of allowing it to build the osmotic pressure difference as has been done previously. Net mass-transfer of the solvent from the pure solvent side (Side A) to Side B, or visa versa, is prohibited by continually re-scaling the z-directional velocities of the solvent molecules so that the net migration across the membrane is zero. However, the net flux is determined from those velocities, and the chemical potential driving force producing the net migration is eliminated by increasing or decreasing the Side B volume in proportion to the observed flux. In this way, chemical potential equilibrium across the membrane is achieved, while net mass-transfer is prohibited. Once equilibrium is achieved, net mass-transfer becomes zero, so the prohibition of net mass-transfer can be discontinued. However, as one might suspect, we have found that continuing to re-zero tiny random fluctuations in net mass-transfer after equilibrium is achieved does not affect the simulation results. This method is analogous to periodically stopping a ball from rolling down an incline while adjusting the angle of inclination to remove the driving force. The ball is briefly allowed to roll, its velocity is measured, the angle is adjusted in proportion to the measured velocity, and the procedure is repeated until the ball no longer rolls when released. If we want the ball to stop rolling (establish equilibrium), we must eliminate both the velocity and the acceleration. We eliminate the velocity by periodically stopping the ball and we eliminate the acceleration by adjusting the angle of inclination. Establishment of equilibrium would be much more difficult if we were to let the ball continue to roll while adjusting the angle of inclination. Likewise, establishment of equilibrium in OMD is facilitated by re-zeroing the net flux of the permeable component while seeking to null the chemical potential driving force. After the establishment of equilibrium, any additional adjustments that are based on random noise are inconsequential. During fast equilibration, the net flux is determined in the following manner. The right half of Side A is grouped with the left half of Side B into “Region 1,” and the left half of Side A is grouped with the right half of Side B into “Region 2” (see Fig. 1). Then, at each time step, the z-directional velocities of permeable molecules in Region 2 are summed and subtracted from the sum of the z-directional velocities of permeable molecules in Region 1. This sum is divided by the number of permeable molecules to produce the net velocity per permeable molecule, φ. k vz,k − l vz,l φ= (7) Npermeable where Npermeable is the total number of permeable molecules in the simulation, k the summation index for permeable molecules in Region 1, and l is the summation index for permeable molecules in
58
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
Fig. 1. The OMD simulation cell is divided by a semi-permeable membrane into Side A (left of the central membrane) and Side B (right of the central membrane). Permeable molecule z velocities from Region 1 (grey area) are added to the negative of the permeable molecule z velocities from Region 2 (white area) to determine the net flux of solvent. The z-direction box length of Side B, LzB , is gradually adjusted until the net flux is brought to zero. Meanwhile, flux is artificially prohibited by re-scaling permeable molecule z velocities.
Region 2. The value of φ is subtracted from the z-directional velocities of all permeable molecules in the simulation cell. By adjusting these velocities, net mass-transfer of solvent between Sides A and B is artificially prohibited, but the chemical potential driving force that produced the flux still must be nulled. Before being reset for the next time step, the recorded value of φ, which is proportional to the instantaneous flux, is added to the time-averaged flux accumulator, Φ, which continues accumulation for an arbitrary, pre-determined length of time (a period of 100 time steps was used in this work). A net chemical potential driving force can be detected above the natural, noisy, random fluctuations of the flux only if a very large sample is obtained. Time-averaging helps, but is insufficient. A sufficiently-strong flux signal can be obtained by performing multiple simulations in parallel and averaging the value of Φ between them. Details about this parallel processing technique are given in Section 2.8. At the end of the accumulation cycle, the grand average for Φ is used to determine the new value for LzB , the cell length of side B, according to the proportional control scheme: LzBnew = LzBold (1 + K )
(8)
where K is an arbitrary gain constant that was chosen to be 0.05 s/m for this work. Whenever LzB is changed, the z-positions of the molecules on Side B are scaled proportionally to fit the new LzB . By changing LzB , the volume on Side B of the simulation cell is adjusted, which in turn increases or decreases the pressure on Side B, and eliminates the net chemical potential driving force. In practice, very gradual and continuous changes of LzB (∼0.01% every 100 time steps) achieve chemical potential equilibrium after only several thousand time steps. In addition to dramatic reductions in equilibration time, this technique allows the user to hold the mixture–side composition constant during the simulation so that equilibrium results are produced at the desired composition and the pressure difference across the membrane remains as the only variable affecting the calculation of the activity coefficient. It is evident when osmotic equilibrium between Sides A and B has been established by the convergence of the pressure difference across the membrane to a constant value (see Fig. 2).
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
59
Fig. 2. Time progression of a typical OMD simulation (Run 10 of Section 2.9). From top to bottom, the ratio of the Side B volume to the Side A volume (left axis), the osmotic pressure difference (right axis), and the mixture mole fraction (left axis) are shown. Chemical potential equilibrium is achieved within roughly 20,000 time steps, as evidenced by the constant values of the plotted variables.
2.8. System size considerations and parallel processing Ideally, one would perform OMD simulations using a system of only a few hundred molecules in order to minimise CPU time requirements. Unfortunately, large instantaneous fluctuations of pressures and mole fractions require the use of a much larger system. One might argue that a small system run over a large amount of time would be preferable to a large system run over a small amount of time. However, a large spatial sample of the desired phase space is much better than a large temporal sample due to the long-time persistence of local anomalies in the local fluid structure. Also, a large system can be broken down into many small sub-systems, or members of the same canonical (NVT) ensemble, and run simultaneously on a parallel processor or on a bank of networked machines. We have chosen a small system size of 200–1000 molecules for each independent ensemble member in the simulations that we have performed. Depending upon the desired accuracy, we have used 10–200 members (Nmembers ) of the ensemble that can be performed in parallel. The parallel simulation procedure that we have adopted is as follows. 1. Each member of the canonical ensemble is moved forward a fixed amount of time (tcycle ), independent of the other ensemble members. Each ensemble member shares a common N, V, T, and simulation cell geometry (cell length in the x- and y-directions, L; Side A cell length in the z-direction, LzA ; and Side B cell length in the z-direction, LzB ).
60
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
2. Once each of the Nmembers of the ensemble reaches tcycle , each of the needed properties, Posm , xi , and the flux (Φ), is computed as the ensemble average. 3. At this point, L is scaled to control the pure solvent density on Side A to the desired liquid density, and LzB is scaled according to Eq. (8). 4. The new L and the new LzB are applied to each of the Nmembers of the ensemble, molecular positions are re-scaled to the new box geometry, and a new cycle begins. Chemical potential equilibrium is achieved once the variables of interest reach constant values, such as clearly occurs after 20,000 time steps in Fig. 2. Property collection for activity coefficient determination can then begin. 2.9. Methodology validation The pure solvent provides an ideal validation of the methodology since activity coefficients produced when all molecules are identical must be unity at all compositions. In our benchmark runs, identical molecules were tagged as either solvent or solute molecules, and simulations were performed at various compositions. We performed a set of simulations using the methodology described in the earlier sections to predict the activity coefficients across the composition range for a pure, simple non-polar fluid of LJ spheres similar to liquid methane. To save computer time for this test case, a radial cutoff for the LJ interactions was set at 6 Å and the intermolecular LJ potential was shifted up in energy so that there was no discontinuity at the cutoff. The LJ σ was set at 3.73 Å, and the LJ ε/k was set at 147.96 K. A set of 200 ensemble members was moved forward through time according to the scheme outlined in Section 2.8, with each ensemble member consisting of 1000 LJ spheres. The simulations were performed at 111.7 K and a pure solvent density of 22 mol/l, with a time step size of 3 fs. All simulations were run for at least 20,000 time steps before property collection began. Each data point required approximately 2 days of CPU time on one CPU of an SGI Origin 2000 supercomputer. Results of this test are given in Table 1. The simulations yielded an activity coefficient prediction very near unity for all 10 of the independent test cases. Wilson equation parameters that were regressed to fit this data are Λ12 = Λ21 = 0.999. Infinite-dilution activity coefficients calculated using the resultant Table 1 OMD test case simulation results Run 1 2 3 4 5 6 7 8 9 10
1 RT
P0 +P
x1
Posm (MPa)
−ln x1
−
0.9512 0.9029 0.8534 0.8054 0.7558 0.7084 0.6594 0.6104 0.5622 0.5141
1.03 2.11 3.23 4.46 5.75 7.18 8.63 10.25 11.98 13.90
0.0500 0.1022 0.1585 0.2164 0.2800 0.3448 0.4165 0.4936 0.5759 0.6654
−0.0502 −0.1027 −0.1570 −0.2167 −0.2786 −0.3470 −0.4162 −0.4930 −0.5746 −0.6646
P0
V˜1 dP
ln γ 1 −0.0001 −0.0005 0.0016 −0.0003 0.0014 −0.0022 0.0003 0.0005 0.0013 0.0007
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
61
parameters are both 1.003. These numbers suggest the level of precision that can be expected in activity coefficient prediction using the system parameters and methodology outlined in the earlier sections. A perfect prediction would have given ln γi = 0 for all ten runs shown in Table 1, Λ12 = Λ21 = 1, and infinite-dilution activity coefficients of unity. 3. OMD results The OMD method has previously been used to predict activity coefficients for mixtures of LJ spheres [30] and for mixtures of structured molecules made up of homogeneous LJ spheres [43]. In this work, we establish the utility of the OMD method for the prediction of activity coefficients for model fluids containing heterogeneous site interactions and polar groups represented by distributed point charges. We begin by comparing OMD-calculated activity coefficients with GEMC-derived activity coefficients using the same molecular models. We then apply the OMD method to the prediction of activity coefficients for six industrially-significant binary mixtures and compare the results to experimentally-measured activity coefficients for those mixtures. Next, we show that modification of key cross-parameters in the molecular models for the acetone/chloroform mixture can greatly enhance agreement with experiment, suggesting the need for improved molecular models. Finally, we demonstrate the ability to predict LLE by using OMD to calculate the LLE phase dome for the n-hexane/methanol binary. 3.1. Comparison with GEMC calculations L´ısal et al. [50] simulate phase equilibrium for complex fluid mixtures using a modified GEMC approach, which they call reaction GEMC or RGEMC. They simulate phase equilibrium for binaries involving isobutene, methanol, and methyl tert-butyl ether by incorporating pure component vapour pressure data into the phase equilibrium calculations for the mixtures. Here we compare their results for methanol and isobutene with OMD simulations. We performed the OMD simulations for the methanol/isobutene binary at 323.15 K using the molecular models given by L´ısal et al. [50]. In order to directly compare with OMD-predicted activity coefficients, P–x–y data from Table 6 of their work were used to calculate activity coefficients according to the VLE equations: γi =
Pyi , Pi,sat xi
i = 1, 2
(9)
where Psat indicates the pure component vapour pressure [51] at the system temperature. These activity coefficients, along with OMD-calculated and experimentally-measured activity coefficients for the methanol/isobutene binary are plotted in Fig. 3. Although both simulation methods yield good qualitative agreement with the experimental values, the two simulation methods do not agree quantitatively. Discrepancies between experimentally-determined and simulation-calculated activity coefficients are easily explained in terms of molecular model inadequacies. However, since the same molecular models were used in both simulations, discrepancies between the OMD and RGEMC calculations are indicative of differences between the two methodologies. In an effort to resolve the observed differences, we tested the results for thermodynamic consistency by the integral and derivative tests. Thermodynamic consistency tests enable examination of activity
62
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
Fig. 3. Comparison between OMD-calculated (䊏, 䊉), RGEMC-calculated (×), and experimentally-measured (䊐, 䊊) activity coefficients for the methanol (1)/isobutene (2) binary at 323.15 K.
coefficient data sets for self-consistency, independent of the molecular model used. To check for internal thermodynamic consistency, the predicted activity coefficients for each component were fit as a function of composition using a three-parameter Redlich–Kister expansion: ln γ1 = (x1 − 1)2 (A − B + C + 4Bx1 − 8Cx1 + 12Cx21 )
(10)
ln γ2 = x12 (A − 3B + 5C + 4Bx1 − 16Cx1 + 12Cx21 )
(11)
where A, B, and C are the variables to be regressed for each component of each set of activity coefficient data. Once fit to this model, each set of data was tested for thermodynamic consistency by the integral test: M − N × 100 (12) % error = M +N where
x
M= N=
ln
γ1 dx γ2
(13)
ln
γ1 dx γ2
(14)
0
x
1
and x is defined as the liquid mole fraction at which γ1 =0 ln γ2
(15)
When subjected to this test, the experimentally-measured data set yielded 1.8% error, the OMD-derived
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
63
data set yielded 2.8% error, and the RGEMC-derived data set yielded 97.8% error. The extremely high error for the RGEMC-derived data set is partially due to the fact that activity coefficient data was not available across the entire composition range. The other binary pairs reported by L´ısal et al. [50] similarly failed the integral test (>40% error), even though their coverage of the composition range was better. The derivative test for thermodynamic consistency was also performed at individual compositions rather than across the entire composition range, and therefore, provided a more definitive test for incomplete data sets. Percent error for the derivative test was defined using the same expression as for the integral test (Eq. (12)), except with the following definitions for M and N: M = x1
d ln γ1 dx1
(16)
N = x2
d ln γ2 dx2
(17)
Again using the data smoothed with Eqs. (10) and (11), the experimental data yielded an average absolute value of error of 5.2% with a maximum value of 17.3% for the derivative test. The OMD-derived data yielded 3.3% average absolute value of error and 8.7% maximum error. Over the range where activity coefficient data was available for the RGEMC data set (x1 < 0.3192), the average absolute value of error for the derivative test was 53% with a maximum value of 78.5%. These results and those shown in Fig. 3 suggest that the OMD method results agree with the RGEMC results within the latter method’s uncertainties as evidenced by its poor thermodynamic consistency between γ 1 and γ 2 . On the other hand, we take the relatively good thermodynamic consistency of the γ i values obtained from the OMD method as strong validation of its accuracy. 3.2. Six representative binary mixtures We use the OMD method to predict activity coefficients for model fluids representing six classes of liquid binary mixtures that have industrial significance: methanol/n-hexane, n-hexane/n-pentane, methanol/water, chloroform/acetone, n-hexane/chloroform, and methanol/chloroform. The represented classes of binaries are A/NP, NP/NP, A/A, NA/NA, NP/NA, and A/NA, respectively where NP is a non-polar component, A an associating component, and NA is a non-associating component. Molecular models for methanol [52], water [53], chloroform [54], acetone [54], n-hexane [55], and n-pentane [55] were obtained from the literature. Values for the model parameters are given in Table 2. The repulsion and dispersion potentials were represented in the models with united-atom, pair-wise additive, site–site, LJ potentials located at all heavy-atom nuclei. Lorenz–Berthelot (LB) combining rules were used to obtain values for all interactions between heterogeneous interaction sites. Point charges located at nuclear centres were used to model the permanent dipole interactions for water, methanol, chloroform, and acetone. The method of auxiliary bond constraints was used to model the planar acetone molecule [56]. In all cases, bond lengths and angles were fixed at their equilibrium positions. Molecular models for alkanes generally include torsional potentials that affect the angle between two sites in the alkane chain that have two other sites between them. The torsional potentials for n-hexane and n-pentane were modelled as a function of dihedral angle N with the Ryckaert–Bellemans cosine
64
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
Table 2 Model fluid parameters Group
σ (Å)
ε/k (K)
q (e)
Intermolecular interactions –O (acetone) –CH3 (acetone) –C (acetone) –CH (chloroform) –Cl (chloroform) –CH3 (n-hexane) –CH2 (n-hexane) –CH3 (methanol) –O (methanol) –H (methanol) –CH3 (n-pentane) –CH2 (n-pentane) –O (water) –H (water)
2.960 3.910 3.750 3.800 3.470 3.905 3.905 3.740 3.030 0.000 3.905 3.905 3.169 0.000
105.750 80.570 52.870 40.285 151.068 88.070 59.380 105.200 86.500 0.000 88.070 59.380 78.178 0.000
−0.4240 +0.0620 +0.3000 +0.4200 −0.1400 0.0000 0.0000 +0.2650 −0.7000 +0.4350 0.0000 0.0000 −0.8476 +0.4238
a0
a1
a2
a3
a4
a5
−1578
−368
3156
−3788
Torsional potential (n-hexane and n-pentane) 1116 1462
series [57]: utors (φ) = ai cosi (φ) k i=0 5
(18)
and the coefficients shown in Table 2. All three dihedral angles in n-hexane and both dihedral angles in n-pentane were assumed to be equivalent and given by the Ryckaert–Bellemans values. Once the OMD simulation was properly equilibrated and the raw data subsequently collected, the activity coefficient of the solvent at the mixture composition was calculated according to Eq. (1). Each simulation produced a single value of the activity coefficient for the solvent at a specific composition. We have performed eight OMD simulations for each component of the six binary mixtures. Activity coefficients independently obtained from the simulations for both components along the entire composition range were then correlated by regressing parameters in a smoothing equation to best fit the OMD simulation data. Here we use the Wilson equations as the smoothing equations. This procedure correlates simultaneously the activity coefficients for both components in a manner that rigorously demands perfect thermodynamic consistency, where the unfit data sets have already shown reasonably good consistency. Figs. 4–9 show OMD-predicted and experimentally-measured activity coefficients for the six systems investigated in this study. In each case, ln γ 1 and ln γ 2 are plotted versus mole fraction of component 1. OMD-predicted activity coefficients closely follow experimentally-measured activity coefficients for the methanol/n-hexane system and for the n-pentane/n-hexane system, and the prediction method clearly distinguishes between the nearly ideal n-pentane/n-hexane system and the highly non-ideal
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
65
Fig. 4. Activity coefficient comparison for the methanol (1)/n-hexane (2) system at 35 ◦ C. Open symbols represent experimentally-measured activity coefficients [63], solid symbols represent OMD calculations, and the lines represent the OMD data fitted to Wilson equations.
Fig. 5. Activity coefficient comparison for the n-pentane (1)/n-hexane (2) system at 25 ◦ C [64]. Symbols are the same as in Fig. 4. Please note the small scale which makes the agreement appear to be worse than it is.
66
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
Fig. 6. Activity coefficient comparison for the acetone (1)/chloroform (2) system at 25 ◦ C [65]. Symbols are the same as in Fig. 4. The key cross-parameter, (ε/k)CH–O = 65 K, was set according to the LB rule.
methanol/n-hexane system. This is significant since no empirical data about the binary mixtures was used to perform this calculation. The results are entirely predictive, using only the simple molecular models and the LB combining rule. Agreement between simulated and experimental results vary for the six systems from excellent for methanol/n-hexane and n-pentane/n-hexane to good for chloroform/n-hexane and chloroform/methanol
Fig. 7. Activity coefficient comparison for the chloroform (1)/n-hexane (2) system at 25 ◦ C [66]. Symbols are the same as in Fig. 4.
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
67
Fig. 8. Activity coefficient comparison for the chloroform (1)/methanol (2) system at 35 ◦ C [67]. Symbols are the same as in Fig. 4.
to poor for acetone/chloroform and methanol/water. The agreement between the OMD values (solid points) and the Wilson equation fits (solid lines) indicates the good thermodynamic consistency of the data. The thermodynamic consistency of the OMD predictions, along with the methodology validation shown in Section 2.9 suggests that the disagreement between experimentally-measured activity coefficients and OMD-predicted activity coefficients lies in the molecular models. This is not surprising, considering
Fig. 9. Activity coefficient comparison for the methanol (1)/water (2) system at 25 ◦ C [68]. Symbols are the same as in Fig. 4.
68
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
the simplistic empirical models that were used. Especially suspect are the LB-derived cross-parameters that were used to model the cross-interactions between heterogeneous pairs of interacting sites. The inadequacy of the LB combining rules has been shown on numerous occasions [58,59]. 3.3. Cross-parameter modification for the acetone/chloroform binary The sensitivity of the simulated values to the LB assumption was examined for the acetone/chloroform system. As this system represents the classic example of a binary in which hydrogen-bonding occurs in the mixture but not in the pure components, it is not surprising that the LB combining rules fail to give even the correct sign for ln γ i in Fig. 6. This hydrogen bonding between the oxygen atom of the acetone and the hydrogen atom of the chloroform molecule explains the negative deviation from Roault’s law that occurs in the acetone/chloroform binary. The acetone/chloroform activity coefficients predicted by the OMD simulations shown in Fig. 6 show positive deviation from Roault’s law, suggesting that the hydrogen bonding effect was dramatically underrepresented by the molecular model crossparameters. As the results are expected to be most sensitive to this oxygen–hydrogen cross interaction, we have looked at two different values for the σ CH–O and εCH–O parameters. In the first case, we use ab initio calculations at the MP2/6-311+G (2df, 2pd) level to compute unbonded dimer energies between acetone and chloroform. Counterpoise-corrected, super-molecule energies obtained in a manner similar to that used by Rowley and co-workers [59,60] were then used to regress the values of σCH–O = 3.00 Å and (ε/k)CH–O = 690 K while using the LB combining rule for all other cross-interactions. The results shown in Fig. 10 are seen to correctly show negative deviations from Roault’s law, but this value of
Fig. 10. Activity coefficient comparison for the acetone (1)/chloroform (2) system at 25 ◦ C. Open symbols represent experimentally-measured activity coefficients [65], solid symbols and solid lines represent OMD calculations and the corresponding Wilson equations fits, respectively for the set of simulations where (ε/k)CH–O was set at 228 K. The dashed lines represent the OMD data fitted to Wilson equations given (ε/k)CH–O = 687 K.
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
69
ε clearly over predicts the effect of hydrogen bonding. For comparison purposes, we have also picked an intermediate value for the εCH–O term. Shown in Fig. 10 are the simulation results when σCH–O = 3.38 Å and (ε/k)CH–O = 228 K are used. As can be seen in Fig. 10, this value of εCH–O produces activity coefficients in good agreement with experimental measurements. Obviously, the OMD method is only as accurate as the interaction models used to represent the mixture. These results suggest that the main limitation in using OMD simulations to predict activity coefficients is the parameterisation of the cross-interactions. Too little attention has been given to this aspect of model development in the literature with too much reliance on combining rules such as LB. 3.4. LLE prediction Perhaps even more challenging than prediction of vapour–liquid equilibrium (VLE) is the prediction of liquid–liquid equilibrium (LLE). For LLE prediction, insertion methods would require the insertion of molecules into two dense liquid phases, further exacerbating insertion difficulties. Since activity coefficients are predicted directly by OMD simulations, and the fast equilibration method eliminates the need for molecules from one-phase to naturally insert through diffusion, no such problem arises when applying OMD to LLE prediction. We have performed OMD simulations to predict LLE for the methanol/n-hexane binary at five temperatures: −20, 0, 20, 40, and 60 ◦ C. The molecular models given in Table 2, along with the LB cross-parameters, were used. Since the Wilson equations are inadequate for LLE prediction, we used the T–K–Wilson equations to smooth the resultant simulated activity coefficients. The OMD-predicted activity coefficients along with the T–K–Wilson fit are shown in Fig. 11 for the −20 ◦ C case. Although there is no apparent problem in doing OMD simulations in the two-phase region (cf. Fig. 4), here we
Fig. 11. OMD-calculated activity coefficient prediction for the methanol (1)/n-hexane (2) system at −20 ◦ C. Symbols are the same as in Fig. 4.
70
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
Fig. 12. Methanol activity plotted against n-hexane activity for the methanol (1)/n-hexane (2) system at −20 ◦ C. The intersection of the plot with itself indicates the presence of two liquid phases in equilibrium.
derive the phase equilibrium from simulations performed only at compositions near pure component to emphasise that only a few simulated points in the one-phase region need be made. Given the activity coefficients of both components across the composition range at the given temperature, the equilibrium compositions of the coexisting liquid phases can be found graphically (as the self-crossing point) in Fig. 12, or alternatively, by solving the isoactivity equations simultaneously: xi γi = xi γi ,
i = 1, 2
(19)
where the primes indicate the two separate liquid phases. Equilibrium compositions were found for all five temperatures and are compared with experimentally-measured LLE data [61] for the methanol/n-hexane binary in Fig. 13. As can be seen, the LLE compositions are in reasonably good agreement with the measured solubility data when using the same molecular model that produced good VLE predictions (Fig. 4). A well-known difficulty with empirical models is their inadequacy in representing LLE and VLE with the same model and/or model parameters [62]. In this respect, the simulations outperform empirical equations. The consolute temperature is not, however, correctly predicted because of its sensitivity to the intermolecular potential model and because neither classical empirical models nor our finite-domain
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
71
Fig. 13. OMD-predicted LLE (䊉) compared with three sets of experimental LLE data (䉱, 䉬, ×) [61] for the methanol/n-hexane binary.
simulations can produce the correct critical exponent in the near-critical region where the correlation length diverges rapidly. 4. Conclusion The OMD method has been refined and applied to the calculation of activity coefficients for complex model chemical mixtures. The methodology has been shown to produce thermodynamically-consistent results, and agrees with the GEMC method of VLE prediction within the accuracy limitations suggested by the thermodynamic consistency test of the GEMC results. OMD-calculated activity coefficients were compared with experimentally-measured activity coefficients for six representative binary mixtures, yielding good results for the systems expected to be insensitive to values of cross-parameters in the potential model. Results for the acetone/chloroform system disagreed even in the sign of ln γ i . However, modification of the cross-interaction energy between the hydrogen bonding atoms brought the simulated results into good agreement with experimental values. LLE predictions were also made using the OMD method. Despite difficulties associated with the prediction of the consolute temperature, reasonably good agreement with experimental solubility data was achieved. The OMD method enables activity coefficient calculation by computer simulation without the need for costly molecule insertions. Activity coefficients can be computed at state points other than phase boundaries with good internal thermodynamic consistency. The method is especially amenable to parallelisation. Improvement of molecular models, especially cross-interactions for binary mixtures, should be pursued to further the effort to accurately predict activity coefficients by computer simulation.
72
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
References [1] D. Frenkel, in: G. Giccotti, W.G. Hoover (Eds.), Molecular Dynamics Simulation of Statistical Mechanical Systems, North-Holland, Amsterdam, 1986, p. 151. [2] M. Mezei, Mol. Phys. 61 (1987) 565. [3] M. Mezei, Mol. Phys. 87 (1989) 4881. [4] B. Guillot, Y. Guissani, Mol. Phys. 54 (1985) 455. [5] J.M. Haile, Fluid Phase Equil. 26 (1986) 103. [6] A.A. Chialvo, J. Chem. Phys. 92 (1990) 673. [7] B. Widom, J. Chem. Phys. 39 (1963) 2808. [8] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1987. [9] G.L. Deitrick, L.E. Scriven, H.T. Davis, J. Chem. Phys. 90 (1989) 2370. [10] Y. Tamai, H. Tanaka, K. Nakanishi, Fluid Phase Equil. 104 (1995) 363. [11] A. Pohorille, M.A. Wilson, J. Chem. Phys. 104 (1996) 3760. [12] P. Jedlovszky, M. Mezei, J. Am. Chem. Soc. 122 (2000) 5125. [13] H. Meirovitch, J. Phys. A 15 (1982) 735. [14] H.C. Oettinger, Macromolecules 18 (1985) 93. [15] S. Kumar, I. Szleifer, A.Z. Panagiotopoulos, Phys. Rev. Lett. 66 (1991) 2935. [16] J.J. de Pablo, M. Laso, U.W. Suter, J. Chem. Phys. 96 (1992) 6157. [17] D. Frenkel, B. Smit, Mol. Phys. 75 (1992) 983. [18] Y. Sheng, A.Z. Panagiotopoulos, D.P. Tassios, AIChE J. 41 (1995) 2306. [19] M. Müller, W. Paul, J. Chem. Phys. 100 (1994) 719. [20] N.B. Wilding, M. Müller, J. Chem. Phys. 101 (1994) 4324. [21] M. Wolfgardt, J. Baschnagel, K. Binder, J. Chem. Phys. 103 (1995) 7166. [22] J.R. Errington, A.Z. Panagiotopoulos, J. Chem. Phys. 111 (1999) 9731. [23] K.S. Shing, K.E. Gubbins, Mol. Phys. 46 (1982) 1109. [24] C.H. Bennett, J. Comput. Phys. 22 (1976) 245. [25] A.M. Ferrenberg, R.H. Swendsen, Phys. Rev. Lett. 61 (1988) 2635. [26] A.M. Ferrenberg, R.H. Swendsen, Phys. Rev. Lett. 63 (1989) 1195. [27] K. Kiyohara, K.E. Gubbins, A.Z. Panagiotopoulos, J. Chem. Phys. 106 (1997) 3338. [28] J.G. Powles, B. Holtz, W.A.B. Evans, J. Chem. Phys. 101 (1994) 7804. [29] R.L. Rowley, T.D. Shupe, M.W. Shuck, Mol. Phys. 82 (1994) 841. [30] R.L. Rowley, M.W. Shuck, J.C. Perry, Mol. Phys. 86 (1995) 125. [31] R.L. Rowley, M. Henrichsen, Fluid Phase Equil. 137 (1997) 75. [32] P. Bryk, A. Patrykiejew, O. Pizio, S. Sokolowski, Mol. Phys. 90 (1997) 483. [33] D.A. Kofke, P.T. Cummings, Mol. Phys. 92 (1997) 973. [34] D.A. Kofke, P.T. Cummings, Fluid Phase Equil. 151 (1998) 41. [35] K.S. Shing, K.E. Gubbins, Mol. Phys. 46 (1982) 1109. [36] K.-K. Han, J.H. Cushman, D.J. Diestler, J. Chem. Phys. 93 (1990) 5167. [37] N.G. Parsonage, J. Chem. Soc. Faraday Trans. 91 (1995) 2971. [38] N.G. Parsonage, J. Chem. Soc. Faraday Trans. 92 (1996) 1129. [39] N.G. Parsonage, Mol. Phys. 89 (1996) 1133. [40] G.C. Boulougouris, I.G. Economou, D.N. Theodorou, Mol. Phys. 96 (1999) 905. [41] P.S. Crozier, R.L. Rowley, D. Henderson, J. Chem. Phys. 113 (2000) 9202. [42] P.S. Crozier, R.L. Rowley, N.B. Holladay, D. Henderson, D.D. Busath, Phys. Rev. Lett. 86 (2001) 2467. [43] E.S. Eriksen, Activity Coefficients and Excess Properties of Non-Polar Binary Mixtures as Predicted by Osmotic Molecular Dynamics, Brigham Young University, Provo, UT, 1999, pp. 31–35. [44] R. Edberg, D.J. Evans, G.P. Morriss, J. Chem. Phys. 84 (1986) 6933. [45] R.L. Rowley, J.F. Ely, Mol. Phys. 72 (1981) 831. [46] M. Henrichsen, Application of Osmotic Molecular Dynamics to Structured Molecules, Brigham Young University, Provo, UT, 1996, pp. 31–34. [47] R.W. Hockney, J.W. Eastwood, Computer Simulation Using Particles, Institute of Physics, Bristol, 1988.
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73 [48] [49] [50] [51] [52] [53] [54] [55] [56] [57] [58] [59] [60] [61] [62] [63] [64] [65] [66] [67] [68]
73
M. Deserno, C. Holm, J. Chem. Phys. 109 (1998) 7678. M. Deserno, C. Holm, J. Chem. Phys. 109 (1998) 7694. M. L´ısal, W.R. Smith, I. Nezbeda, J. Phys. Chem. B 103 (1999) 10496. S. Ung, M.F. Doherty, Chem. Eng. Sci. 50 (1995) 23. M.E. van Leeuwen, B. Smit, J. Phys. Chem. 99 (1995) 1831. H.J. Berendsen, J.R. Grigera, T.P. Straatsma, J. Phys. Chem. 91 (1987) 6269. W.L. Jorgensen, J.M. Briggs, M.L. Contreras, J. Phys. Chem. 94 (1990) 1683. W.L. Jorgensen, J.D. Madura, C.J. Swenson, J. Am. Chem. Soc. 106 (1984) 6638. G. Ciccotti, M. Ferrario, J.-P. Ryckaert, Mol. Phys. 47 (1982) 1253. J.P. Ryckaert, A. Bellemans, Chem. Phys. Lett. 30 (1975) 123. J. Delhommelle, P. Millié, Mol. Phys. 99 (2001) 619. R.L. Rowley, T. Pakkanen, J. Chem. Phys. 110 (1999) 3368. R.L. Rowley, Y. Yang, T. Pakkanen, J. Chem. Phys. 114 (2001) 6058. J.M. Sørensen, W. Arlt, in: D. Behrens, R. Eckermann (Eds.), Liquid–liquid Equilibrium Data Collection, Binary Systems (Part 1), Vol. V, DECHEMA, Frankfurt, Germany, 1979, p. 92. J.M. Prausnitz, R.N. Lichtenthaler, E.G. de Azevedo, Molecular Thermodynamics of Fluid Phase Equilibria, Third Edition, Prentice-Hall, Upper Saddle River, NJ, 1999, p. 795. J. Gmehling, U. Onken, in: D. Behrens, R. Eckermann (Eds.), Vapour–Liquid Equilibrium Data Collection, Organic Hydroxy Compounds: Alcohols (Part 2a), Vol. I, DECHEMA, Frankfurt, Germany, 1977, p. 255. J. Gmehling, U. Onken, W. Arlt, in: D. Behrens, R. Eckermann (Eds.), Vapour–Liquid Equilibrium Data Collection, Aliphatic Hydrocarbons, C4 –C6 (Part 6a), Vol. I, DECHEMA, Frankfurt, Germany, 1980, p. 122. J. Gmehling, U. Onken, J.R. Rarey, in: R. Eckermann, G. Kreysa (Eds.), Vapour–Liquid Equilibrium Data Collection, Ketones (Part 3b), Vol. I (Suppl. 1), DECHEMA, Frankfurt, Germany, 1993, p. 12. J. Gmehling, U. Onken, W. Arlt, in: D. Behrens, R. Eckermann (Eds.), Vapour–Liquid Equilibrium Data Collection, Aliphatic Hydrocarbons, C4 –C6 (Part 6a), Vol. I, DECHEMA, Frankfurt, Germany, 1980, p. 423. J. Gmehling, U. Onken, in: D. Behrens, R. Eckermann (Eds.), Vapour–Liquid Equilibrium Data Collection, Organic Hydroxy Compounds: Alcohols (Part 2a), Vol. I, DECHEMA, Frankfurt, Germany, 1977, p. 17. J. Gmehling, U. Onken, J.R. Rarey-Nies, in: D. Behrens, R. Eckermann (Eds.), Vapour–Liquid Equilibrium Data Collection, Aqueous Systems (Part 1b), Vol. I (Suppl. 2), DECHEMA, Frankfurt, Germany, 1988, p. 29.
Activity coefficient prediction by osmotic molecular dynamics Paul S. Crozier∗ , Richard L. Rowley Department of Chemical Engineering, Brigham Young University, Provo, UT 84602, USA Received 6 March 2001; accepted 24 September 2001
Abstract The osmotic molecular dynamics method (OMD) is used to calculate activity coefficients for vapour–liquid equilibria (VLE) and liquid–liquid equilibria (LLE) predictions. The previously reported OMD methodology is refined and applied to mixtures of polar, structured molecular fluids. Other computer simulation approaches to phase equilibria prediction are discussed briefly, and comparison to recent Gibbs Ensemble Monte Carlo (GEMC) results is made. OMD-predicted activity coefficients are compared to experimentally-measured activity coefficients for six industrially-significant binary mixtures (methanol/n-hexane, n-hexane/n-pentane, chloroform/acetone, chloroform/methanol, methanol/water, chloroform/n-hexane). Molecular model inadequacies, especially cross-parameters between unlike molecules, are shown. A single cross-parameter for the acetone/chloroform binary is modified to produce good agreement with experimentally-measured activity coefficients. Also, OMD-derived LLE predictions are produced for the methanol/n-hexane system and compared with experimentally-measured LLE data. © 2002 Elsevier Science B.V. All rights reserved. Keywords: Computer simulation; Methods of calculation; Molecular interactions; Vapour–liquid equilibria; Liquid–liquid equilibria; Chemical potential
1. Introduction For decades, scientists and engineers have sought to predict activity coefficients and excess thermodynamic properties based upon molecular-level chemical characteristics. As computer speed and availability continue to increase, predicting thermodynamic properties by simulating substances on the molecular-level has gained acceptance. Thermodynamic integration [1–6], particle insertion [7,8], particle growing [9–22], distribution-histogram methods [23–27], and osmotic methods [28–32] are some of the more common methods currently used for the calculation of chemical potential by computer simulation. Kofke and Cummings [33,34] offer excellent reviews of the various methods and explain the inherent ∗
Corresponding author. Present address: Sandia National Laboratiories, P.O. Box 5800, MS 0316 Albuquerque, NM 87185-0316, USA. Tel.: +1-505-845-9714; fax: +1-505-284-5451. E-mail address: [email protected] (P.S. Crozier). 0378-3812/02/$ – see front matter © 2002 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 8 - 3 8 1 2 ( 0 1 ) 0 0 7 3 4 - 8
54
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
error in particle deletion schemes [35–40]. Insertion methods are widely used, the Gibbs Ensemble Monte Carlo (GEMC) method being the most widely accepted for VLE prediction. In this study, we refine the osmotic molecular dynamics method (OMD) [29–31] that can be used to calculate chemical potential and phase equilibrium. Unlike GEMC methods, OMD does not require costly particle insertions, and is not limited to the calculation of chemical potentials at phase boundaries. Using the same molecular models, we compare OMD results with recent GEMC phase equilibrium predictions. In addition, we compare OMD activity coefficient predictions to experimentally-measured activity coefficients for VLE and LLE systems. 2. OMD methodology 2.1. OMD thermodynamics As discussed in [29–31], activity coefficients can be determined by osmosis experiments in which a pure liquid solvent is separated from a liquid mixture of solvent and solute by a semi-permeable membrane, permeable only to solvent molecules. The activity coefficient is calculated according to the thermodynamic identity: 1 P0 +Posm ˜ ln γi = −ln xi − Vi dP (1) RT P0 where γ i is the activity coefficient of the permeable component (solvent) at a mole fraction xi in the mixture, P0 the pressure on the pure solvent side of the membrane, and P0 + P osm is the pressure on the liquid mixture side of the membrane. The molar volume of solvent, V˜i , is left inside the integral and considered to be a function of pressure. This molar volume versus pressure relationship can be determined by simulation as described in [30]. 2.2. Membrane interactions We used an artificial fixed-in-space membrane that is permeable only to solvent molecules and of the same functionality as that given by Rowley et al. [29]. Rowley et al. used the OMD method to calculate activity coefficients for Lennard–Jones (LJ) fluid mixtures, and modelled the interaction between the membrane and solute particles with one-dimensional, LJ potentials truncated at the potential-well minimum and shifted up in energy by ε, the LJ energy interaction parameter. This provided a continuous, purely repulsive potential based on the distance of the particle from the membrane. Their method was modified in this study to accommodate heterogeneous interaction sites on larger structured molecules. Each site i on solute (impermeable) molecules was set to interact with the membrane according to the following force formula: σi,i 6 σi,i 6 Fi (z) = 24εi,i (2) − 1 , for z < 21/6 σi,i 2 z z and Fi (z) = 0,
for z ≥ 21/6 σi,i
(3)
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
55
where z is given by z = 21/6 σi,i − d
(4)
and d is the perpendicular distance between the molecule’s centre of mass and the membrane surface. The value of d is defined to be positive only after the centre of mass of the solvent molecule has passed through the membrane, and is negative otherwise (no interaction). Although the value of the activity coefficient obtained from OMD simulations is independent of the form of the purely repulsive potential chosen for the solute-membrane interactions (the membrane must simply be semi-permeable), using a soft repulsive force like Eq. (2) avoids impulses that can occur with harder potentials during the finite time step associated with the numerical integration of the equations of motion. The form of Eq. (4) is chosen so that a molecule uniformly feels the repulsion when it crosses the membrane surface. The repulsion is exerted on all sites simultaneously, regardless of the molecule’s orientation. 2.3. Slot definition The simulation cell must be of adequate length in the z-direction (perpendicular to the membrane surface) to allow the formation of a fluid region that is isolated from, and independent of, the artificial membrane. Even though the chosen membrane described in the earlier sections is relatively undisruptive of the natural liquid structure, a volume of fluid within a fixed distance from the membrane is excluded from property calculations to ensure that only fluid whose structure is not perturbed by the membrane surface is used to calculate the activity coefficient. But it is desirable to choose small values for the z-direction box dimensions in order to minimise the time required for diffusion and equilibration. For the purposes of this work, the central “slot” region of unperturbed fluid was defined as all molecules with centre of mass at least 10 Å away from the membrane surface. Simulation studies [41,42] of net-neutral systems that include long-range Coulumbic interactions, similar to the present study, suggest that a 10 Å interfacial layer adequately shields even long-range interactions. Benchmark tests, discussed in the following sections, confirm that this distance is sufficient for the molecules of this study. Smaller distances were found to introduce error into the calculation. 2.4. Thermostatting Precise temperature control of the system is crucial to obtaining accurate activity coefficient predictions. Thermostatting of the entire simulation cell to the set-point temperature does not guarantee that the individual slot regions will be held at precisely the desired temperature during mass-transfer [43], especially since the membrane artificially adds energy to the system in the excluded volume regions. We independently thermostat the slot and adjacent exclusion volumes of each half-cell to ensure uniform temperature at the set-point throughout. The population of molecules pertaining to any given slot or exclusion volume is updated every 100 time steps, at which point the net momentum of that population of molecules in each direction is set to zero by velocity re-scaling, followed by an additional velocity re-scaling that fixes the population at the set-point temperature. A Gaussian thermostat [44,45] is used to hold each population’s temperature constant between updates. The zeroing of net momenta ensures that a net drift of the population that can accumulate in OMD simulations does not produce an erroneously high temperature measurement [46].
56
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
Temperature is calculated from the translational motion using 2 i mi vi T = k(3N − NRDF)
(5)
where T is the temperature, mi and v i are the mass and velocity of molecule i, k the Boltzmann constant, N the number of molecules in the population for a molecular thermostat or the number of atoms for an atomistic thermostat, and NRDF is the number of removed degrees of freedom (d.f.). One d.f. is removed due to the thermostat. In a regular MD simulation, momentum is conserved in all three directions for the thermostatted population; so three additional d.f. must be removed. However, in an OMD simulation, the thermostatted population is continually changing; new populations have not been previously thermostatted to the set-point temperature. Also, momentum is not conserved for the thermostatted populations since the molecules of each slot and exclusion volume interact with the fixed, rigid membrane and/or molecules outside that region. No d.f. were removed due to the bond and angle constraints because molecular temperatures, rather than atomic temperatures, were computed. Hence, for the simulations of this work, NRDF = 1 in Eq. (5). 2.5. Osmotic pressure calculation The system temperature has a strong effect on the system pressure and on the osmotic pressure difference. Prediction of activity coefficients by Eq. (1) is extremely sensitive to the precise determination of the equilibrium osmotic pressure difference, Posm . Although this pressure difference can be found by taking the difference between the pressure found on the mixture side of the membrane (Side A) and the pressure found on the solvent side of the membrane (Side B), each being calculated using the virial equation, this procedure is computationally expensive and yields large instantaneous fluctuations in the calculation of Posm . We have found that this method of finding Posm is especially prone to error if discontinuous molecular models (e.g. cutoff fluids) are used. A better way to determine Posm is by simply finding the average force that the membrane exerts to hold impermeable molecules on the mixture side and then dividing that force by the membrane surface area: Fmem Posm = (6) Amem We use this method of finding Posm throughout this work. Theoretical justification for the use of Eq. (6) is based upon a simple equilibrium force balance argument: the force per unit area that artificially maintains the pressure gradient must clearly be equivalent to the osmotic pressure difference. In practice, the osmotic pressure difference is computed as the time-average of the instantaneous force per unit area exerted by the membrane. 2.6. Force calculations The molecular models used in this work include LJ and Coulumbic interactions. The LJ interactions and real-space contributions to the Coulumbic interactions were truncated at r cut = 10 Å. Long-range Coulumbic interactions were computed in reciprocal space using the particle–particle–particle–mesh (P3 M) methodology that is described in detail elsewhere [47–49]. Reciprocal space forces were determined
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
57
by ik differentiation (differentiation in Fourier space). The value of α, which determines the partition of the calculation between real and reciprocal space, was recalculated at each volume change according to the prescription given by Deserno and Holm [49] in order to optimise calculation accuracy. The P3 M optimal influence function was also updated at each volume change. 2.7. Fast equilibration Previously, the primary difficulty associated with the OMD method has been the difficulty of establishing chemical potential equilibrium across the membrane due to slow mass-transfer rates. In this work, it was found that the problem of slow mass-transfer could be greatly mitigated by pre-emptively changing the volume (and, thus, the pressure) on the mixture side (Side B) to null out the net migration of solvent. In essence, we apply the appropriate osmotic pressure to null the molecular flux instead of allowing it to build the osmotic pressure difference as has been done previously. Net mass-transfer of the solvent from the pure solvent side (Side A) to Side B, or visa versa, is prohibited by continually re-scaling the z-directional velocities of the solvent molecules so that the net migration across the membrane is zero. However, the net flux is determined from those velocities, and the chemical potential driving force producing the net migration is eliminated by increasing or decreasing the Side B volume in proportion to the observed flux. In this way, chemical potential equilibrium across the membrane is achieved, while net mass-transfer is prohibited. Once equilibrium is achieved, net mass-transfer becomes zero, so the prohibition of net mass-transfer can be discontinued. However, as one might suspect, we have found that continuing to re-zero tiny random fluctuations in net mass-transfer after equilibrium is achieved does not affect the simulation results. This method is analogous to periodically stopping a ball from rolling down an incline while adjusting the angle of inclination to remove the driving force. The ball is briefly allowed to roll, its velocity is measured, the angle is adjusted in proportion to the measured velocity, and the procedure is repeated until the ball no longer rolls when released. If we want the ball to stop rolling (establish equilibrium), we must eliminate both the velocity and the acceleration. We eliminate the velocity by periodically stopping the ball and we eliminate the acceleration by adjusting the angle of inclination. Establishment of equilibrium would be much more difficult if we were to let the ball continue to roll while adjusting the angle of inclination. Likewise, establishment of equilibrium in OMD is facilitated by re-zeroing the net flux of the permeable component while seeking to null the chemical potential driving force. After the establishment of equilibrium, any additional adjustments that are based on random noise are inconsequential. During fast equilibration, the net flux is determined in the following manner. The right half of Side A is grouped with the left half of Side B into “Region 1,” and the left half of Side A is grouped with the right half of Side B into “Region 2” (see Fig. 1). Then, at each time step, the z-directional velocities of permeable molecules in Region 2 are summed and subtracted from the sum of the z-directional velocities of permeable molecules in Region 1. This sum is divided by the number of permeable molecules to produce the net velocity per permeable molecule, φ. k vz,k − l vz,l φ= (7) Npermeable where Npermeable is the total number of permeable molecules in the simulation, k the summation index for permeable molecules in Region 1, and l is the summation index for permeable molecules in
58
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
Fig. 1. The OMD simulation cell is divided by a semi-permeable membrane into Side A (left of the central membrane) and Side B (right of the central membrane). Permeable molecule z velocities from Region 1 (grey area) are added to the negative of the permeable molecule z velocities from Region 2 (white area) to determine the net flux of solvent. The z-direction box length of Side B, LzB , is gradually adjusted until the net flux is brought to zero. Meanwhile, flux is artificially prohibited by re-scaling permeable molecule z velocities.
Region 2. The value of φ is subtracted from the z-directional velocities of all permeable molecules in the simulation cell. By adjusting these velocities, net mass-transfer of solvent between Sides A and B is artificially prohibited, but the chemical potential driving force that produced the flux still must be nulled. Before being reset for the next time step, the recorded value of φ, which is proportional to the instantaneous flux, is added to the time-averaged flux accumulator, Φ, which continues accumulation for an arbitrary, pre-determined length of time (a period of 100 time steps was used in this work). A net chemical potential driving force can be detected above the natural, noisy, random fluctuations of the flux only if a very large sample is obtained. Time-averaging helps, but is insufficient. A sufficiently-strong flux signal can be obtained by performing multiple simulations in parallel and averaging the value of Φ between them. Details about this parallel processing technique are given in Section 2.8. At the end of the accumulation cycle, the grand average for Φ is used to determine the new value for LzB , the cell length of side B, according to the proportional control scheme: LzBnew = LzBold (1 + K )
(8)
where K is an arbitrary gain constant that was chosen to be 0.05 s/m for this work. Whenever LzB is changed, the z-positions of the molecules on Side B are scaled proportionally to fit the new LzB . By changing LzB , the volume on Side B of the simulation cell is adjusted, which in turn increases or decreases the pressure on Side B, and eliminates the net chemical potential driving force. In practice, very gradual and continuous changes of LzB (∼0.01% every 100 time steps) achieve chemical potential equilibrium after only several thousand time steps. In addition to dramatic reductions in equilibration time, this technique allows the user to hold the mixture–side composition constant during the simulation so that equilibrium results are produced at the desired composition and the pressure difference across the membrane remains as the only variable affecting the calculation of the activity coefficient. It is evident when osmotic equilibrium between Sides A and B has been established by the convergence of the pressure difference across the membrane to a constant value (see Fig. 2).
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
59
Fig. 2. Time progression of a typical OMD simulation (Run 10 of Section 2.9). From top to bottom, the ratio of the Side B volume to the Side A volume (left axis), the osmotic pressure difference (right axis), and the mixture mole fraction (left axis) are shown. Chemical potential equilibrium is achieved within roughly 20,000 time steps, as evidenced by the constant values of the plotted variables.
2.8. System size considerations and parallel processing Ideally, one would perform OMD simulations using a system of only a few hundred molecules in order to minimise CPU time requirements. Unfortunately, large instantaneous fluctuations of pressures and mole fractions require the use of a much larger system. One might argue that a small system run over a large amount of time would be preferable to a large system run over a small amount of time. However, a large spatial sample of the desired phase space is much better than a large temporal sample due to the long-time persistence of local anomalies in the local fluid structure. Also, a large system can be broken down into many small sub-systems, or members of the same canonical (NVT) ensemble, and run simultaneously on a parallel processor or on a bank of networked machines. We have chosen a small system size of 200–1000 molecules for each independent ensemble member in the simulations that we have performed. Depending upon the desired accuracy, we have used 10–200 members (Nmembers ) of the ensemble that can be performed in parallel. The parallel simulation procedure that we have adopted is as follows. 1. Each member of the canonical ensemble is moved forward a fixed amount of time (tcycle ), independent of the other ensemble members. Each ensemble member shares a common N, V, T, and simulation cell geometry (cell length in the x- and y-directions, L; Side A cell length in the z-direction, LzA ; and Side B cell length in the z-direction, LzB ).
60
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
2. Once each of the Nmembers of the ensemble reaches tcycle , each of the needed properties, Posm , xi , and the flux (Φ), is computed as the ensemble average. 3. At this point, L is scaled to control the pure solvent density on Side A to the desired liquid density, and LzB is scaled according to Eq. (8). 4. The new L and the new LzB are applied to each of the Nmembers of the ensemble, molecular positions are re-scaled to the new box geometry, and a new cycle begins. Chemical potential equilibrium is achieved once the variables of interest reach constant values, such as clearly occurs after 20,000 time steps in Fig. 2. Property collection for activity coefficient determination can then begin. 2.9. Methodology validation The pure solvent provides an ideal validation of the methodology since activity coefficients produced when all molecules are identical must be unity at all compositions. In our benchmark runs, identical molecules were tagged as either solvent or solute molecules, and simulations were performed at various compositions. We performed a set of simulations using the methodology described in the earlier sections to predict the activity coefficients across the composition range for a pure, simple non-polar fluid of LJ spheres similar to liquid methane. To save computer time for this test case, a radial cutoff for the LJ interactions was set at 6 Å and the intermolecular LJ potential was shifted up in energy so that there was no discontinuity at the cutoff. The LJ σ was set at 3.73 Å, and the LJ ε/k was set at 147.96 K. A set of 200 ensemble members was moved forward through time according to the scheme outlined in Section 2.8, with each ensemble member consisting of 1000 LJ spheres. The simulations were performed at 111.7 K and a pure solvent density of 22 mol/l, with a time step size of 3 fs. All simulations were run for at least 20,000 time steps before property collection began. Each data point required approximately 2 days of CPU time on one CPU of an SGI Origin 2000 supercomputer. Results of this test are given in Table 1. The simulations yielded an activity coefficient prediction very near unity for all 10 of the independent test cases. Wilson equation parameters that were regressed to fit this data are Λ12 = Λ21 = 0.999. Infinite-dilution activity coefficients calculated using the resultant Table 1 OMD test case simulation results Run 1 2 3 4 5 6 7 8 9 10
1 RT
P0 +P
x1
Posm (MPa)
−ln x1
−
0.9512 0.9029 0.8534 0.8054 0.7558 0.7084 0.6594 0.6104 0.5622 0.5141
1.03 2.11 3.23 4.46 5.75 7.18 8.63 10.25 11.98 13.90
0.0500 0.1022 0.1585 0.2164 0.2800 0.3448 0.4165 0.4936 0.5759 0.6654
−0.0502 −0.1027 −0.1570 −0.2167 −0.2786 −0.3470 −0.4162 −0.4930 −0.5746 −0.6646
P0
V˜1 dP
ln γ 1 −0.0001 −0.0005 0.0016 −0.0003 0.0014 −0.0022 0.0003 0.0005 0.0013 0.0007
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
61
parameters are both 1.003. These numbers suggest the level of precision that can be expected in activity coefficient prediction using the system parameters and methodology outlined in the earlier sections. A perfect prediction would have given ln γi = 0 for all ten runs shown in Table 1, Λ12 = Λ21 = 1, and infinite-dilution activity coefficients of unity. 3. OMD results The OMD method has previously been used to predict activity coefficients for mixtures of LJ spheres [30] and for mixtures of structured molecules made up of homogeneous LJ spheres [43]. In this work, we establish the utility of the OMD method for the prediction of activity coefficients for model fluids containing heterogeneous site interactions and polar groups represented by distributed point charges. We begin by comparing OMD-calculated activity coefficients with GEMC-derived activity coefficients using the same molecular models. We then apply the OMD method to the prediction of activity coefficients for six industrially-significant binary mixtures and compare the results to experimentally-measured activity coefficients for those mixtures. Next, we show that modification of key cross-parameters in the molecular models for the acetone/chloroform mixture can greatly enhance agreement with experiment, suggesting the need for improved molecular models. Finally, we demonstrate the ability to predict LLE by using OMD to calculate the LLE phase dome for the n-hexane/methanol binary. 3.1. Comparison with GEMC calculations L´ısal et al. [50] simulate phase equilibrium for complex fluid mixtures using a modified GEMC approach, which they call reaction GEMC or RGEMC. They simulate phase equilibrium for binaries involving isobutene, methanol, and methyl tert-butyl ether by incorporating pure component vapour pressure data into the phase equilibrium calculations for the mixtures. Here we compare their results for methanol and isobutene with OMD simulations. We performed the OMD simulations for the methanol/isobutene binary at 323.15 K using the molecular models given by L´ısal et al. [50]. In order to directly compare with OMD-predicted activity coefficients, P–x–y data from Table 6 of their work were used to calculate activity coefficients according to the VLE equations: γi =
Pyi , Pi,sat xi
i = 1, 2
(9)
where Psat indicates the pure component vapour pressure [51] at the system temperature. These activity coefficients, along with OMD-calculated and experimentally-measured activity coefficients for the methanol/isobutene binary are plotted in Fig. 3. Although both simulation methods yield good qualitative agreement with the experimental values, the two simulation methods do not agree quantitatively. Discrepancies between experimentally-determined and simulation-calculated activity coefficients are easily explained in terms of molecular model inadequacies. However, since the same molecular models were used in both simulations, discrepancies between the OMD and RGEMC calculations are indicative of differences between the two methodologies. In an effort to resolve the observed differences, we tested the results for thermodynamic consistency by the integral and derivative tests. Thermodynamic consistency tests enable examination of activity
62
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
Fig. 3. Comparison between OMD-calculated (䊏, 䊉), RGEMC-calculated (×), and experimentally-measured (䊐, 䊊) activity coefficients for the methanol (1)/isobutene (2) binary at 323.15 K.
coefficient data sets for self-consistency, independent of the molecular model used. To check for internal thermodynamic consistency, the predicted activity coefficients for each component were fit as a function of composition using a three-parameter Redlich–Kister expansion: ln γ1 = (x1 − 1)2 (A − B + C + 4Bx1 − 8Cx1 + 12Cx21 )
(10)
ln γ2 = x12 (A − 3B + 5C + 4Bx1 − 16Cx1 + 12Cx21 )
(11)
where A, B, and C are the variables to be regressed for each component of each set of activity coefficient data. Once fit to this model, each set of data was tested for thermodynamic consistency by the integral test: M − N × 100 (12) % error = M +N where
x
M= N=
ln
γ1 dx γ2
(13)
ln
γ1 dx γ2
(14)
0
x
1
and x is defined as the liquid mole fraction at which γ1 =0 ln γ2
(15)
When subjected to this test, the experimentally-measured data set yielded 1.8% error, the OMD-derived
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
63
data set yielded 2.8% error, and the RGEMC-derived data set yielded 97.8% error. The extremely high error for the RGEMC-derived data set is partially due to the fact that activity coefficient data was not available across the entire composition range. The other binary pairs reported by L´ısal et al. [50] similarly failed the integral test (>40% error), even though their coverage of the composition range was better. The derivative test for thermodynamic consistency was also performed at individual compositions rather than across the entire composition range, and therefore, provided a more definitive test for incomplete data sets. Percent error for the derivative test was defined using the same expression as for the integral test (Eq. (12)), except with the following definitions for M and N: M = x1
d ln γ1 dx1
(16)
N = x2
d ln γ2 dx2
(17)
Again using the data smoothed with Eqs. (10) and (11), the experimental data yielded an average absolute value of error of 5.2% with a maximum value of 17.3% for the derivative test. The OMD-derived data yielded 3.3% average absolute value of error and 8.7% maximum error. Over the range where activity coefficient data was available for the RGEMC data set (x1 < 0.3192), the average absolute value of error for the derivative test was 53% with a maximum value of 78.5%. These results and those shown in Fig. 3 suggest that the OMD method results agree with the RGEMC results within the latter method’s uncertainties as evidenced by its poor thermodynamic consistency between γ 1 and γ 2 . On the other hand, we take the relatively good thermodynamic consistency of the γ i values obtained from the OMD method as strong validation of its accuracy. 3.2. Six representative binary mixtures We use the OMD method to predict activity coefficients for model fluids representing six classes of liquid binary mixtures that have industrial significance: methanol/n-hexane, n-hexane/n-pentane, methanol/water, chloroform/acetone, n-hexane/chloroform, and methanol/chloroform. The represented classes of binaries are A/NP, NP/NP, A/A, NA/NA, NP/NA, and A/NA, respectively where NP is a non-polar component, A an associating component, and NA is a non-associating component. Molecular models for methanol [52], water [53], chloroform [54], acetone [54], n-hexane [55], and n-pentane [55] were obtained from the literature. Values for the model parameters are given in Table 2. The repulsion and dispersion potentials were represented in the models with united-atom, pair-wise additive, site–site, LJ potentials located at all heavy-atom nuclei. Lorenz–Berthelot (LB) combining rules were used to obtain values for all interactions between heterogeneous interaction sites. Point charges located at nuclear centres were used to model the permanent dipole interactions for water, methanol, chloroform, and acetone. The method of auxiliary bond constraints was used to model the planar acetone molecule [56]. In all cases, bond lengths and angles were fixed at their equilibrium positions. Molecular models for alkanes generally include torsional potentials that affect the angle between two sites in the alkane chain that have two other sites between them. The torsional potentials for n-hexane and n-pentane were modelled as a function of dihedral angle N with the Ryckaert–Bellemans cosine
64
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
Table 2 Model fluid parameters Group
σ (Å)
ε/k (K)
q (e)
Intermolecular interactions –O (acetone) –CH3 (acetone) –C (acetone) –CH (chloroform) –Cl (chloroform) –CH3 (n-hexane) –CH2 (n-hexane) –CH3 (methanol) –O (methanol) –H (methanol) –CH3 (n-pentane) –CH2 (n-pentane) –O (water) –H (water)
2.960 3.910 3.750 3.800 3.470 3.905 3.905 3.740 3.030 0.000 3.905 3.905 3.169 0.000
105.750 80.570 52.870 40.285 151.068 88.070 59.380 105.200 86.500 0.000 88.070 59.380 78.178 0.000
−0.4240 +0.0620 +0.3000 +0.4200 −0.1400 0.0000 0.0000 +0.2650 −0.7000 +0.4350 0.0000 0.0000 −0.8476 +0.4238
a0
a1
a2
a3
a4
a5
−1578
−368
3156
−3788
Torsional potential (n-hexane and n-pentane) 1116 1462
series [57]: utors (φ) = ai cosi (φ) k i=0 5
(18)
and the coefficients shown in Table 2. All three dihedral angles in n-hexane and both dihedral angles in n-pentane were assumed to be equivalent and given by the Ryckaert–Bellemans values. Once the OMD simulation was properly equilibrated and the raw data subsequently collected, the activity coefficient of the solvent at the mixture composition was calculated according to Eq. (1). Each simulation produced a single value of the activity coefficient for the solvent at a specific composition. We have performed eight OMD simulations for each component of the six binary mixtures. Activity coefficients independently obtained from the simulations for both components along the entire composition range were then correlated by regressing parameters in a smoothing equation to best fit the OMD simulation data. Here we use the Wilson equations as the smoothing equations. This procedure correlates simultaneously the activity coefficients for both components in a manner that rigorously demands perfect thermodynamic consistency, where the unfit data sets have already shown reasonably good consistency. Figs. 4–9 show OMD-predicted and experimentally-measured activity coefficients for the six systems investigated in this study. In each case, ln γ 1 and ln γ 2 are plotted versus mole fraction of component 1. OMD-predicted activity coefficients closely follow experimentally-measured activity coefficients for the methanol/n-hexane system and for the n-pentane/n-hexane system, and the prediction method clearly distinguishes between the nearly ideal n-pentane/n-hexane system and the highly non-ideal
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
65
Fig. 4. Activity coefficient comparison for the methanol (1)/n-hexane (2) system at 35 ◦ C. Open symbols represent experimentally-measured activity coefficients [63], solid symbols represent OMD calculations, and the lines represent the OMD data fitted to Wilson equations.
Fig. 5. Activity coefficient comparison for the n-pentane (1)/n-hexane (2) system at 25 ◦ C [64]. Symbols are the same as in Fig. 4. Please note the small scale which makes the agreement appear to be worse than it is.
66
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
Fig. 6. Activity coefficient comparison for the acetone (1)/chloroform (2) system at 25 ◦ C [65]. Symbols are the same as in Fig. 4. The key cross-parameter, (ε/k)CH–O = 65 K, was set according to the LB rule.
methanol/n-hexane system. This is significant since no empirical data about the binary mixtures was used to perform this calculation. The results are entirely predictive, using only the simple molecular models and the LB combining rule. Agreement between simulated and experimental results vary for the six systems from excellent for methanol/n-hexane and n-pentane/n-hexane to good for chloroform/n-hexane and chloroform/methanol
Fig. 7. Activity coefficient comparison for the chloroform (1)/n-hexane (2) system at 25 ◦ C [66]. Symbols are the same as in Fig. 4.
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
67
Fig. 8. Activity coefficient comparison for the chloroform (1)/methanol (2) system at 35 ◦ C [67]. Symbols are the same as in Fig. 4.
to poor for acetone/chloroform and methanol/water. The agreement between the OMD values (solid points) and the Wilson equation fits (solid lines) indicates the good thermodynamic consistency of the data. The thermodynamic consistency of the OMD predictions, along with the methodology validation shown in Section 2.9 suggests that the disagreement between experimentally-measured activity coefficients and OMD-predicted activity coefficients lies in the molecular models. This is not surprising, considering
Fig. 9. Activity coefficient comparison for the methanol (1)/water (2) system at 25 ◦ C [68]. Symbols are the same as in Fig. 4.
68
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
the simplistic empirical models that were used. Especially suspect are the LB-derived cross-parameters that were used to model the cross-interactions between heterogeneous pairs of interacting sites. The inadequacy of the LB combining rules has been shown on numerous occasions [58,59]. 3.3. Cross-parameter modification for the acetone/chloroform binary The sensitivity of the simulated values to the LB assumption was examined for the acetone/chloroform system. As this system represents the classic example of a binary in which hydrogen-bonding occurs in the mixture but not in the pure components, it is not surprising that the LB combining rules fail to give even the correct sign for ln γ i in Fig. 6. This hydrogen bonding between the oxygen atom of the acetone and the hydrogen atom of the chloroform molecule explains the negative deviation from Roault’s law that occurs in the acetone/chloroform binary. The acetone/chloroform activity coefficients predicted by the OMD simulations shown in Fig. 6 show positive deviation from Roault’s law, suggesting that the hydrogen bonding effect was dramatically underrepresented by the molecular model crossparameters. As the results are expected to be most sensitive to this oxygen–hydrogen cross interaction, we have looked at two different values for the σ CH–O and εCH–O parameters. In the first case, we use ab initio calculations at the MP2/6-311+G (2df, 2pd) level to compute unbonded dimer energies between acetone and chloroform. Counterpoise-corrected, super-molecule energies obtained in a manner similar to that used by Rowley and co-workers [59,60] were then used to regress the values of σCH–O = 3.00 Å and (ε/k)CH–O = 690 K while using the LB combining rule for all other cross-interactions. The results shown in Fig. 10 are seen to correctly show negative deviations from Roault’s law, but this value of
Fig. 10. Activity coefficient comparison for the acetone (1)/chloroform (2) system at 25 ◦ C. Open symbols represent experimentally-measured activity coefficients [65], solid symbols and solid lines represent OMD calculations and the corresponding Wilson equations fits, respectively for the set of simulations where (ε/k)CH–O was set at 228 K. The dashed lines represent the OMD data fitted to Wilson equations given (ε/k)CH–O = 687 K.
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
69
ε clearly over predicts the effect of hydrogen bonding. For comparison purposes, we have also picked an intermediate value for the εCH–O term. Shown in Fig. 10 are the simulation results when σCH–O = 3.38 Å and (ε/k)CH–O = 228 K are used. As can be seen in Fig. 10, this value of εCH–O produces activity coefficients in good agreement with experimental measurements. Obviously, the OMD method is only as accurate as the interaction models used to represent the mixture. These results suggest that the main limitation in using OMD simulations to predict activity coefficients is the parameterisation of the cross-interactions. Too little attention has been given to this aspect of model development in the literature with too much reliance on combining rules such as LB. 3.4. LLE prediction Perhaps even more challenging than prediction of vapour–liquid equilibrium (VLE) is the prediction of liquid–liquid equilibrium (LLE). For LLE prediction, insertion methods would require the insertion of molecules into two dense liquid phases, further exacerbating insertion difficulties. Since activity coefficients are predicted directly by OMD simulations, and the fast equilibration method eliminates the need for molecules from one-phase to naturally insert through diffusion, no such problem arises when applying OMD to LLE prediction. We have performed OMD simulations to predict LLE for the methanol/n-hexane binary at five temperatures: −20, 0, 20, 40, and 60 ◦ C. The molecular models given in Table 2, along with the LB cross-parameters, were used. Since the Wilson equations are inadequate for LLE prediction, we used the T–K–Wilson equations to smooth the resultant simulated activity coefficients. The OMD-predicted activity coefficients along with the T–K–Wilson fit are shown in Fig. 11 for the −20 ◦ C case. Although there is no apparent problem in doing OMD simulations in the two-phase region (cf. Fig. 4), here we
Fig. 11. OMD-calculated activity coefficient prediction for the methanol (1)/n-hexane (2) system at −20 ◦ C. Symbols are the same as in Fig. 4.
70
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
Fig. 12. Methanol activity plotted against n-hexane activity for the methanol (1)/n-hexane (2) system at −20 ◦ C. The intersection of the plot with itself indicates the presence of two liquid phases in equilibrium.
derive the phase equilibrium from simulations performed only at compositions near pure component to emphasise that only a few simulated points in the one-phase region need be made. Given the activity coefficients of both components across the composition range at the given temperature, the equilibrium compositions of the coexisting liquid phases can be found graphically (as the self-crossing point) in Fig. 12, or alternatively, by solving the isoactivity equations simultaneously: xi γi = xi γi ,
i = 1, 2
(19)
where the primes indicate the two separate liquid phases. Equilibrium compositions were found for all five temperatures and are compared with experimentally-measured LLE data [61] for the methanol/n-hexane binary in Fig. 13. As can be seen, the LLE compositions are in reasonably good agreement with the measured solubility data when using the same molecular model that produced good VLE predictions (Fig. 4). A well-known difficulty with empirical models is their inadequacy in representing LLE and VLE with the same model and/or model parameters [62]. In this respect, the simulations outperform empirical equations. The consolute temperature is not, however, correctly predicted because of its sensitivity to the intermolecular potential model and because neither classical empirical models nor our finite-domain
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
71
Fig. 13. OMD-predicted LLE (䊉) compared with three sets of experimental LLE data (䉱, 䉬, ×) [61] for the methanol/n-hexane binary.
simulations can produce the correct critical exponent in the near-critical region where the correlation length diverges rapidly. 4. Conclusion The OMD method has been refined and applied to the calculation of activity coefficients for complex model chemical mixtures. The methodology has been shown to produce thermodynamically-consistent results, and agrees with the GEMC method of VLE prediction within the accuracy limitations suggested by the thermodynamic consistency test of the GEMC results. OMD-calculated activity coefficients were compared with experimentally-measured activity coefficients for six representative binary mixtures, yielding good results for the systems expected to be insensitive to values of cross-parameters in the potential model. Results for the acetone/chloroform system disagreed even in the sign of ln γ i . However, modification of the cross-interaction energy between the hydrogen bonding atoms brought the simulated results into good agreement with experimental values. LLE predictions were also made using the OMD method. Despite difficulties associated with the prediction of the consolute temperature, reasonably good agreement with experimental solubility data was achieved. The OMD method enables activity coefficient calculation by computer simulation without the need for costly molecule insertions. Activity coefficients can be computed at state points other than phase boundaries with good internal thermodynamic consistency. The method is especially amenable to parallelisation. Improvement of molecular models, especially cross-interactions for binary mixtures, should be pursued to further the effort to accurately predict activity coefficients by computer simulation.
72
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73
References [1] D. Frenkel, in: G. Giccotti, W.G. Hoover (Eds.), Molecular Dynamics Simulation of Statistical Mechanical Systems, North-Holland, Amsterdam, 1986, p. 151. [2] M. Mezei, Mol. Phys. 61 (1987) 565. [3] M. Mezei, Mol. Phys. 87 (1989) 4881. [4] B. Guillot, Y. Guissani, Mol. Phys. 54 (1985) 455. [5] J.M. Haile, Fluid Phase Equil. 26 (1986) 103. [6] A.A. Chialvo, J. Chem. Phys. 92 (1990) 673. [7] B. Widom, J. Chem. Phys. 39 (1963) 2808. [8] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1987. [9] G.L. Deitrick, L.E. Scriven, H.T. Davis, J. Chem. Phys. 90 (1989) 2370. [10] Y. Tamai, H. Tanaka, K. Nakanishi, Fluid Phase Equil. 104 (1995) 363. [11] A. Pohorille, M.A. Wilson, J. Chem. Phys. 104 (1996) 3760. [12] P. Jedlovszky, M. Mezei, J. Am. Chem. Soc. 122 (2000) 5125. [13] H. Meirovitch, J. Phys. A 15 (1982) 735. [14] H.C. Oettinger, Macromolecules 18 (1985) 93. [15] S. Kumar, I. Szleifer, A.Z. Panagiotopoulos, Phys. Rev. Lett. 66 (1991) 2935. [16] J.J. de Pablo, M. Laso, U.W. Suter, J. Chem. Phys. 96 (1992) 6157. [17] D. Frenkel, B. Smit, Mol. Phys. 75 (1992) 983. [18] Y. Sheng, A.Z. Panagiotopoulos, D.P. Tassios, AIChE J. 41 (1995) 2306. [19] M. Müller, W. Paul, J. Chem. Phys. 100 (1994) 719. [20] N.B. Wilding, M. Müller, J. Chem. Phys. 101 (1994) 4324. [21] M. Wolfgardt, J. Baschnagel, K. Binder, J. Chem. Phys. 103 (1995) 7166. [22] J.R. Errington, A.Z. Panagiotopoulos, J. Chem. Phys. 111 (1999) 9731. [23] K.S. Shing, K.E. Gubbins, Mol. Phys. 46 (1982) 1109. [24] C.H. Bennett, J. Comput. Phys. 22 (1976) 245. [25] A.M. Ferrenberg, R.H. Swendsen, Phys. Rev. Lett. 61 (1988) 2635. [26] A.M. Ferrenberg, R.H. Swendsen, Phys. Rev. Lett. 63 (1989) 1195. [27] K. Kiyohara, K.E. Gubbins, A.Z. Panagiotopoulos, J. Chem. Phys. 106 (1997) 3338. [28] J.G. Powles, B. Holtz, W.A.B. Evans, J. Chem. Phys. 101 (1994) 7804. [29] R.L. Rowley, T.D. Shupe, M.W. Shuck, Mol. Phys. 82 (1994) 841. [30] R.L. Rowley, M.W. Shuck, J.C. Perry, Mol. Phys. 86 (1995) 125. [31] R.L. Rowley, M. Henrichsen, Fluid Phase Equil. 137 (1997) 75. [32] P. Bryk, A. Patrykiejew, O. Pizio, S. Sokolowski, Mol. Phys. 90 (1997) 483. [33] D.A. Kofke, P.T. Cummings, Mol. Phys. 92 (1997) 973. [34] D.A. Kofke, P.T. Cummings, Fluid Phase Equil. 151 (1998) 41. [35] K.S. Shing, K.E. Gubbins, Mol. Phys. 46 (1982) 1109. [36] K.-K. Han, J.H. Cushman, D.J. Diestler, J. Chem. Phys. 93 (1990) 5167. [37] N.G. Parsonage, J. Chem. Soc. Faraday Trans. 91 (1995) 2971. [38] N.G. Parsonage, J. Chem. Soc. Faraday Trans. 92 (1996) 1129. [39] N.G. Parsonage, Mol. Phys. 89 (1996) 1133. [40] G.C. Boulougouris, I.G. Economou, D.N. Theodorou, Mol. Phys. 96 (1999) 905. [41] P.S. Crozier, R.L. Rowley, D. Henderson, J. Chem. Phys. 113 (2000) 9202. [42] P.S. Crozier, R.L. Rowley, N.B. Holladay, D. Henderson, D.D. Busath, Phys. Rev. Lett. 86 (2001) 2467. [43] E.S. Eriksen, Activity Coefficients and Excess Properties of Non-Polar Binary Mixtures as Predicted by Osmotic Molecular Dynamics, Brigham Young University, Provo, UT, 1999, pp. 31–35. [44] R. Edberg, D.J. Evans, G.P. Morriss, J. Chem. Phys. 84 (1986) 6933. [45] R.L. Rowley, J.F. Ely, Mol. Phys. 72 (1981) 831. [46] M. Henrichsen, Application of Osmotic Molecular Dynamics to Structured Molecules, Brigham Young University, Provo, UT, 1996, pp. 31–34. [47] R.W. Hockney, J.W. Eastwood, Computer Simulation Using Particles, Institute of Physics, Bristol, 1988.
P.S. Crozier, R.L. Rowley / Fluid Phase Equilibria 193 (2002) 53–73 [48] [49] [50] [51] [52] [53] [54] [55] [56] [57] [58] [59] [60] [61] [62] [63] [64] [65] [66] [67] [68]
73
M. Deserno, C. Holm, J. Chem. Phys. 109 (1998) 7678. M. Deserno, C. Holm, J. Chem. Phys. 109 (1998) 7694. M. L´ısal, W.R. Smith, I. Nezbeda, J. Phys. Chem. B 103 (1999) 10496. S. Ung, M.F. Doherty, Chem. Eng. Sci. 50 (1995) 23. M.E. van Leeuwen, B. Smit, J. Phys. Chem. 99 (1995) 1831. H.J. Berendsen, J.R. Grigera, T.P. Straatsma, J. Phys. Chem. 91 (1987) 6269. W.L. Jorgensen, J.M. Briggs, M.L. Contreras, J. Phys. Chem. 94 (1990) 1683. W.L. Jorgensen, J.D. Madura, C.J. Swenson, J. Am. Chem. Soc. 106 (1984) 6638. G. Ciccotti, M. Ferrario, J.-P. Ryckaert, Mol. Phys. 47 (1982) 1253. J.P. Ryckaert, A. Bellemans, Chem. Phys. Lett. 30 (1975) 123. J. Delhommelle, P. Millié, Mol. Phys. 99 (2001) 619. R.L. Rowley, T. Pakkanen, J. Chem. Phys. 110 (1999) 3368. R.L. Rowley, Y. Yang, T. Pakkanen, J. Chem. Phys. 114 (2001) 6058. J.M. Sørensen, W. Arlt, in: D. Behrens, R. Eckermann (Eds.), Liquid–liquid Equilibrium Data Collection, Binary Systems (Part 1), Vol. V, DECHEMA, Frankfurt, Germany, 1979, p. 92. J.M. Prausnitz, R.N. Lichtenthaler, E.G. de Azevedo, Molecular Thermodynamics of Fluid Phase Equilibria, Third Edition, Prentice-Hall, Upper Saddle River, NJ, 1999, p. 795. J. Gmehling, U. Onken, in: D. Behrens, R. Eckermann (Eds.), Vapour–Liquid Equilibrium Data Collection, Organic Hydroxy Compounds: Alcohols (Part 2a), Vol. I, DECHEMA, Frankfurt, Germany, 1977, p. 255. J. Gmehling, U. Onken, W. Arlt, in: D. Behrens, R. Eckermann (Eds.), Vapour–Liquid Equilibrium Data Collection, Aliphatic Hydrocarbons, C4 –C6 (Part 6a), Vol. I, DECHEMA, Frankfurt, Germany, 1980, p. 122. J. Gmehling, U. Onken, J.R. Rarey, in: R. Eckermann, G. Kreysa (Eds.), Vapour–Liquid Equilibrium Data Collection, Ketones (Part 3b), Vol. I (Suppl. 1), DECHEMA, Frankfurt, Germany, 1993, p. 12. J. Gmehling, U. Onken, W. Arlt, in: D. Behrens, R. Eckermann (Eds.), Vapour–Liquid Equilibrium Data Collection, Aliphatic Hydrocarbons, C4 –C6 (Part 6a), Vol. I, DECHEMA, Frankfurt, Germany, 1980, p. 423. J. Gmehling, U. Onken, in: D. Behrens, R. Eckermann (Eds.), Vapour–Liquid Equilibrium Data Collection, Organic Hydroxy Compounds: Alcohols (Part 2a), Vol. I, DECHEMA, Frankfurt, Germany, 1977, p. 17. J. Gmehling, U. Onken, J.R. Rarey-Nies, in: D. Behrens, R. Eckermann (Eds.), Vapour–Liquid Equilibrium Data Collection, Aqueous Systems (Part 1b), Vol. I (Suppl. 2), DECHEMA, Frankfurt, Germany, 1988, p. 29.