The free energy of the metastable supersaturated vapor via restricted ensemble simulations

Pressure, excess chemical potential, and excess free energy, with respect to ideal gas data at different densities of the supersaturated Lennard-Jones particle vapor at the reduced temperature 0.7 are obtained by the restricted canonical ensemble Monte Carlo simulation method [D. S. Corti and P. Debenedetti, Chem. Eng. Sci. 49, 2717 (1994)]. The excess free energy values depend upon the constraints imposed on the system with local minima exhibited for densities below the spinodal density and monotonic variation for densities larger than the spinodal density. The results are compared with a molecular dynamics simulation study [A. Linharton et al., J. Chem. Phys. 122, 144506 (2005)] on the same system. The current study verifies the conclusion drawn by the simulation work that clustering of Lennard-Jones atoms exists even in the vicinity of spinodal. Our method gives an alternative to molecular dynamic simulations for the determination of equilibrium properties of a metastable fluid, especially close to the spinodal, and does not require a very large system to carry out the simulation.     

Monte Carlo simulations of critical cluster sizes and nucleation rates of water

We have calculated the critical cluster sizes and homogeneous nucleation rates of water at temperatures and vapor densities corresponding to experiments by Wolk and Strey [J. Phys. Chem B 105, 11683 (2001)]. The calculations have been done with an expanded version of a Monte Carlo method originally developed by Vehkamaki and Ford [J. Chem. Phys. 112, 4193 (2000)]. Their method calculates the statistical growth and decay probabilities of molecular clusters. We have derived a connection between these probabilities and kinetic condensation and evaporation rates, and introduce a new way for the calculation of the work of formation of clusters. Three different interaction potential models of water have been used in the simulations. These include the unpolarizable SPC/E [J. Phys. Chem. 91, 6269 (1987)] and TIP4P [J. Chem. Phys. 79, 926 (1983)] models and a polarizable model by Guillot and Guissani [J. Chem. Phys. 114, 6720 (2001)]. We show that TIP4P produces critical cluster sizes and a temperature and vapor density dependence for the nucleation rate that agree well with the experimental data, although the magnitude of nucleation rate is constantly overestimated by a factor of 2 x 10(4). Guissani and Guillot's model is somewhat less successful, but both the TIP4P and Guillot and Guissani models are able to reproduce a much better experimental temperature dependency of the nucleation rate than the classical nucleation theory. Using SPC/E results in dramatically too small critical clusters and high nucleation rates. The water models give different average binding energies for clusters. We show that stronger binding between cluster molecules suppresses the decay probability of a cluster, while the growth probability is not affected. This explains the differences in results from different water models.     

Study of thermal properties of the metastable supersaturated vapor with the integral equation method

Pressure, excess chemical potential, and excess free energy data for different densities of the supersaturated argon vapor at reduced temperatures from 0.7 to 1.2 are obtained by solving the integral equation with perturbation correction to the radial distribution function [F. Lado, Phys. Rev. 135, A1013 (1964)]. For those state points where there is no solution, the integral equation is solved with the interaction between argon atoms modeled by Lennard-Jones potential plus a repulsive potential with one controlling parameter, alpha exp(-rsigma) and in the end, all the thermal properties are mapped back to the alpha=0 case. Our pressure data and the spinodal obtained from the current method are compared with a molecular dynamics simulation study [A. Linhart et al., J. Chem. Phys. 122, 144506 (2005)] of the same system.     

Thermodynamic properties of van der Waals fluids from Monte Carlo simulations and perturbative Monte Carlo theory

Computer simulations have been performed for fluids with van der Waals potential, that is, hard spheres with attractive inverse power tails, to determine the equation of state and the excess energy. On the other hand, the first- and second-order perturbative contributions to the energy and the zero- and first-order perturbative contributions to the compressibility factor have been determined too from Monte Carlo simulations performed on the reference hard-sphere system. The aim was to test the reliability of this "exact" perturbation theory. It has been found that the results obtained from the Monte Carlo perturbation theory for these two thermodynamic properties agree well with the direct Monte Carlo simulations. Moreover, it has been found that results from the Barker-Henderson [J. Chem. Phys. 47, 2856 (1967)] perturbation theory are in good agreement with those from the exact perturbation theory.     

Chemical potential perturbation: a method to predict chemical potentials in periodic molecular simulations

A new method, called chemical potential perturbation (CPP), has been developed to predict the chemical potential as a function of density in periodic molecular simulations. The CPP method applies a spatially varying external force field to the simulation, causing the density to depend upon position in the simulation cell. Following equilibration the homogeneous (uniform or bulk) chemical potential as a function of density can be determined relative to some reference state after correcting for the effects of the inhomogeneity of the system. We compare three different methods of approximating this correction. The first method uses the van der Waals density gradient theory to approximate the inhomogeneous Helmholtz free energy density. The second method uses the local pressure tensor to approximate the homogeneous pressure. The third method uses the Triezenberg-Zwanzig definition of surface tension to approximate the inhomogeneous free energy density. If desired, the homogeneous pressure and Helmholtz free energy can also be predicted by the new method, as well as binodal and spinodal densities of a two-phase fluid region. The CPP method is tested using a Lennard-Jones (LJ) fluid at vapor, liquid, two-phase, and supercritical conditions. Satisfactory agreement is found between the CPP method and an LJ equation of state. The efficiency of the CPP method is compared to that for Widom's method under the tested conditions. In particular, the new method works well for dense fluids where Widom's method starts to fail.     

Escorted free energy simulations

We describe a strategy to improve the efficiency of free energy estimates by reducing dissipation in nonequilibrium Monte Carlo simulations. This strategy generalizes the targeted free energy perturbation approach [C. Jarzynski, Phys. Rev. E 65, 046122 (2002)] to nonequilibrium switching simulations, and involves generating artificial, "escorted" trajectories by coupling the evolution of the system to updates in external work parameter. Our central results are: (1) a generalized fluctuation theorem for the escorted trajectories, and (2) estimators for the free energy difference ΔF in terms of these trajectories. We illustrate the method and its effectiveness on model systems.     

Investigation of the salting out of methane from aqueous electrolyte solutions using computer simulations

We calculate the excess chemical potential of methane in aqueous electrolyte solutions of NaCl using Monte Carlo computer simulations. In a recent work [Docherty et al. J. Chem. Phys. 2006, 125, 074510], we presented a new potential model for methane in water which is capable of describing accurately the excess chemical potential of methane in pure water over a range of temperatures, a quantity that can be related to the solubility and which is commonly used to study the hydrophobic effect. Here, we use the same potential model for the water-methane interactions and investigate the effect of added salt on the chemical potential of methane in the solution. The methane molecules are modeled as single Lennard-Jones (LJ) interaction sites, and the water molecules are modeled with the TIP4P/2005 model. A correcting factor of chi = 1.07 for the energetic Berthelot (geometric) combining rule of the methane-water interaction is also used, which mimics the polarization of methane in water. We consider NaCl as the salt and treat the ions with the Smith and Dang model (i.e., as charged LJ interaction sites). Ion-water, ion-ion, and ion-methane interactions are treated using Lorentz-Berthelot combining rules. In addition, the Coulombic potential is used to model charge-charge interactions which are calculated using the Ewald sum. We have carried out isobaric-isothermal (NpT) simulations to determine the equilibrium densities of the solutions. The simulation data is in excellent agreement with experimental densities of aqueous NaCl solutions of different concentration. Hydration numbers are also obtained and found to be in agreement with reported data. Canonical (NVT) simulations at the averaged densities are then performed using the Widom test-particle insertion method to obtain the excess chemical potential of methane in the saline solutions. An increase in the chemical potential of methane, corresponding to a salting out effect, is observed when salt is added to the solution. We investigate different concentrations and ion sizes. An overprediction of the salting out effect as compared with experimental data is observed, which we believe is due to the polarizing effect of the ions in the solution, which is not taken into account by the model. We also find a direct correlation between the increase in the chemical potential and the packing fraction of the solution and argue that the main cause of the observed salting out effect (as represented by an increase in the excess chemical potential) is the increase in the packing fraction of the solutions due to the added salt. Together, with this, we put forward an argument toward explaining the anomalous Hofmeister effect of Li(+).     

Investigation of excess adsorption, solvation force, and plate-fluid interfacial tension for Lennard-Jones fluid confined in slit pores

The excess Helmholtz free energy functional is formulated in terms of a modified fundamental measure theory [Y. X. Yu and J. Z. Wu, J. Chem. Phys. 117, 10156 (2002)] for a short ranged repulsion and a first-order mean-spherical approximation theory [Y. P. Tang, J. Chem. Phys. 118, 4140 (2003)] for a long ranged attraction. Within the framework of the density functional theory, the density profile, excess adsorption, solvation force, and plate-fluid interfacial tension of a Lennard-Jones fluid confined in slit pores are predicted, and the results agree well with the simulation data. The phase equilibria inside the slit pores are determined according to the requirement that temperature, chemical potential, and grand potential in coexistence phases should be equal, and the plate-fluid interfacial tensions at equilibrium states are predicted consequently.     

Phase coexistence in polydisperse multi-Yukawa hard-sphere fluid: high temperature approximation

High temperature approximation (HTA) is used to describe the phase behavior of polydisperse multi-Yukawa hard-sphere fluid mixtures. It is demonstrated that in the frames of the HTA the model belongs to the class of "truncatable free energy models," i.e., the models with thermodynamical properties (Helmholtz free energy, chemical potential, and pressure) defined by the finite number of generalized moments. Using this property we were able to calculate the complete phase diagram (i.e., cloud and shadow curves as well as binodals) and size distribution functions of the coexisting phases of several different models of polydisperse fluids. In particular, we consider polydisperse one-Yukawa hard-sphere mixture with factorizable Yukawa coefficients and polydisperse Lennard-Jones (LJ) mixture with interaction energy parameter and/or size polydispersity. To validate the accuracy of the HTA we compare theoretical results with previously published results of more advanced mean spherical approximation (MSA) for the one-Yukawa model and with the Monte Carlo (MC) computer simulation results of [Wilding et al. J. Chem. Phys. 121, 6887 (2004); Phys. Rev. Lett. 95, 155701 (2005)] for the LJ model. We find that overall predictions of the HTA are in reasonable agreement with predictions of the MSA and MC, with the accuracy range from semiquantitative (for the phase diagram) to quantitative (for the size distribution functions).     

Numerical verification of the generalized Crooks nonequilibrium work theorem for non-Hamiltonian molecular dynamics simulations

The generalized Crooks theorem (GCT) for deterministic non-Hamiltonian molecular dynamics simulations [Phys. Rev. E 75, 050101 (2007)] connects the probabilities of nonequilibrium realizations switching the system between two thermodynamic states, to the partition functions of these states. In comparison to the "classical" Crooks nonequilibrium work theorem [J. Stat. Phys. 90, 1481 (1998)], which deals with realizations involving only mechanical work, the GCT also accounts for additional work resulting from changes of the intensive and extensive thermodynamic variables of the system. In this article we present a numerical verification of the GCT using a Lennard-Jones fluid model where two particles are subject to a time-dependent external potential. Moreover, in order to switch the system between different thermodynamic states, the temperature and the pressure (or volume), which are controlled through the Martyna-Tobias-Klein equations of motion [J. Chem. Phys. 101, 4177 (1994)], are also varied externally. The free energy difference between states characterized by different distances of the target particles is evaluated using both a standard methodology (pair radial distribution functions) and the GCT. In order to exploit the various options provided by the GCT approach, i.e., the possibility of temperature/pressure/volume changes during the realizations, the free energy difference is recovered via arbitrary thermodynamic cycles. In all tests, the GCT is quantitatively verified.     

Communication: Thermodynamics of water modeled using ab initio simulations

We regularize the potential distribution framework to calculate the excess free energy of liquid water simulated with the BLYP-D density functional. Assuming classical statistical mechanical simulations at 350 K model the liquid at 298 K, the calculated free energy is found in fair agreement with experiments, but the excess internal energy and hence also the excess entropy are not. The utility of thermodynamic characterization in understanding the role of high temperatures to mimic nuclear quantum effects and in evaluating ab initio simulations is noted.     

Molecular dynamics simulation of potassium along the liquid-vapor coexistence curve

The applicability of pair potential functions to liquid alkali metals is questionable. On the one hand, some recent reports in the literature suggest the validity of two-parameter pair-wise additive Lennard-Jones (LJ) potentials for liquid alkali metals. On the other hand, there are some reports suggesting the inaccuracy of pair potential functions for liquid metals. In this work, we have performed extensive molecular dynamics simulations of vapor-liquid phase equilibria in potassium to check the validity of the proposed LJ potentials and to improve their accuracy by changing the LJ exponents and taking into account the temperaturedependencies of the potential parameters. We have calculated the orthobaric liquid and vapor densities of potassium using LJ (12–6), LJ (8.5–4) and LJ (5–4), effective pair potential energy functions. The results show that using an LJ (8.5–4) potential energy function with temperature-independent parameters, ε and σ, is inadequate to account for the vapor-liquid coexistence properties of potassium. Taking into account the temperature-dependencies of the LJ parameters, ε(T) and σ(T), we obtained the densities of coexisting liquid and vapor potassium in a much better agreement with experimental data. Changing the magnitude of repulsive and attractive contributions to the potential energy function shows that a two-parameter LJ (5–4) potential can well reproduce the densities of liquid and vapor potassium. The results show that LJ (5–4) potential with temperature-dependent parameters produces the densities of liquid and vapor potassium more accurately, compared to the results obtained using LJ (12–6) and LJ (8.5–4) potential energy functions.     

Overcoming free energy barriers using unconstrained molecular dynamics simulations

Association of unconstrained molecular dynamics (MD) and the formalisms of thermodynamic integration and average force [Darve and Pohorille, J. Chem. Phys. 115, 9169 (2001)] have been employed to determine potentials of mean force. When implemented in a general MD code, the additional computational effort, compared to other standard, unconstrained simulations, is marginal. The force acting along a chosen reaction coordinate xi is estimated from the individual forces exerted on the chemical system and accumulated as the simulation progresses. The estimated free energy derivative computed for small intervals of xi is canceled by an adaptive bias to overcome the barriers of the free energy landscape. Evolution of the system along the reaction coordinate is, thus, limited by its sole self-diffusion properties. The illustrative examples of the reversible unfolding of deca-L-alanine, the association of acetate and guanidinium ions in water, the dimerization of methane in water, and its transfer across the water liquid-vapor interface are examined to probe the efficiency of the method.     

An application of coupled reference interaction site model/molecular dynamics to the conformational analysis of the alanine dipeptide

We present an application of our recently proposed coupled reference interaction site model (RISM) molecular dynamics (MD) solvation free energy methodology [Freedman and Truong, Chem. Phys. Lett. 381, 362 (2003); J. Chem. Phys. 121, 2187 (2004)] to study the conformational stability of alanine dipeptide in aqueous solution. In this methodology, radial distribution functions obtained from a single MD simulation are substituted into a RISM expression for solvation free energy. Consequently, iterative solution of the RISM equation is not needed. The relative solvation free energies of seven different conformations of the alanine dipeptide in aqueous solution are calculated. Results from the coupled RISM/MD methodology are in good agreement with those from earlier simulations using the accurate free energy perturbation approach, showing that the alphaR conformation is most stabilized by solution. This study establishes a framework for applying this coupled RISM/MD method to larger biological systems.     

Computational studies on thermodynamic properties, effective diameters, and free volume of argon using an ab initio potential

A quantum mechanical derived ab initio interaction potential for the argon dimer was tested in molecular simulations to reproduce the thermophysical properties of the vapor-liquid phase equilibria using the Gibbs ensemble Monte Carlo simulations as well as the liquid and supercritical equation of state using the NVT Monte Carlo simulations. The ab initio interaction potential was taken from the literature. A recently developed theory [R. Laghaei et al., J. Chem. Phys. 124, 154502 (2006)] was used to compute the effective diameters of argon in fluid phases and the results were subsequently applied in the generic van der Waals theory to compute the free volume of argon. The calculated densities of the coexisting phases, the vapor pressure, and the equation of state show excellent agreement with experimental values. The effective diameters and free volumes of argon are given over a wide range of densities and temperatures. An empirical formula was used to fit the effective diameters as a function of density and temperature. The computed free volume will be used in future investigations to calculate the transport properties of argon.     

The nature of the calculation of the pressure in molecular simulations of continuous models from volume perturbations

We consider some fundamental aspects of the calculation of the pressure from simulations by performing volume perturbations. The method, initially proposed for hard-core potentials by Eppenga and Frenkel [Mol. Phys.52, 1303 (1984)] and then extended to continuous potentials by Harismiadis et al. [J. Chem. Phys. 105, 8469 (1996)], is based on the numerical estimate of the change in Helmholtz free energy associated with the perturbation which, in turn, can be expressed as an ensemble average of the corresponding Boltzmann factor. The approach can be easily generalized to the calculation of components of the pressure tensor and also to ensembles other than the canonical ensemble. The accuracy of the method is assessed by comparing simulation results obtained from the volume-perturbation route with those obtained from the usual virial expression for several prototype fluid models. Monte Carlo simulation data are reported for bulk fluids and for inhomogeneous systems containing a vapor-liquid interface.     

Heat capacity effects associated with the hydrophobic hydration and interaction of simple solutes: a detailed structural and energetical analysis based on molecular dynamics simulations

We examine the SPCE [H. J. C. Berendsen et al., J. Chem. Phys. 91, 6269 (1987)] and TIP5P [M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys 112, 8910 (2000)] water models using a temperature series of molecular-dynamics simulations in order to study heat-capacity effects associated with the hydrophobic hydration and interaction of xenon particles. The temperature interval between 275 and 375 K along the 0.1-MPa isobar is studied. For all investigated models and state points we calculate the excess chemical potential for xenon employing the Widom particle insertion technique. The solvation enthalpy and excess heat capacity is obtained from the temperature dependence of the chemical potentials and, alternatively, directly by Ewald summation, as well as a reaction field based method. All three methods provide consistent results. In addition, the reaction field technique allows a separation of the solvation enthalpy into solute/solvent and solvent/solvent parts. We find that the solvent/solvent contribution to the excess heat capacity is dominating, being about one order of magnitude larger than the solute/solvent part. This observation is attributed to the enlarged heat capacity of the water molecules in the hydration shell. A detailed spatial analysis of the heat capacity of the water molecules around a pair of xenon particles at different separations reveals that even more enhanced heat capacity of the water located in the bisector plane between two adjacent xenon atoms is responsible for the maximum of the heat capacity found for the desolvation barrier distance, recently reported by Shimizu and Chan [J. Am. Chem. Soc. 123, 2083 (2001)]. The about 60% enlarged heat capacity of water in the concave part of the joint xenon-xenon hydration shell is the result of a counterplay of strengthened hydrogen bonds and an enhanced breaking of hydrogen bonds with increasing temperature. Differences between the two models with respect to the heat capacity in the xenon-xenon contact state are attributed to the different water model bulk heat capacities, and to the different spatial extension of the structure effect introduced by the hydrophobic particles. Similarities between the different states of water in the joint xenon-xenon hydration shell and the properties of stretched water are discussed.     

Pressure calculation in hybrid particle-field simulations

In the framework of a recently developed scheme for a hybrid particle-field simulation techniques where self-consistent field (SCF) theory and particle models (molecular dynamics) are combined [J. Chem. Phys. 130, 214106 (2009)], we developed a general formulation for the calculation of instantaneous pressure and stress tensor. The expressions have been derived from statistical mechanical definition of the pressure starting from the expression for the free energy functional in the SCF theory. An implementation of the derived formulation suitable for hybrid particle-field molecular dynamics-self-consistent field simulations is described. A series of test simulations on model systems are reported comparing the calculated pressure with those obtained from standard molecular dynamics simulations based on pair potentials.     

Monte Carlo simulation study of droplet nucleation

A new rigorous Monte Carlo simulation approach is employed to study nucleation barriers for droplets in Lennard-Jones fluid. Using the gauge cell method we generate the excess isotherm of critical clusters in the size range from two to six molecular diameters. The ghost field method is employed to compute the cluster free energy and the nucleation barrier with desired precision of (1-2)kT. Based on quantitative results obtained by Monte Carlo simulations, we access the limits of applicability of the capillarity approximation of the classical nucleation theory and the Tolman equation. We show that the capillarity approximation corrected for vapor nonideality and liquid compressibility provides a reasonable assessment for the size of critical clusters in Lennard-Jones fluid; however, its accuracy is not sufficient to predict the nucleation barriers for making practical estimates of the rate of nucleation. The established dependence of the droplet surface tension on the droplet size cannot be approximated by the Tolman equation for small droplets of radius less than four molecular diameters. We confirm the conclusion of ten Wolde and Frenkel [J. Chem. Phys. 109, 9901 (1998)] that integration of the normal component of the Irving-Kirkwood pressure tensor severely underestimates the nucleation barriers for small clusters.