Metastable and Unstable Fluid States by Molecular Simulation

Nov 7, 1996 - ... is to sample by the Metropolis Monte Carlo method [Metropolis et al., ..... give an indication of the agreement between laboratory e...
0 downloads 0 Views 226KB Size
Download PDF
4336

Ind. Eng. Chem. Res. 1996, 35, 4336-4341

Metastable and Unstable Fluid States by Molecular Simulation B. S. Watson, R. A. Greenkorn, and K. C. Chao* School of Chemical Engineering, Purdue University, West Lafayette, Indiana 47907-1283

A method of determination of thermodynamic properties of metastable and unstable fluid states by molecular simulation is presented. As illustrative examples spinodal and binodal states are calculated for the Lennard-Jones fluid, ethane, and n-butane using the optimized potentials for liquid simulations for the representation of the real molecules. Introduction Metastable, unstable, and spinodal states are homogeneous fluid states at intermediate densities between the saturated liquid and the saturated vapor. These states are the poorest understood of fluid states. Experimental investigation of the metastable states is difficult though, in principle, possible with expanded effort. Spinodal states have rarely been observed. The unstable states do not exist physically and cannot be observed. These states are, nevertheless, important to the understanding of fluids. The metastable and spinodal states are the states of interest in investigations of nucleation kinetics. Unstable and metastable states are routinely used in integration for enthalpy, free energy, entropy, etc., of liquids in thermodynamic calculations by equation of state. Equation of state is used to calculate metastable and spinodal states in nucleation kinetics studies, but the equations of state are of unknown reliability for this purpose since they are fitted to experimental data in the stable states and extrapolated to the metastable and spinodal states. Efforts have been made to calculate the metastable and unstable states by molecular simulation. The excessive fluctuation caused by molecular clustering has made the simulation difficult or impossible. Results have been obtained only for a few model molecules with the use of specialized calculational devices. Wood (1968) examined the metastable and unstable states of the Lennard-Jones fluid at the isotherm kT/ ) 1.079 by using a cell of 32 molecules. It is not possible for the separate phases to be formed with such a small number of molecules. Although the fluctuations remained manageable, the uncertainties are large. Hansen and Verlet (1969) repressed large fluctuations in the metastable and unstable states by subdividing the simulation cell into boxes and setting upper and lower limits for the number of molecules allowed in a box. However, this method of restraining density fluctuations proved ineffective near the critical point. Torrie and Valleau (1974) carried out perturbation calculations to obtain the free energy of Lennard-Jones fluids at subcritical states by using the purely repulsive molecules u ) 4(σ/r)12 as the base fluid. The 12th power repulsive molecules, like hard spheres, do not tend to form heterogeneous phases. The method requires the free energy of the base fluid to be known. Poole et al. (1993) calculated spinodal states of subcooled liquid water at subfreezing temperature by using molecular dynamic calculations at short times of nanoseconds or shorter. In this work, the residual Helmholtz free energy at metastable and unstable states is found with the * Author to whom correspondence is addressed. Phone: (317) 463-3270. Fax: (317) 494-0805. email: [email protected].

S0888-5885(96)00126-1 CCC: $12.00

temperature perturbation method (Watson and Chao, 1992) of Monte Carlo molecular simulation. The residual Helmholtz energy at (N, V, T) is determined by summing the logarithms of the ensemble averages of the Boltzmann exponential function at 2T, 4T, 8T, ...ad infinitum. By not performing any molecular simulation at the T of interest but starting at 2T and proceeding to higher temperatures, molecular clustering and excessive fluctuations are avoided, thus allowing rigorous molecular calculations of metastable and unstable fluid states. The spinodal and binodal states are determined from the spectrum of isothermal states thus obtained. Results are reported for the Lennard-Jones fluid, ethane, and n-butane, in this paper. The optimized potentials for liquid simulations (OPLS) of Jorgensen et al. (1983, 1984) are used in this work to represent the molecules of the real fluids. These are atom-to-atom potentials interacting by the Lennard-Jones function plus partial electric charges interacting by Coulombic law. We refer the reader to Jorgensen et al. for parameter values and structure details of the real molecules. Since no experimental data are available or expected soon for comparison with our calculated unstable and metastable states, we make a comparison of the calculated saturated states with experimental data to give an indication of the agreement of the simulated with the real fluid property. We report in the following comparisons for vapor pressure, saturated liquid volume, and saturated vapor volume for the fluids studied in this work. The method of calculation being rigorous, the agreement, or lack of it, between the calculation and experiment depends on the molecular model of the simulation and secondarily on the numerical accuracy of the calculation. The OPLS potentials that are used in this work will be seen to give simulation results that vary to some extent from one fluid to another but are generally good descriptions of the molecules. Temperature Perturbation Method The first step in the temperature perturbation method for the determination of residual Helmholtz energy Ar in the canonical ensemble is to sample by the Metropolis Monte Carlo method [Metropolis et al., 1953] the Boltzmann exponential factor at twice the temperature of interest and to obtain the ensemble average denoted by 〈 〉.

〈exp[-U/k2T]〉2T )

∫exp(-U/2kT) exp(-U/2kT) drN ) ∫exp(-U/2kT) drN ∫exp(-U/kT) drN ∫exp(-U/2kT) drN

© 1996 American Chemical Society

(1)

Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996 4337

The integrals in the numerator and denominator on the right-hand side of eq 1 are the configurational integrals at temperatures of T and 2T, respectively. Taking the logarithm of eq 1 gives

( ) ( )

〈

( )〉

ArT Ar2T -U ) -ln exp kT k2T k2T

(2)

2T

Repeating this procedure at temperatures of 4T, 8T, ... and summing the equations leads to

()( ) ArT

Ar2mT

-

kT

〈 [ ]〉

m

ln exp ∑ i)1

)-

m

k2 T

-U

k2iT

(3) 2iT

At a sufficiently large value of m, i.e., at a sufficiently high temperature, the second term on the left-hand side of eq 3 approaches either zero for molecules without a hard core or a constant value for molecules with a hard core. Theoretical results for the residual Helmholtz energy of the hard core are available for a number of hard cores. By successively sampling and averaging exp(-U/kT) for a N-molecule cell at 2T, 4T, ..., summing the series, and combining with the high-temperature Ar value, if any, according to eq 3, the Ar at T is obtained. The Boltzmann factors in eq 3 are difficult to compute because U, the total configurational energy of N molecules, is large compared to kT. This causes the Boltzmann factor to be either nearly zero or excessively large. This problem can be avoided if the Boltzmann factor is based on the energy of a single molecule rather than N molecules. The potential energy of N molecules is expressed as the sum of the energies of the single molecules:

ui ∑ i)1

(4)

Substitute eq 4 in the Boltzmann factor

〈 ( )〉 〈 ( )〉 N

-U

exp

)

kT

exp ∏ i)1

-ui kT

(5)

The average of the products of uncorrelated statistical quantities is equal to the product of the averages of the quantities. If the exp(ui/kT) in eq 5 are uncorrelated, then the right-hand side of eq 5 is:

〈 ( )〉 ∏〈 ( )〉 N

-U

exp

exp

)

kT

i)1

()( ) arT

ar2mT

-

kT

〈 ( )〉

∞

ln exp ∑ i)1

)-

k2mT

-u

k2iT

(7)

2iT

Equation 7 expresses single-particle sampling as the method of implementation of the temperature perturbation method of eq 3. It is the method of this work. However, it is not a general method. Thus, the neighbor particles in a solid are not transient like in a fluid but are fixed and always share interaction energy. The energies of neighbors are correlated. For the onedimensional Ising model used by Talbot et al. (1993) to show the approximate nature of single-particle sampling, the energies of neighbors have been found to have a correlation coefficient of 0.5 (Bereolos et al., 1997). The factorization of eq 6 does not apply, and eq 7 does not hold in this case except in an approximate sense. In new applications of single-particle sampling, the absence of correlation of the particle energies should be established. Lennard-Jones (L-J) Fluid

N

U)

That the correlation between the exponential factors of eq 6 is negligibly small in fluid states is confirmed in the course of the molecular simulation of the LennardJones molecules. Correlation between the energies of two molecules is evaluated for all pairs in the simulation cell. Sampling pairs confined in a cell should greatly exaggerate the correlation from that of pairs in a bulk fluid. Nevertheless, no significant nonzero correlation is found to exist between any two molecules (Bereolos et al., 1996). Since the molecules are indistinguishable, the subscript i in eq 6 can be dropped and the continuous product replaced by the Nth power of exp[-u/kT], eq 3, expressed in single-particle form is

-ui kT

(6)

The L-J fluid has been widely used in example calculations of theoretical methods. While the potential function is simple, the L-J fluid gives excellent representation for the noble gases argon, krypton, and xenon and, to a lesser extent, also for methane. To promote convergence of eq 7 at high temperatures for L-J molecules, a hard sphere of d ) 0.7σ was inserted in the molecular center where σ is the zero potential distance of separation of two molecules. The potential function is unchanged at all molecular distances r > 0.7σ. The change in potential at r < 0.7σ is at such a high value of uij that the probability of its occurrence at moderate temperatures is negligible. For a molecule with a hard core, eq 7 becomes

arT

Are the factors exp(-ui/kT) in eq 6 for different molecules uncorrelated? Any two molecules are generally not together but separated at a distance. There is no correlation for the fluctuation of energies of molecules at a distance in a fluid. The radial distribution function and pair correlation function approach the asymptotic no correlation behavior at microscopic distances between a pair of molecules. Although correlation exists in the microscopically close range, being at close range is a transient state for any two molecules. The probability of two molecules coming together to a close range is insignificant. The energies of molecules in a fluid are not correlated as they are moved in the simulation process. Equation 6 is true for fluids.

kT

[] ar

)

kT

∞

hc

-

〈 [ ]〉

ln exp ∑ i)1

-u

k2iT

(8)

2iT

Figure 1 shows the radial distribution function at T* ) 1.0 of a L-J fluid and the same function of the hardcore augmented L-J fluid. The two functions coincide, which shows the hard-core augmented L-J fluid behaves the same as the L-J fluid at ordinary temperatures. It takes about 12 terms in eq 8 to achieve convergence for the augmented molecules; many more terms are required for the unaugmented molecules. The Carnahan-Starling (1969) equation of state is used to obtain the hard-core Helmholtz energy in eq 8. Free-energy isotherms were calculated for the L-J fluid at 12 temperatures from kT/ ) 0.8 up to 1.32. At

4338 Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996 Table 1. Residual Helmholtz Energy of the Lennard-Jones Fluid F*

Figure 1. Radial distribution function of the Lennard-Jones fluid at T* ) 1.0 and F* ) 0.7: (4) with no hard core; (b) with a hard core of d* ) 0.7.

Figure 2. Residual Helmholtz energy of the Lennard-Jones fluid in the metastable and unstable states.

each temperature, the residual Helmholtz energy is obtained using eq 8 at four to eight densities. At each of these densities, Monte Carlo simulations in the canonical (NVT) ensemble were run to evaluate the first 16 terms of eq 8. For each term simulations were performed with 216 molecules, applying a periodic boundary condition at the cell edges to negate surface effects. Each simulation was run for 3000 cycles, where each cycle consisted of an attempt to move each molecule. To avoid bias caused by the choice of the initial configuration, results from the first 50 cycles were not included in the computation of the averages used in eq 8. For each simulation, the potential cutoff limit was set at half the box length. To account for molecular interactions at distances greater than the cutoff limit, the molecular density beyond the cutoff distance is assumed to be uniform. The interaction energy between a molecule and the molecules outside the cutoff distance is computed by integrating the product of the density

ar/kT

F*

ar/kT

0.300 0.500

T* ) 0.800 -1.689 0.700 -2.833 0.800

-3.674 -3.910

0.200 0.400

T* ) 0.928 -1.009 0.600 -1.862 0.700

-2.560 -2.791

0.200 0.400

T* ) 1.000 -0.879 0.600 -1.627 0.700

-2.225 -2.403

0.200 0.500

T* ) 1.058 -0.795 0.600 -1.757 0.700

-1.991 -2.141

0.200 0.400

T* ) 1.100 -0.742 0.600 -1.371 0.800

-1.841 -1.944

0.400 0.500

T* ) 1.150 -1.262 0.600 -1.496 0.700

-1.679 -1.778

0.200 0.400 0.600 0.650

T* ) 1.200 -0.629 0.725 -1.164 0.750 -1.527 0.800 -1.569

-1.600 -1.587 -1.532

0.200 0.300 0.400 0.500

T* ) 1.250 -0.583 0.600 -0.842 0.700 -1.067 0.800 -1.253

-1.388 -1.440 -1.364

0.200 0.300 0.400

T* ) 1.280 -0.556 0.450 -0.801 0.500 -1.016 0.600

-1.109 -1.186 -1.311

0.200 0.300 0.400

T* ) 1.300 -0.541 0.450 -0.773 0.500 -0.981 0.600

-1.071 -1.150 -1.260

0.100 0.200 0.300

T* ) 1.310 -0.277 0.400 -0.532 0.450 -0.765 0.500

-0.967 -1.051 -1.127

0.200 0.300

T* ) 1.320 -0.524 0.400 -0.753 0.500

-0.949 -1.106

and the potential function from the cutoff distance to infinity. A factor of 1/2 is then applied to the integral to prevent double counting. Table 1 presents the results. The simulated isothermal point values are fitted with a cubic polynomial in F that passes through the origin. For some of the high-temperature LennardJones isotherms, the residual free energy was calculated at a relatively large number of densities to verify that the cubic polynomial gives a good representation of the residual free energy. Establishing the adequacy of the cubic polynomial with the Lennard-Jones fluid permitted the free energy to be calculated at a smaller number of densities for the ethane and n-butane isotherms. The simulated data and the fitted polynomials for several temperatures are shown in Figure 2 for the L-J molecules. The total Helmholtz energy was obtained by adding the ideal gas free energy to the residual value. One isotherm, at a reduced temperature of kT/ ) 1.00 obtained from the fitting polynomial, is shown in Figure 3 as a function of V. A strong maximum is observed in the unstable region where the second derivative with respect to V is negative; the temperature is below the critical.

Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996 4339 Table 2. Binodal and Spinodal States of the Lennard-Jones Fluida binodal state

spinodal state

V/(Nσ3)

a

liq side

gas side

kT/

liq

gas

Pσ3/

V/(Nσ3)

Pσ3/

V/(Nσ3)

Pσ3/

0.800 0.928 1.000 1.058 1.100 1.150 1.200 1.250 1.280 1.300

1.22 1.21 1.37 1.48 1.68 1.72 1.93 2.14 2.33 2.59

159 53.5 31.2 21.5 16.9 13.5 10.6 7.93 5.77 4.95

0.0049 (0.0047) 0.0157 (0.0148) 0.0275 (0.0248) 0.0399 (0.0358) 0.0508 (0.0457) 0.0642 (0.0598) 0.0801 (0.0769) 0.101 (0.0972) 0.120 (0.1111) 0.129 (0.121)

1.64 1.65 1.82 1.94 2.16 2.18 2.40 2.58 2.71 2.90

-7.56 -0.378 -0.238 -0.143 -0.080 -0.038 0.0201 0.0729 0.112 0.1265

11.4 9.77 8.46 7.60 7.16 6.47 6.02 5.33 4.52 4.12

0.0367 0.0462 0.0575 0.0673 0.0753 0.0868 0.0971 0.111 0.123 0.1312

Values in parentheses are from Lofti et al. (1992).

Figure 3. Stable, metastable, and unstable regions of the Lennard-Jones fluid at T* ) 1.0.

The stable, metastable, and unstable states of a fluid at a particular temperature are determined by the Helmholtz energy as a function of volume. A negative second derivative (∂2A/∂V2)T indicates an unstable state. A homogeneous fluid in an unstable state splits into two phases. The two points at which (∂2A/∂V2)T is zero are the spinodal states. These represent the theoretical limit of superheating of a liquid or subcooling of a vapor. The saturation of binodal states of the liquid and vapor phases are found as the two points on the A(V) curve that share a common tangent. The vapor pressure is obtained from the slope of the common tangent. Fluids with volumes greater than the saturated vapor volume or volumes less that the saturated liquid volume are in stable states. Homogeneous fluids with volumes between the binodal and spinodal are metastable. The spinodal and binodal states for a Lennard-Jones fluid obtained from the free-energy isotherms are given in Table 2. Figure 4 shows as points the spinodal and saturation volumes obtained in this work. The saturation curves determined by Lofti and co-workers (1992) are also shown for comparison. Ethane Ethane was studied at four temperatures, 150, 200, 250, and 280 K, based on the OPLS potential function for ethane of Jorgensen et al. (1984). Rotation of the

Figure 4. Saturated and spinodal volumes of the Lennard-Jones fluid. The saturation envelope curve is drawn from data in Lofti et al. (1992): (4) saturated volume; this work; (9) spinodal volume; this work.

ethane molecule was achieved by generating a random unit vector. The angle between this vector and the original orientation of the molecule was measured. If this angle was greater than a threshold value, a new random unit vector was generated and the process was repeated until a new orientation with an angle smaller than the threshold value was encountered. An attempt was made to accept the new orientation by calculating the energy change caused by the rotation and using the Metropolis criteria to determine if the rotation would actually be accepted. The value of the threshold angle is periodically adjustable to keep the rate of acceptance close to 50%. A hard dumbbell with lobes of diameter 0.7σ was inserted into the molecules to facilitate convergence. Tildesley and Streett’s (1980) equation of state for hard dumbbells was used to obtain (ar/kT)∞. For each temperature, the residual free energy was calculated at three or four densities. As with the Lennard-Jones fluid, a third-order polynomial was used to fit the data. Sixteen terms in the series in eq 8 were computed to achieve convergence. The NVT simulation for each term used a simulation cell of 180 molecules carried out for 1500 cycles. The results of the first 50 cycles were ignored to prevent the initial configuration from biasing

4340 Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996 Table 6. Residual Helmholtz Energy of n-Butane ar/kT

F/(mol/L)

Figure 5. Saturated and spinodal volumes of ethane. The saturation envelope curve is drawn from data in Perry and Chilton (1973): (4) saturated volume; this work; (9) spinodal volume; this work. Table 3. Residual Helmholtz Energy of Ethane ar/kT

F/(mol/L)

ar/kT

F/(mol/L)

6.176 12.35

T ) 150 K -3.028 18.53 -5.159 21.56

6.176 12.35

T ) 200 K -1.858 18.53 -3.280

-6.805 -6.791

2.790 5.579

T ) 250 K -1.7688 9.848 -3.467

-5.154

2.790 5.579

T ) 300 K -1.368 9.848 -2.5191

-3.569

2.790 4.479

T ) 350 K -1.0213 9.848 -1.8398

-2.4797

2.790 5.579

T ) 400 K -0.7754 9.848 -1.402

-1.6538

Table 7. Saturated Volumes and Vapor Pressures of n-Butanea T/K

Vl/(L/mol)

Vv/(L/mol)

P/bar

250 300 350 400

0.09445 0.09893 0.1171 0.1568

47.95 8.610 2.383 0.9906

0.428 (0.392) 2.72 (2.581) 10.2 (9.46) 24.0 (25.0)

a

Values in parentheses are from Perry and Chilton (1973).

Table 8. Spinodal States of n-Butane -4.070

T/K

V1/(L/mol)

P1/bar

V2/(L/mol)

P2/bar

6.176 12.35

T ) 250 K -1.252 18.53 -2.184

-2.444

0.1287 0.1357 0.1564 0.1969

-343 -181 -66.9 0.520

1.274 0.9813 0.7058 0.5308

8.26 12.5 19.8 30.0

6.176 12.35

T ) 280 K -1.006 18.53 -1.740

250 300 350 400

-1.781

Table 4. Saturated Volumes and Vapor Pressures of Ethanea T/K

Vl × 102/(L/mol)

Vv/(L/mol)

P/bar

150 200 250 280

4.902 5.214 6.938 8.905

119.3 7.893 1.435 0.7338

0.104 (0.0967) 2.026 (2.174) 12.39 (13.01) 24.61 (28.06)

a

ar/kT

F/(mol/L)

Values in parentheses are from Perry and Chilton (1973).

Table 5. Spinodal States of Ethane T/K

V1(L/mol)

P1/atm

V2/(L/mol)

P2/atm

150 200 250 280

0.06854 0.07237 0.09177 0.1034

-447 -273 -93.5 -27.8

0.9796 0.6163 0.4075 0.3267

6.26 13.1 24.8 34.7

the final average. As was done with the Lennard-Jones simulations, the cutoff distance was chosen to be half the cell length, and the interactions beyond the cutoff length were accounted for by assuming uniform density beyond the cutoff distance. The computed residual free energies are listed in Table 3. The binodal and spinodal volumes are presented in Figure 5, where the saturated volumes from the literature (Perry and Chilton, 1973) are also presented as curves for comparison. The volume and pressure of the binodal states at the simulated temperatures are given in Table 4, where the pressure is compared with the vapor pressure from the literature (Perry and Chilton, 1973) to give an indication of the agreement between laboratory experiment and molecular simulation of this work. Table 5 lists the calculated spinodal states of ethane. Since experimental data are not available for the spinodal states, comparison for the binodal states is the closest indication that can be obtained.

n-Butane n-Butane was simulated at 250, 300, 350, and 400 K based on the OPLS potential of Jorgensen et al. (1984). The CH2-CH2 bond of butane was rotated similarly to the CH3-CH3 bond of the ethane molecule. In addition, the methyl groups were rotated around the central bond while maintaining a constant CH3-CH2-CH2 angle. The molecule was not explicitly required to remain in a stable trans or one of the two metastable cis configurations, but the high energy of the other unstable configuration ensured that the molecules spent most of their time in a cis or trans configuration. A hard sphere of diameter 0.7σ was placed on each carbon atom to facilitate convergence of eq 8 which was attained with the use of 16 terms. The canonical average in each of these terms was computed from a canonical (NVT) simulation using 180 molecules in a cell with over 1000 cycles, with the results from the initial 50 cycles not included in the final average. Half the cell length was used as the interaction cutoff distance, and interactions beyond the cutoff distance were integrated by assuming a uniform density beyond the cutoff distance. Boublik et al.’s (1990) hard-chain equation of state was used to obtain the free energy of the hard-core fluid. Table 6 presents the residual Helmholtz energy of n-butane at the four temperatures studied. Table 7 presents the saturated volumes and vapor pressures derived from the simulated Helmholtz energy. Vapor pressure data from the literature (Perry and Chilton, 1973) are included in Table 7 for comparison with the simulated values. Table 8 presents the spinodal states. The volumes of the saturated vapor, the saturated liquid, and the spinodal states of this work are shown as points in Figure 6, where the saturated vapor and liquid volumes from the literature are shown as curves for comparison. Generally good agreement is obtained

Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996 4341

advantage of not performing any calculations at the temperature of interest. It can be expected that the same advantage will prevail when the method is applied to the critical state and near critical states where clustering and fluctuation can be overwhelming to defeat simulations at the temperature of interest. Literature Cited

Figure 6. Saturated and spinodal volumes of n-butane. The saturation envelope curve is drawn from data in Perry and Chilton (1973): (4) saturated volume; this work; (9) spinodal volume; this work.

between the simulated and the experimental values except for the saturated vapor volume at 400 K. Conclusion and Discussion Helmholtz free energies of metastable and unstable fluid states are obtained by molecular simulation with the temperature perturbation method for the LennardJones fluid and for two real fluids for which a good molecular potential function is available. Pressurevolume-temperature (PVT) properties are derived from the simulated free energies. Comparison of the derived PVT properties with literature values confirms the validity of the method applied to metastable and unstable states. The method being rigorous, results of increasing accuracy can be expected as improved molecular potential functions are developed and as numerical methods for the derivation of PVT properties are improved. The simulation method of this work is unique in that the Monte Carlo calculations start at a temperature equal to two times the temperature of interest and subsequent calculations are performed at ever higher temperatures. No calculations are required at the temperature of interest. The success of the method at unstable and metastable states demonstrates the special

Bereolos, P.; Talbot, J.; Chao, K. C. Simulation of Free Energy Without Particle Insertion in the NPT Ensemble. Mol. Phys., submitted for publication. Boublik, T.; Vega, C.; Diaz-Pena, M. J. Chem. Phys. 1990, 93, 130. Carnahan, N. F.; Starling, K. E. Equation of State for Nonattracting Rigid-Spheres. J. Chem. Phys. 1969, 51, 635-6. Hansen, J.; Verlet, L. Phase Transitions of the Lennard-Jones System. Phys. Rev. 1969, 184, 151. 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. Jorgensen, W. L.; Madur, J. D.; Swenson, C. J. Optimized Intermolecular Potential Functions for Liquid Hydrocarbons. J. Am. Chem. Soc. 1984, 106, 6638. Lofti, A.; Vrabec, J.; Fischer, J. Vapor Liquid Equilibria of the Lennard-Jones Fluid from the NpT plus Test Particle Method. Mol. Phys. 1992, 76, 1319. Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. Equation of State Calculations by Fast Competing Machines. J. Chem. Phys. 1953, 21, 1087. Perry, R. H., Chilton, C. H., Eds. Chemical Engineersapo Handbook; McGraw-Hill Book Co.: New York, 1973. Poole, P. H.; Sciortino, F.; Essmann, U.; Stanley, H. E. Spinodal States of Water. Phys. Rev. E 1993, 48, 3799. Talbot, J.; Bereolos, P.; Chao, K. C. Estimation of Free Energy via Single Particle Sampling in Monte Carlo Simulations. J. Chem. Phys. 1993, 98, 1531. Tildesley, D. J.; Streett, W. B. Mol. Phys. 1980, 41, 85. Torrie, G. T.; Valleau, J. P. Monte Carlo Free Energy Estimates using Non-Boltzmann Sampling: Application to the Sub-critical Lennard-Jones Fluid. Chem. Phys. Lett. 1974, 28, 578. Watson, B. S.; Chao, K. C. A Method of Molecular Simulation of Free Energy. J. Chem. Phys. 1992, 96, 9046. Wood, W. W. Monte Carlo Studies of Simple Liquid Models. In Physics of Simple Liquids; Temperley, H. N. V., Rowlinson, J. S., Rushbrook, G. S., Eds.; North Holland: Amsterdam, The Netherlands, 1968; pp 205-214.

Received for review March 1, 1996 Revised manuscript received September 3, 1996 Accepted September 6, 1996X IE960126Y Abstract published in Advance ACS Abstracts, October 15, 1996. X