(455e) Studying Thermophysical Properties with Molecular Dynamics

Vapor Liquid Equilibrium has been
studied by molecular simulation since the mid 1970's [1,2] with Monte Carlo methods being the most prevalent [3]. Monte Carlo methods have some difficulties when studying difficult systems [4] such as dense phases, so techniques have been developed to study VLE with molecular dynamics. In order to simulate the two phase envelope with molecular dynamics within the same simulation cell, two new techniques have been developed. One involves a temperature quench from above the critical point to a temperature that corresponds to a temperature inside the two phase envelope. The second uses a volume expansion to force the system density to an unstable value along the line of rectilinear diameters. Using these methods to determine the coexisting densities, one can determine many of the thermophysical properties of interest. The properties of interest are: the coexisting densities, the vapor pressure, the enthalpy of vaporization, the critical density, critical pressure, and the critical pressure. Another advantage to using molecular dynamics is the investigation of the interfacial molecule orientation.

In order to study the properties, a
comparison of the two techniques had to be conducted to determine the most
efficient method. Unparameterized methane was simulated to obtain the
coexisting densities. The simulations were carried out in the canonical
ensemble (NVT) with the same number of molecules (8000), for the same length of
equilibration (2 ns). Temperature Quench Molecular Dynamics (TQMD) was
conducted as described as Gelb and Muller [4], and the Volume Expansion Molecular Dynamics (VEMD) was conducted as described by Pamies et al. [5]. VEMD starts with an equilibrated liquid in a cubic simulation cell, and the simulation cell is then expanded to 2.5 times the initial starting configuration. TQMD starts as an equilibrated vapor at a temperature of 220K, and then the temperature was suddenly lowered to the temperature of interest. The system then separates into vapor and liquid phases separated by an interface. The initial starting geometry was 1x1x 2.5 to force the interface into the preferential geometry for the density sampling routine. The coexisting densities were found by cutting the simulation volume into slabs and finding the density profile as a function of the z-axis. This histogram
was then inverted to give the frequency of a given density. This inverted
histogram was used to determine the liquid, vapor, and interfacial densities.
In order to not include interfacial densities with the bulk densities, a
fraction (0.50) of the maximums above and below the average was used to
determine the maximum vapor and minimum liquid densities. Preliminary data for
a polyatomic case, ethylene glycol, will also be shown. To simulate the
polyatomic case the UFF potential[6] was used with multi timescale functionality (RESPA)[7].  

The results were compared to an
accepted equation of state for a Lennard Jones fluid (LJEOS) [8] and to experimental data [9,10]. The comparison was made using the experimental reduced properties from the critical properties determined for each method. The average error of the liquid phase density between the VEMD simulations and the LJEOS is 0.6%.  The average error of the liquid phase density between the VEMD simulations and the experimental data is 1.7%.  The average error of the liquid phase density between the TQMD simulations and the LJEOS is 4.5% and between the TQMD simulations and the experimental data is 2.6%.  As for the vapor phase densities, the average error between the VEMD simulations and the LJEOS is 8.1% and between the VEMD simulations and the experimental data is 11.4%.  For the vapor phase densities, the average error between the TQMD simulations and the LJEOS is 9.2% and between the TQMD simulations and the experimental data is 14.5%. The critical temperature of each method was determined by fitting a power law equation [11] to the difference to the bulk densities as a function of temperature. The critical density was determined by a linear fit of the average density with respect to temperature [12].For VEMD, the resulting critical temperature is T_c =197.6K ± 0.4K.  The values for the critical density were found using the critical temperature and the resulting critical density is r_c = 149.5 kg/m³ ± 0.2
kg/m³. The TQMD model was fit to the same set of equations to give a
critical temperature of T_c =197.5K ± 0.8K and a critical density of r_c = 157.9 kg/m³ ± 0.4
kg/m³  The LJEOS has a critical temperature of T_c =
1.316*e (202.7 K) and a critical
density of r_c = 0.304
molecules/s³ (149.5 kg/m³).
The resulting errors for the critical density and temperature versus the LJEOS
are 2.4% and 0.3% respectively for the VEMD model, and the errors for the TQMD
model are 2.5% for the critical temperature and 5.3% for the critical density.
The errors of the critical temperature and density when compared to
experimental data are 3.7% and 8.0% respectively for the VEMD model. The errors
for the TQMD model when compared to experimental data are 3.6% for the critical
temperature and 2.6% for the critical density. The reported errors in the
properties are from taking into account one standard deviation with the
temperature and density. The results show that there is comparable precision
between the two methods.

To determine the vapor pressure,
single phase simulations were conducted on the vapor densities obtained from
the two phase simulations from the VEMD model. The average relative error of
the vapor pressure between the VEMD simulations and the LJEOS is 3.0%, and the
average relative error between the VEMD simulations and experimental data is
12.4%. The largest deviations of the vapor pressure as compared to experimental
data were at low values of reduced temperature due to the small values obtained.
To determine the critical pressure the vapor pressure values were fit to an
Antoine Equation. The resulting critical pressure is 46.0 bar ± 0.4 bar and the
resulting error for the critical pressure is 11% compared to the LJEOS and 0.1%
when compared to experiment.

The enthalpy of vaporization was
determined by running single phase liquid simulations at the obtained vapor
pressure for a given temperature in the isothermal-isobaric ensemble (NPT). All
of the single phase simulations calculated the enthalpy associated with each
phase, the enthalpy of vaporization was found by subtracting the vapor values
from the liquid values. There is also good agreement with experimental values [10] for methane at all values of the reduced temperature with an average relative error of 2.2%.

In order to compare the
computational efficiency of the two methods, the values of the coexisting
densities were obtained for each step in the simulation. The resulting data was
then fit to simple exponential curves to determine the resulting relaxation
times. For all of the temperatures investigated the VEMD model was faster to
the equilibrated densities. This is especially true at low temperatures where
the equilibration was orders of magnitude faster, but as the temperature
approached the critical temperature the times became closer together.

A new inverted density histogram
technique was tested, and gives good results when compared to the experimental
data and the equation of state. An algorithm is presented to find all of the
properties of interest. The relative efficiency of both TQMD and VEMD was
compared, and VEMD equilibrates faster than TQMD at all temperatures of
interest. Preliminary data was shown for a polyatomic case.

            [1]        A. J. C. Ladd,L. V. Woodcock,
Molecular Physics 36 (1978) 611.

            [2]        G. A. Chapela,G. Saville,S. M. Thompson,J.
S. Rowlinson, Journal of the Chemical Society-Faraday Transactions Ii 73 (1977)
1133.

            [3]        A. Z. Panagiotopoulos, J. Phys.:Condens.
Matter 12 (2000) R25.

            [4]        L. D. Gelb,E. A. Muller, Fluid Phase Equilibria
203 (2002) 1.

            [5]        J. C. Pamies,C. McCabe,P. T. Cummings,L.
F. Vega, Molecular Simulation 29 (2003) 463.

            [6]        A. K. Rappe^',C. J. Casewit,K.
S. Colwell,W. A. G. III,W. M. Skiff, J. Am. Chem. Soc 114 (1992) 10024.

            [7]        M. E. Tuckerman,B. J. Berne,G. J. Martyna,
J. Chem. Phys 97 (1992) 1990.

            [8]        J. K. Johnson,J. A. Zollweg,K. E. Gubbins,
Molecular Physics 78 (1993) 591.

            [9]        U. Setzmann,W. Wagner, J. Phys. Chem.
Ref. Data 20 (1991) 1061.

            [10]      V. Majer,V. Svoboda, Enthalpies of
Vaporization of Organic Compunds: A Critical Review and Data Compilation,
Blackwell Scientific Publications, Oxford, 1985.

            [11]      J. S. Rowlinson,B. Windom, Molecular
Theory of Capillarity, Oxford University Press, New York, 1982.

            [12]      J. S. Rowlinson, Liquids and Liquid
Mixtures, Butterworth & Co, London, 1969.