Molecular Modeling of Fluid-Phase Equilibria Using an Isotropic

Jul 19, 2003 - In a recent paper,1 we suggested the use of temperature-quench ..... the initial constant-volume quench produces a three-phase system, ...
0 downloads 0 Views 152KB Size
Download PDF
Ind. Eng. Chem. Res. 2003, 42, 4123-4131

4123

Molecular Modeling of Fluid-Phase Equilibria Using an Isotropic Multipolar Potential Erich A. Mu 1 ller* Departamento de Termodina´ mica y Feno´ menos de Transferencia, Universidad Simo´ n Bolı´var, Caracas 1080, Venezuela

Lev D. Gelb† Department of Chemistry, Washington University, St. Louis, Missouri 63130-1134

We apply a simplified intermolecular potential function to model the fluid-phase behavior of realistic multicomponent polar systems. The isotropic multipolar potential includes dispersion, repulsion, and spherically averaged multipolar contributions and has a simple isotropic mathematical form with only two adjustable parameters. These energy and size parameters may be fitted to experimental pure-component properties such as liquid densities or critical points. These potentials, coupled with large-scale temperature-quench molecular dynamics simulations, yield thermophysical data such as the temperature-density diagram of benzene and naphthalene, the high-pressure PvT properties of carbon dioxide, the vapor-liquid equilibria of 1,2-dichloroethane/cyclohexane, and the liquid-liquid equilibria of the sulfolane/n-octane/ benzene system. Despite the explicit omission of both molecular shape and directionality of multipolar interactions from the potential, the results are in quantitative agreement with experimental results, suggesting that it is the average energetic interactions that dictate the gross details of fluid phase equilibria. 1. Introduction In a recent paper,1 we suggested the use of temperature-quench molecular dynamics (TQMD) as a probe of multiphase equilibria. In TQMD, a homogeneous system is quenched at a constant volume to a temperature where the thermodynamically stable state consists of multiple phases. The system evolves rapidly into locally equilibrated subdomains, which may be explored by statistical means to obtain densities and/or compositions. The method requires the use of relatively large system sizes, usually of the O(105) particles, so that “bulklike” domains of each phase will be present. Simulations on this scale are time-consuming and may be accelerated by (a) the use of parallel processing (which is straightforward for molecular dynamics,2-4 and several such codes are available5-9) and/or (b) simplification of the intermolecular potential energy functions, the evaluation of which accounts for nearly all of the computer time used by a typical molecular simulation code. This paper focuses on the latter proposal, based on the hypothesis that it is the energetic contributions (as opposed to the entropic) that distinguish the differences in fluid-phase behavior between different polar systems and that detailed atomistic potential models may not be necessary to capture such behavior. For many fluids at relatively high temperatures (>100 K), rapid molecular reorientation and translation destroys both shortrange order and orientational correlation, so a detailed * To whom correspondence should be addressed. Tel.: +58 (212) 9063749. Fax: +58 (212) 9063743. E-mail: emuller@ usb.ve. † Tel.: (314) 935-5026. Fax: (314) 935-4481. E-mail: gelb@ wuchem.wustl.edu.

model may not be needed. We suggest that isotropic potentials may be sufficient under these conditions. In the first part of this paper, we present an overview and a sample derivation of an effective intermolecular potential of this type, and in the latter part, we employ it to model the phase equilibria of several multicomponent polar systems. 2. Intermolecular Potential for Polar Molecules A brief discussion is given here of commonly used semiempirical intermolecular functions for the calculation of fluid-phase thermodynamic and transport properties; however, the reader is referred to the standard literature10-12 and recent reviews13,14 for a more thorough discussion. 2.1. General Treatment. The simplest empirical potentials useful in modeling condensed phases are pairwise-additive isotropic central-force functions with both attractive and repulsive terms. One of the most commonly used is the Lennard-Jones (LJ) potential, which sums a soft repulsive term Urepulsion and an attractive (dispersion) term Udispersion, and is given by

[( ) ( ) ]

ULJ ij ) 4ij

σij rij

12

-

σij rij

6

(1)

where σ and , the size and energy parameters, may be regressed from experimental data, e.g., second virial coefficients or gas-phase viscosity data,15 or from a corresponding states principle fit to the critical temperature and density. The LJ model has been used extensively and successfully for evaluation of phase equilibrium and transport properties of both polar and nonpolar molecules. The model is simple and useful, but it does

10.1021/ie030033y CCC: $25.00 © 2003 American Chemical Society Published on Web 07/19/2003

4124

Ind. Eng. Chem. Res., Vol. 42, No. 17, 2003

not reproduce the energy curves of simple fluids, and it cannot treat in an explicit fashion (a) molecular shape, (b) charge distributions, (c) intramolecular interactions, or (d) second-order effects, such as induction or multibody interactions. More realistic potential models assign a spherical site to each atom or group of atoms in a molecule. Additionally, one may place partial point charges on selected locations:

Uij )

qiqj

+ Udispersion )+∑ ∑(Urepulsion ij ij 4π0rij

(2)

where the summations run over all possible pairs of sites on each molecule (noting that the charge sites need not coincide with molecular centers), qi is a partial charge assigned to a particular site, and 0 is the vacuum permittivity. Last, one may include bond bending, stretching, and vibration. The resulting classes of potentials are frequently used to represent realistic systems; examples are the CHARMM,16 AMBER,17 MMFF,18 and COMPASS19 force fields. Many of the necessary potential parameters may be fit to ab initio calculations [e.g., refs 20 and 21], though further tuning of dispersion terms to match experimental data is generally required. Highly detailed force fields are difficult to parametrize and computationally expensive to use. One may attempt to simplify the picture by lumping together atomic centers, in what are loosely known as “united atom” models. The success of these models in reproducing the phase behavior of alkanes22,23 is an example of how averaging some effects (such as induction forces and intramolecular vibrations) into simple “effective” potentials is sufficient for some applications. Consider the benzene molecule. A detailed potential might explicitly treat 12 atoms, corresponding bonds, and the electronic charge distribution. A recently proposed24 model unites each carbon atom with its corresponding hydrogen and averages electrostatic effects into six identical isotropic sites, one for each C-H group. The evaluation of the intermolecular potential among a single pair of benzene molecules is thus reduced from calculation of (12)2 atom-atom terms (plus charge-distribution interactions) to (6)2 terms, which, in turn, reduces by a factor of 4 the computational effort required for a simulation. Further reduction could be accomplished by reduction to a single spherical site; however, the resulting model is too simple to reproduce even the basic temperaturedensity diagram (cf. the dashed curve in Figure 4). A common approach to describing electrostatic interactions is to represent the charge distribution with a centrally located multipole expansion25,26

Uij ) Unonpolar + Upermanent moments + Uinduced moments µQ ) (Urepulsion + Udispersion ) + (Uµµ ij ij ij + Uij +

UQQ ij

+ ...) +

(UµR ij

+

UQR ij

+ ...) (3)

where the superscripts µµ, µQ, and QQ refer to the dipole-dipole, dipole-quadrupole, and quadrupolequadrupole interactions, respectively, while the induced moments µR and QR correspond to the dipole-induced dipole and quadrupole-induced quadrupole interactions, respectively, and R is the polarizability. With some exceptions, notably water, polarizability is nor-

mally neglected. For example, for an axially symmetric molecule such as benzene or carbon dioxide, the corre27 sponding quadrupole-quadrupole energy UQQ ij is

3QiQj QQ f (θi,θj,φij) 4rij5

UQQ ij )

(4)

where f QQ is a function of the molecular orientation (thus of the Euler angles θi and φi) of both interacting molecules, given by

f QQ ) 1 - 5 cos2 θi - 5 cos2 θj + 17 cos2 θi cos2 θj + 2 sin2 θi sin2 θj cos2 φij 16 sin θi cos θi sin θj cos θj cos φij (5) Some simplicity is gained by the multipolar expansion, but even when the attraction and repulsion of a molecule such as benzene is modeled with a single LJ site, the potential requires the calculation of the Euler angles θi and φi25 for each intermolecular interaction. 2.2. Isotropic Multipolar Potential (IMP). One can further simplify the intermolecular potentials by averaging over molecular orientations,28 yielding purely isotropic terms in the potential. This approach excludes orientational correlation from the treatment of a fluid, but at reasonably high temperatures, this may be an acceptable approximation. As an example, let us consider the quadrupolar contribution, given by eqs 4 and 5. For a single pair of molecules, the total potential energy Uij will be a function of both the distance and the orientation angles. The configurational integral, Z, for the pair interaction may be written as

Z)

V 8π

∫∫exp(-βUQQ ij ) dr dω

(6)

where rj is the relative coordinate vector, rj ) rjj - rji, and dω ) sin θi cos θj dθi dθj dφij, V is the system volume, and β ) 1/kT. If the potential function were isotropic, Γij, then the configurational integral would have a simpler form, i.e.

∫∫exp(-βΓQQ j ij ) dr

Z)V

(7)

Inspecting the last two equations, one finds that, to be equivalent,

ΓQQ ij ) -

∫

1 1 exp(-βUQQ ln ij ) dω β 8π

[

]

(8)

One may substitute eq 4 into eq 8 and expand the exponential into an infinite series (provided T is not too low and Q not too high) to obtain

ΓQQ ij ) -

[ (

1 1 ln β 8π

∫ 1-

) ]

3QiQj QQ Qi2Qj2 2 QQ2 βf + βf - ... dω 4rij5 4rij10 (9)

where upon performance of the integral, which is now independent of rij, becomes

ΓQQ ij

(

)

7Qi2Qj2 2 1 ) - ln 1 + β + ... β 5r 10 ij

(10)

Ind. Eng. Chem. Res., Vol. 42, No. 17, 2003 4125

which itself can be approximated by

7βQi2Qj2

ΓQQ ij ≈ -

(11)

5rij10

Similar procedures may be applied to the other types of multipolar contributions, leading to28 2

Γµµ ij ≈ 2

ΓµQ ij ≈ -

UKeesom ij

2

βµi µj

(12)

3rij6 2

2

(13)

2rij8

[( ) ( ) ] [ σij rij

12

β

-

σij rij

µi2µj2 3rij6

+

6

(µi2Rj + Riµj2)

(14)

6

rij

σij rij

6

µi2Qj2 + Qi2µj2

+

]

7Qi2Qj2 5rij10

(15)

For mixtures, the Lorentz-Berthelot mixing rules are usually employed to generate unlike-pair interactions; ij ) xij and σij ) (σi + σj)/2, where i and j refer to the pure-component values. Note that the multipolar terms require no empirical cross-parameters for mixtures, and if the Lorentz-Berthelot rules are used for the dispersion part of the potential, multicomponent phase equilibria calculations require no additional parameters other than pure-component quantities, i and σi. The potential, its parameters, and other macroscopic properties may be made dimensionless by taking U* ) U/1, r* ) r/σ1, (µ*)2 ) µ2/1σ13, (Q*)2 ) Q2/1σ1,5 R* ) R/σ13 , T* ) kT/1, and P* ) Pσ13/1 where the subscript 1 refers to one component of the mixture, or to the single component in the case of a simple fluid. It can be seen that the IMP is dependent only on the intermolecular distance rij; the last three terms of eq 15 are no longer functions of orientations. The IMP are not potential functions in the strict sense but rather an “effective” force field because they depend explicitly on temperature, an ensemble property. However, this rarely is a problem, as shall be seen later, for most simulations are performed under isothermal conditions, and because the quantities of interest are pressures, densities, and mole fractions, the nontrivial temperature derivatives of the free energy need not concern us. Equivalently, eqs 11-14 are the Helmholtz free energy of orientation for a molecular pair. Starting from eq 4, one may calculate the canonical ensemble angle-averaged energy,29 which gives a functional form similar to eq 11 (actually eq 11 multiplied

µi2µj2

-β

3rij6

)

[( ) ( ) ] σµµ rij

12

-

σµµ rij

6

(16)

where the temperature dependence of the modified energy and size parameters, µµ and σµµ, is given by

(σµµ)6 ) σ6/F, µµ ) F 2 F)1+

-

2rij8

[( ) ( ) ] 12

2

In summary, this IMP is the sum of isotropic dispersion and repulsion terms and equations of the type (11)(14). The actual number of terms taken will depend on the particular system. If we neglect polarizabilities, consider only dipolar and quadrupolar fluids, and use a LJ form for the nonpolar interaction, the IMP takes the form of

Uij ) 4ij

) 4ij

σij rij

4µµ

β(µi Qj + Qi µj )

ΓµR ij ≈ -

by a factor of 2). However, this and other approaches fail to give the correct configurational integral.28 The IMP may be interpreted as a modified LJ potential with temperature-dependent parameters. For a pure dipolar fluid, the IMP potential has been previously described12 as the Keesom potential, UKeesom , ij written25 as

βµi2µj2 12µµ(σµµ)6

(17) (18)

Also, if one is willing to go even further, one may fit the LJ parameters to functions of temperature as well. The approach works in an engineering sense;30 however, for higher multipoles, this mapping onto a LJ model is not possible because of the different scalings of the r-8 and r-10 terms in eq 15. The real incongruence of these approaches becomes apparent when simulating multicomponent systems. It is well-known that the formal condition for convergence of the multipolar expansion is that the distance between two interacting molecules be larger than the sum of the distances from the center of each molecule to the furthest part of its charge distribution. In practical terms, this convergence sphere may be larger than the van der Waals radii of molecules, especially for nonspherical molecules. This issue applies to both realistic potential models that explicitly treat electrostatics with truncated multipole expansions and IMP. In addition to this, the approximations and series expansions that produce the IMP equations are themselves only valid for relatively high temperatures and low values of the multipole moments. IMP should thus be understood as a rational way to obtain a functional form for an effective intermolecular potential for some species. An acceptable description of strongly polar molecules within this framework may require treatment of the multipole moments as adjustable parameters instead of use of experimental or ab initio values. Such a fitting is standard practice even for sophisticated models of polar fluids, e.g., the various potentials for liquid water.31 The general approach is expected to work in a satisfactory manner only if the density, mole fraction, or similar data are sought. Clearly, liquid crystal structures and other systems where orientational order is important cannot be modeled using IMPs. On the other hand, the simplicity of the potentials make them excellent for use in comprehensive phase equilibria studies. 3. Simulation Method 3.1. TQMD. The TQMD method for calculating phase equilibria was presented in a previous paper.1 The basis of the method is summarized in Figure 1. A purecomponent liquid-vapor system is considered in the

4126

Ind. Eng. Chem. Res., Vol. 42, No. 17, 2003

Figure 1. Temperature T vs density F diagram for a pure fluid. Point a corresponds to the initial one-phase system. During the simulation, the temperature is dropped to a lower value [point (b)], placing the system in an unstable state. The fluid spontaneously separates into two phases in coexistence with each other, a dense liquid and a vapor.

figure, but the method is entirely general. One starts with a one-phase system [point (a)] with a given number of particles N, volume V, and temperature T. The temperature is controlled by means of a “thermostat” affecting the equations of motion. While both N and V are kept constant, the temperature of the system is lowered abruptly by adjusting the thermostat [point (b)]. The target temperature must be such that the resulting state point lies within the spinodal envelope, e.g., at conditions that are both mechanically and thermodynamically unstable. As the system relaxes at the new state point, particles will aggregate to form domains of liquid and of vapor. The connectivity and morphology of these phases will depend on their volume fractions. In a binary system, there is an additional degree of freedom. By specification of the overall composition, temperature, and density, the pressure is given by the mechanical and thermodynamic equilibria. However, most experimental data are given in terms of constantpressure T-x-y data or constant-temperature P-x-y data. Furthermore, ternary or higher systems are presented as phase compositions at given pressure and temperature, and a temperature quench at constant volume is not sufficient to target both a final temperature and pressure. This can be accomplished either by a quench to a position held at constant temperature and pressure or by a temperature quench at constant volume followed by an isothermal-isobaric (NPT) simulation to bring the system to the desired state point. We have found (see Figure 2) that the second recipe appears to be more efficient for dense systems, including liquidliquid and solid-liquid equilibria. If the initial density is chosen such that the vapor phase is present after the quench, the initial constant-volume quench produces a three-phase system, where one phase is a vapor. The separation between the two dense phases is achieved very rapidly as a result of flow. In a quench under NPT conditions, the vapor phase does not appear, and dense phase separation is achieved only through diffusion and takes much longer. 3.2. Simulation Details. We have implemented TQMD both as a serial code that can be run on a workstation and as a spatially decomposed2-4 parallel program suitable for clusters or supercomputers. In the parallel algorithm, each processor is responsible for the storage of data pertaining to particles within its as-

Figure 2. Snapshots of intermediate and final configurations for a binary 50% molar sulfolane (gray spheres) plus n-octane (white spheres) system at 402.35 K. The simulation cell contains 8000 particles. Sphere sizes are scaled with respect to the diameter given by the corresponding values of σ (Table 1). Initially the system is placed a temperature of 2000 K and a density of 7.7 mol/L, corresponding to a well-mixed supercritical one-phase fluid (a). Temperature quench is done isochorically to the desired temperature of 402.35 K via NVT dynamics. Snapshots b-e show the evolution at intervals of 51.8 ps. State f is the final NVT configuration at 2000 ps. At state f, the system is mechanically unstable, registering a negative pressure. Using this last configuration as a starting point, the system is taken to the desired pressure of 101 kPa via NPT dynamics. Snapshots g and h correspond to the evolution at intervals of 518 ps. Snapshot i corresponds to a final configuration after a total of 2000 ps of NVT dynamics followed by 2000 ps of NPT dynamics. The final overall density is 16.8 mol/L. The octane-rich phase appears to be pure; the sulfolane-rich phase has 3% mole fraction n-octane.

signed volume, for the integration of the equations of motion for these particles, for the computation of the forces that these particles exert on each other, and for computation of half of the interactions (forces) between its own particles and those in neighboring domains. Communications between processors are handled using the Message Passing Interface library. This approach is well documented.4 In the calculations performed here, integration of the equations of motion is accomplished with a fifth-order Gear predictor-corrector algorithm.10 The temperature is controlled with the Nose´-Hoover32 thermostat. For NPT simulations, this approach is extended to also couple the volume of the simulation cell to a second external degree of freedom, which drives the system toward an imposed constant pressure.33 After quenching to the desired temperature, the simulation cell presents clusters corresponding to the two coexisting phases. TQMD relies on the fact that a local equilibrium is achieved rapidly, and it is therefore not necessary to run simulations until a single (planar) interface is present in the system. In this heterogeneous scenario, obtaining the coexisting densities and compositions requires suitable analysis. For the generation of probability distributions of densities and mole fractions, the simulation cell is divided into cubic subcells, in each of which the mole fractions of each species and the total

Ind. Eng. Chem. Res., Vol. 42, No. 17, 2003 4127 Table 1. IMP Parameters for the Fluids Used in This Worka substance

/k (K)

σ (Å)

µ (C‚m)

Q (C‚m2)

benzene carbon dioxide cyclohexane 1,2-dichloroethane naphthalene n-octane sulfolane

358.74 215.0 412.0 417.83 530.0 432.22 579.81

5.236 3.748 5.360 4.820 5.5806 4.748 5.230

0 0 0 6.971 × 10-30 0 0 4.57 × 10-29

4.09 × 10-39 1.367 × 10-39 0 0 4.503 × 10-39 0 0

a When available, literature values are used for the multipolar moments. The energy and size parameters are adjusted to saturated liquid densities for multipolar fluids and to critical temperature and density for others.

Figure 3. Temperature T* vs density F* diagram for a fluid with average dipolar attraction of µ* ) x2. Open diamonds are the GEMC results of Sadus;35 solid circles are the results of this work. The uncertainties are comparable with the size of the symbols. The solid line is the phase envelope for the LJ potential (µ* ) 0).

density are calculated and accumulated into histograms. We do not collect histogram data from interfacial regions.1 After the initial temperature quenching, local equilibration is attained rapidly, usually within 10 000 time steps of ∆t* ) ∆t/σxm/ ) 0.005, where m is the mass of the lightest molecule. From this point on, results are averaged from at least five different configurations, each spaced by 1000 time steps. In an effort to improve the statistics without increasing the run length and to ensure that sizable “bulklike” regions of each phase are present, relatively large system sizes are employed. For liquid-vapor coexistence, good results were obtained with N ) 8000 for single-component systems and N ) 32 000 for multicomponent systems. The higher the temperature, the larger the system size that is needed because the liquid-vapor interface becomes broad at high temperature, reducing the volume of the bulk phases and the quality of the data collected. System

sizes of N ) 64 000 were used in extreme cases and for liquid-liquid equilibria, which still complete in a reasonable time if the computation is done with a parallel machine. These calculations have all been performed on a cluster consisting of 20 dual-processor Pentium 550 MHz machines connected via a 2 GB/s Myrinet switch. Unless otherwise noted, we have cut and shifted the IMP potential with a cutoff radius of 5σ (in order to avoid force discontinuities and impulse pressures) and used no long-range corrections; this radius is deemed large enough to allow us to neglect these effects.34 4. Results Figure 3 shows the vapor-liquid equilibrium results for a pure dipolar (Keesom) fluid with µ* ) x2. As a comparison, the solid line corresponds to the underlying LJ fluid (with the same  and σ but µ* ) 0). Upon addition of multipolar attractions, the phase diagram widens considerably and reaches a higher critical temperature. The results obtained using TQMD agree within statistical uncertainties with the GEMC results of Sadus.35 However, TQMD allows the exploration of low-temperature dense states without difficulty. Because IMP is a two-parameter model, one is tempted to use a corresponding states principle to backcalculate the parameters; i.e., if one knows the reduced values of

Figure 4. Temperature vs density F diagram for benzene. Solid circles are simulation results. The solid line is the equation of state of Goodwin.37 The dashed line is the LJ potential (eq 1) with parameters fit to the experimental critical temperature and density (eq 19).

4128

Ind. Eng. Chem. Res., Vol. 42, No. 17, 2003

Figure 5. Saturated density of naphthalene as a function of temperature. Symbols are the TQMD simulations using the IMP potential, squares correspond to liquid densities, and triangles correspond to solid densities. The open symbols correspond to the data points used to fit the potential parameters. Error bars are smaller than the symbols. Dashed and solid lines are the correlated experimental40 liquid and solid data, respectively.

the critical temperature and density, one may calculate the  and σ: 3

 ) kTc/T/c , σ ) xF/c /Fc

(19)

For the LJ fluid, the reduced critical properties are known,36 T/c ≈ 1.31 and F/c ) 0.316, so the calculation is direct. Parameters for nonpolar molecules in Table 1 are found this way. For polar molecules, the procedure is complicated by the fact that critical points are dependent in a nontrivial way to the absolute value of the multipole moment (e.g., Figure 3); therefore, one should have this information beforehand. This not being available, for polar molecules, we have fit the potential parameters to several experimental saturated liquid density points. Figure 4 presents data for pure benzene, a quadrupolar molecule. The energy and size parameters of IMP,  and σ, have been adjusted to the saturated liquid density at 400 K (roughly a reduced temperature of 0.7 with respect to the critical temperature). The experimental25 quadrupole moment is used, and no dipolar moment is considered. The solid line corresponds to a multiparameter equation of state37 which, on the scale of the plot, represents the experimental data exactly. The result is surprisingly good considering the simplicity of the model, i.e., using only two adjustable parameters. We point out that such a quality of fit to experimental data is not possible wth a simple LJ potential because either the high liquid density prediction or the critical point location will then have to be sacrificed. In Figure 5 we show the saturated liquid and solid densities for pure naphthalene. Parameters for naphthalene were fit to some of the high-temperature liquid density data. The results are compared with the correlated experimental data of both the saturated liquid and solid. A failing of the IMP appears at low temperatures where solid phases may appear. Clearly, shortrange repulsion and shape factors will be dominant in determining the structure and density at these condi-

Figure 6. High-pressure temperature vs density diagram for carbon dioxide. Solid circles are simulation results for the phase envelope. Open circles are simulation results for isobars at 2, 10, and 20 MPa. The uncertainties are comparable with the size of the symbols. Lines are the equation of state of Bender.41

tions. In Figure 5, the largest deviations from experimental values are obtained at conditions close to and below the triple point, and the solid-fluid transition is poorly determined. For similar reasons, vapor pressures of solids are not expected to be quantitatively correct. Figure 6 shows a simulated temperature-density diagram for carbon dioxide. For this particular case, because the use of the potential in modeling supercritical fluid extraction is anticipated, the parameters (Table 1) were chosen to faithfully reproduce the critical region.

Ind. Eng. Chem. Res., Vol. 42, No. 17, 2003 4129

Figure 7. Vapor pressures of pure fluids. Squares are simulation results for CO2, triangles are simulation results for 1,2-dichloroethane, and lines are smoothed experimental data.41,42

Figure 8. Vapor composition y1 vs liquid composition x1 for the binary 1,2-dichloroethane (1) + cyclohexane (2) system at 1 atm. Solid symbols are simulations using IMP potentials, the solid line is the smoothed experimental data,42 and the dashed line is the diagonal, placed as a guide to the eye.

Figure 7 shows the vapor pressure of the model along with smoothed experimental results. While for CO2 the prediction is of acceptable quality, it is not so for the dipolar fluid. Vapor pressures of solids, and to a lesser degree of liquids, may be strongly dependent on both the actual shape and energetic anisotropy of the molecules. A better agreement of vapor pressures (to the detriment of the liquid-phase properties) may be obtained if vapor-pressure information is included in the fitting of the parameters. Figure 8 shows a vapor-liquid-phase diagram of a binary mixture involving a polar compound, 1,2-dichloroethane and cyclohexane. The results are plotted in an x-y diagram, where coexisting liquid and vapor compositions are given. Each data point corresponds to a

Figure 9. Equilibrium result for a ternary liquid-liquid equilibrium at 402.35 K and 101 kPa. The top vertex corresponds to pure benzene, the lower left vertex corresponds to pure n-octane, and the lower right vertex corresponds to pure sulfolane. Solid circles are simulation results using IMP results, open squares are experimental results of Lee and Kim,38 and the dashed line is a UNIQUAC correlation with data from ternary experimental results.

different temperature. The system is azeotropic; i.e., the data cross the diagonal, which is typical of a polarnonpolar mixture. The diagram is a prediction because the model parameters are fit only to pure-component data. The pressure-composition and pressure-temperature diagrams obtained are in only qualitative agreement with the experiment (because of the poor fit to vapor pressures) and are not shown here. In Figure 9, we show the liquid-liquid equilibrium of a ternary system in which benzene, a quadrupolar fluid, is mixed with n-octane, a nonpolar component, and sulfolane, an oxygenated thiophene that is a strongly non-hydrogen-bonding polar solvent. The large differences in polarity produce liquid-liquid-phase separation between the sulfolane and the octane (cf. Figure 2). Benzene, which is soluble in both liquids, is separated among the two phases. The results are plotted alongside the experimental results of Lee and Kim,38 and the

4130

Ind. Eng. Chem. Res., Vol. 42, No. 17, 2003

correlation is made by the same authors using the UNIQUAC model. Despite the simplicity of our model, the results are in quantitative agreement. Some details are not well modeled by IMP. Figure 2, corresponding to the binary sulfolane/n-octane system, shows a solidlike structure in the sulfolane phase, probably caused by an incorrectly large attraction. This can be attributed, in part, to a poor fit of the sulfolane parameters, for which very little experimental data are available in the open literature. If one were to use the dipole moment of sulfolane as an adjustable parameter, a better representation might be possible of both the pure fluid and its behavior in the mixture. It is, nevertheless, striking that a simple model like IMP has the ability to model such complex phase behavior. 5. Conclusions We have presented a simple potential model where the shape of the molecule is assumed spherical and the polar directional attractions have been averaged to yield an isotropic functional form. Despite its simplicity, the potential is capable of correctly reproducing the phase behavior, including liquid-liquid separation and liquidvapor equilibria, of a number of nontrivial test systems. By separation of the permanent polar terms from the nonpolar (dispersion) attractions, the parameters ought to be well behaved. This is expected to be of importance when dealing with mixtures where binary interaction parameters or similar artifices are often invoked in order to fit data. The microscopic view of fluid-phase equilibria afforded by molecular simulations provides insight as to why these simplified approaches are so successful in some cases. The IMP potential resembles the “cubic equation of state” approach common in engineering, i.e., the family of the so-called van der Waals equations of state. Today, we fully understand the shortcomings of equations of this latter type: both the repulsive term and the dispersion term are too simplified. However, these equations are the basis of most modern fluid-phase calculations, even though physically sensible alternatives are available, e.g., SAFT.39 Cubic equations of state produce acceptable results with minimal effort. In this sense, the results shown in this paper follow an analogous direction. They illustrate that an important contribution to fluid-phase equilibria is the average energetics. If a full description of the volumetric properties of a fluid, e.g., a precise simultaneous fit of PvT data and the derived heat capacities, free energies, etc., is sought, a detailed structured potential seems necessary; in the same way, a multiparameter equation of state is the current engineering solution to the problem of correlating high-precision data for fluids. As inspection of Figure 7 suggests, vapor-pressure calculations can be in error if the parameters are fit only to liquid densities. Particularly in this case, if no multipolar contribution is applied, IMP reduces to a LJ model that is incapable of simultaneously fitting both vapor pressures and liquid densities of nonspherical molecules. This is mainly due to the sensitivity of the vapor-pressure calculations to molecular details. While this is of lesser importance in some cases such as liquid-liquid equilibrium, it is crucial for vapor-liquid equilibrium. However, the multipolar contribution (eqs 12-14) may be added to a more complex potential, e.g., one which incorporates a chain term or a suitable shape factor. Likewise, shape details, short-range directional

interactions, and association are likely to be important for solidification, liquid-crystalline behavior, interfaces, adsorption, and transport properties. IMP is not expected to be suitable for simulating such phenomena. It is likely that this is also the reason most generalized equations of state often fail when modeling these phenomena. Last, we point out that the TQMD approach is a promising and rapid method for obtaining phase behavior by molecular simulation. The ever-increasing performance of inexpensive computers and the popularity of clustered parallel systems suggest that we are not far from the time that fluid-phase simulations will rival equations of state for engineering purposes. Acknowledgment E.A.M. acknowledges the financial support given by the Agenda Petroleo (Project 97-003530) to this work. Some calculations were performed by Marı´a E. Sua´rez, Susana Curbelo, and Simo´n Albo at the Laboratorio de Computacio´n de Alto Rendimiento, Universidad Simo´n Bolı´var. Literature Cited (1) Gelb, L. G.; Mu¨ller, E. A. Location of phase equilibria by temperature-quench molecular dynamics simulations. Fluid Phase Equilib. 2002, 203, 1. (2) Fincham, D. Parallel Computers and Molecular Simulation. Mol. Simul. 1987, 1, 1. (3) Plimpton, S. Fast parallel algorithms for short-range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1. (4) Heffelfinger, G. S. Parallel atomistic simulations. Comput. Phys. Commun. 2000, 128, 219. (5) http://www.amber.ucsf.edu/amber/amber.html. (6) http://yuri.harvard.edu/. (7) http://www.dl.ac.uk/TCSC/Software/DL_POLY/main.html. (8) http://igc.ethz.ch/gromos/. (9) http://www.earth.ox.ac.uk/∼keith/moldy.html. (10) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, U.K., 1987. (11) Leach, A. R. Molecular Modelling; Longman: Essex, U.K., 1997. (12) Sadus, R. J. Molecular Simulation of Fluids; Elsevier: Amsterdam, The Netherlands, 1999. (13) Price, S. L. Toward more accurate model intermolecular potentials for organic molecules. Rev. Comput. Chem. 2000, 14, 225. (14) Panagiotopoulos, A. Z. Force field development for simulations of condensed phases. AIChE Symp. Ser. 2001, 97, 61. (15) Poling, B. E.; Prausnitz, J. M.; O’Connell, J. P. The properties of gases and liquids, 5th ed.; McGraw-Hill: New York, 2000. (16) MacKerell, A. D.; Wiorkiewiczkuczera, J.; Karplus, M. An all-atom empirical energy function for the simulation of nucleic acids. J. Am. Chem. Soc. 1995, 117, 11946. (17) Cornell, W. D.; Cieplak, P.; Bayly, C. L.; Gould, I. R.; Merz, K. M., Jr.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A second generation force field for the simulation of proteins, nucleic acids and organic molecules. J. Am. Chem. Soc. 1995, 117, 5179. (18) Halgren, T. A. Merck Molecular Force Field. I. Basis, form, scope, parametrization, and performance of MMFF94. J. Comput. Chem. 1996, 17, 491. (19) Sun, H. COMPASS: An ab initio force-field optimized for condensed-phase applications. Overview with details on alkane and benzene compounds. J. Phys. Chem. B 1998, 102, 7338. (20) Engkvist, O.; Astrand, P. O.; Karlstrom, G. Accurate intermolecular potentials obtained from molecular wave functions: Bridging the gap between quantum chemistry and molecular simulations. Chem. Rev. 2000, 100, 4087.

Ind. Eng. Chem. Res., Vol. 42, No. 17, 2003 4131 (21) Rowley, R. L.; Yang, Y.; Pakkanen, T. A. Determination of an ethane intermolecular potential model for use in molecular simulations from ab initio calculations. J. Chem. Phys. 2001, 114, 6058. (22) Martin, M. G.; Siepmann, J. I. Transferable potentials for phase equilibria. 1. United-atom description of n-alkanes. J. Phys. Chem. B 1998, 102, 2569. (23) Toxvaerd, S. Equation of state of alkanes. 2. J. Chem. Phys. 1997, 107, 5197. (24) Errington, J. R.; Panagiotopoulos, A. Z. New intermolecular potential models for benzene and cyclohexane. J. Chem. Phys. 1999, 111, 9731. (25) Gray, C. G.; Gubbins, K. E. Theory of Molecular Fluidss Volume 1: Fundamentals; Clarendon Press: Oxford, U.K., 1984. (26) Luckas, K. Applied Statistical Thermodynamics; SpringerVerlag: Berlin, 1992. (27) Buckingham, A. D. Molecular quadrupole moments. Q. Rev. 1959, 13, 183. (28) Reed, T. M.; Gubbins, K. E. Applied Statistical Mechanics; McGraw-Hill: New York, 1973. (29) Hirschfelder, J. O.; Curtiss, C. F.; Bird, R. B. Molecular Theory of Gases and Liquids; Wiley: New York, 1954. (30) Machado, J. M. V.; Zabaloy, M. S.; Macedo, E. A. Saturated vapor pressure through a modified Lennard-Jones equation of state. Fluid Phase Equilib. 2001, 182, 75. (31) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926. (32) Nose´, S. A molecular dynamics method for simulations in the canonical ensemble. Mol. Phys. 1984, 52, 255. (33) Hoover, W. G. Canonical dynamics: equilibrium phasespace distribution. Phys. Rev. A 1985, 31, 1695. (34) Mecke, M.; Winkelmann, J.; Fischer, M. Molecular dynamics simulation of the liquid-vapor interface: the Lennard-Jones fluid. J. Chem. Phys. 1997, 107, 9264.

(35) Sadus, R. J. Molecular simulation of the vapour-liquid equilibria of pure fluids and binary mixtures containing dipolar components: the effect of Keesom interactions. Mol. Phys. 1996, 87, 979. (36) Shi, W.; Johnson, J. K. Histogram reweighing and finitesize scaling study of the Lennard-Jones fluids. Fluid Phase Equilib. 2001, 187, 171. (37) Goodwin, R. D. Benzene thermophysical properties from 279 to 900 K at pressures to 1000 bar. J. Phys. Chem. Ref. Data 1988, 17, 1541. (38) Lee, S.; Kim, H. Liquid-liquid equilibria for the ternary systems sulfolane + octane + benzene, sulfolane + octane + toluene, and sulfolane + octane + p-xylene at elevated temperatures. J. Chem. Eng. Data 1998, 43, 358. (39) Mu¨ller, E. A.; Gubbins, K. E. Molecular-based equations of state for associating fluids: A review of SAFT and related approaches. Ind. Eng. Chem. Res. 2001, 40, 2193. (40) Daubert, T. E.; Danner, R. P. Physical and Thermodynamic Properties of Pure Chemicals; National Standard References Data System (NSRDS); AIChE: New York, 1992; Vols. 1 and 3. (41) Bender, E. Equations of state exactly representing the phase behavior of pure substances. Proc. Symp. Thermophys. Prop. 1970, 5, 227. (42) Stephan, K.; Hildwein, H. Recommended data of selected compounds and binary mixtures; DECHEMA Chemistry Data Series; DECHEMA: Frankfurt, Germany, 1987; Vol. IV, Parts 1 and 2.

Received for review January 15, 2003 Revised manuscript received May 29, 2003 Accepted June 11, 2003 IE030033Y