The free energy of the metastable supersaturated vapor via restricted ensemble simulations. II. Effects of constraints and comparison with molecular dynamics simulations

Extensive restricted canonical ensemble Monte Carlo simulations [D. S. Corti and P. Debenedetti, Chem. Eng. Sci. 49, 2717 (1994)] were performed. Pressure, excess chemical potential, and excess free energy with respect to ideal gas data were obtained at different densities of the supersaturated Lennard-Jones (LJ) vapor at reduced temperatures from 0.7 to 1.0. Among different constraints imposed on the system studied, the one with the local minimum of the excess free energy was taken to be the approximated equilibrium state of the metastable LJ vapor. Also, a comparison of our results with molecular dynamic simulations [A. Linhart et al., J. Chem. Phys. 122, 144506 (2005)] was made.     

Reconstruction of the free energy in the metastable region using the path ensemble

By quenching into the metastable region of the three-dimensional Ising model, we investigate the paths that the magnetization (energy) takes as a function of time. We accumulate the magnetization (energy) paths into time-dependent distributions from which we reconstruct the free energy as a function of the magnetic field, temperature, and system size. From the reconstructed free energy, we obtain the free-energy barrier that is associated with the transition from a metastable state to the stable equilibrium state. Although mean-field theory predicts a sharp transition between the metastable and the unstable region where the free-energy barrier is zero, the results for the nearest-neighbor Ising model show that the free-energy barrier does not go to zero.     

Thermal properties of the metastable supersaturated vapor of the Lennard-Jones fluid

p, rho, T data of the supersaturated vapor of the Lennard-Jones fluid are obtained by molecular dynamics simulations. The metastable state points are identified before a phase separation takes place. An estimation of the location of the spinodal is given. The results are compared to two theoretically based equations of state and one empirical equation of state which was parametrized also taking into account metastable state points. The pressure obtained by simulation is found to be lower than that from both theoretically based equations of state, which do not account for the inhomogeneous density distribution of the supersaturated vapor.     

Molecular-dynamics simulation of argon nucleation from supersaturated vapor in the NVE ensemble

The possibility to conduct simulations of homogeneous nucleation of argon from a supersaturated vapor phase using a microcanonical or NVE ensemble is evaluated (NVE: number of particles N, volume V, and energy E are constant). In order to initiate a phase separation kinetic energy is removed from the system in one step which transfers the system into a supersaturated state. After this temperature jump the simulation is continued in a NVE ensemble. The simulations are performed for different initial-state points and different temperature jumps. The cluster formation and growth over the course of the adiabatic simulations are analyzed. The progression of the temperature being related to the cluster size in NVE systems is traced. Also the influence of the size of the simulation system is investigated. For a certain range of low supersaturation a dynamic coexistence between two states has been found. Furthermore, the obtained nucleation rates are correlated with two simple functions. By applying the nucleation theorems to these functions the size and excess energy of the critical cluster are estimated. The results are consistent with other theoretical data and experimental data available in the literature.     

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.     

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.     

Toward efficient chemical potential calculations by expanded ensemble simulations; to make the free energy pathway fairly level

A scheme is suggested of how to construct good bias potentials ("balancing factors") to be used in expanded ensemble (EE) calculations of chemical potentials of solutions. A combination of two strategies are used: (i) to use a pathway for particle insertions that avoids large variations in free energy and (ii) to use calculated free energy derivatives to construct a bias potential that makes the pathway fairly level. Only a few very short simulations are needed to accomplish the latter, and then, a full EE simulation is done to obtain the chemical potential. By practical calculations of the chemical potential of benzene, cyclohexane, and benzylamine in water, it is shown that this method is at least equally efficient to the recent adaptive EE (AEE) method by Aberg et al. (J. Chem. Phys. 2004, 120, 3370). Furthermore, the new method provides an alternative strategy that complements existing EE methods.     

Dissolution of solid NaCl nanoparticles embedded in supersaturated water vapor probed by molecular dynamic simulations

The dissolution process for small, on the order of 1000 atoms, crystalline NaCl particles with defects embedded in highly supersaturated water vapor was studied by the molecular dynamics (MD) simulation method. We found that a breakdown of the crystal lattice does not occur unless (1) the thickness of water layer covering the surface of salt particles exceeds several molecular layers and (2) there are a considerable number of defects in the crystal. The collapse of the crystal lattice starts when the amount of water taken up by a salt particle reaches about half ( approximately 50%) of the amount of salt in this particle. The number of defects required to initiate subsequent dissolution of the NaCl crystal on the time scale accessible by our simulations ( approximately 40 ns) is in the range of 10 to 12%. We also report the estimates for the time required to form supersaturated aqueous solutions of NaCl from originally crystalline particles as a function of the number of defects in the crystal.     

Calculating Gibbs free energies of transfer from Gibbs ensemble Monte Carlo simulations

The partitioning of the ternary systems n-pentane/n-heptane/(helium or argon) at ambient conditions is investigated using configurational-bias Monte Carlo simulations in the Gibbs ensemble. The results demonstrate that this approach yields very precise partition constants and free energies of transfer. Simulations are carried out to study the dependence of the n-pentane partitioning with respect to the carrier gas, the system size, and the overall solute concentrations. None of the changes of variables, within the ranges used here, has a significant effect on the alkane partitioning. However, chemical potentials calculated via Widom's ghost particle insertions show a strong number dependence for phases containing relatively few molecules of a given type. This problem originates from the fact that the chemical potential is calculated for a concentration of real particles plus one ghost particle that is systematically larger than the equilibrium concentration. A simple correction term is suggested to account for this problem. Received: 13 May 1998 / Accepted: 18 June 1998 / Published online: 4 September 1998     

Combinational simulations of free energy of polymer solution

It is conceptually proposed that the total entropy of polymer solution is contributed from two distinct parts: the positional and the oomformational. The former can be represented analytically, while the latter can be simulated with the random self-avoiding walk model on the simple cubic lattice for multichain systems. The obtained results indicated that both the conformational entropy and the mixing heat are consistent with the scaling laws wry well.     

Grain-boundary free energy via thermodynamic integration

In a previous publication by Lusk and Beale [Phys. Rev. E 69, 026117 (2004)], fluctuating cell (FC) theory was used to estimate the free energy of symmetric tilt grain boundaries in an assembly of nearly hard disks. The FC method is much faster than the more traditional thermodynamic integration, but the accuracy of the algorithm has not been assessed in association with persistent defect structures. This motivated the present work wherein the FC free energies are compared directly with the data obtained via thermodynamic integration from an Einstein crystal to an assembly of hard disks. This comparison is made over the range of possible misorientations for symmetric tilt boundaries and indicates that the FC method gives quantitatively accurate estimates for grain-boundary free energy. We also demsonstrate that the FC approximation is quantitatively accurate at determining the free-energy contribution of each particle whether in the bulk or the grain boundary. The FC calculation is about two orders of magnitude faster than a full thermodynamic integration. This approach may offer a numerically efficient means of estimating the free energy of persistent defect structures to greater accuracy than is afforded by the quasiharmonic and local harmonic approximations.     

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.     

Kinetics of supersaturated vapor nucleation on molecular condensation nuclei

The kinetics of nucleation is calculated for a supersaturated vapor containing molecular condensation nuclei, that is, foreign molecules able to induce the formation of viable nuclei of a condensed phase by themselves. In contrast to the previous calculation, the possibility of the escape of molecular condensation nuclei from very small clusters containing a few condensed vapor molecules is taken into account. More exact equations are derived for the rate of steady-state nucleation and the concentration of aerosol particles in a quasisteady-state regime of nucleation. The calculation demonstrates that, at a high probability of the escape of a molecular condensation nucleus, the predominating mechanism of cluster formation is the attachment of a molecular condensation nucleus to a cluster formed from vapor molecules rather than their condensation on the nucleus. At the same time, allowances for the possible escape of molecular condensation nuclei from clusters slightly affect the rate of nucleation and the concentration of aerosol particles being formed.     

Minimum free energy pathways and free energy profiles for conformational transitions based on atomistic molecular dynamics simulations

An efficient method for the calculation of minimum free energy pathways and free energy profiles for conformational transitions is presented. Short restricted perturbation-targeted molecular dynamics trajectories are used to generate an approximate free energy surface. Approximate reaction pathways for the conformational change are constructed from one-dimensional line segments on this surface using a Monte Carlo optimization. Accurate free energy profiles are then determined along the pathways by means of one-dimensional adaptive umbrella sampling simulations. The method is illustrated by its application to the alanine "dipeptide." Due to the low computational cost and memory demands, the method is expected to be useful for the treatment of large biomolecular systems.     

Computation of methodology-independent ionic solvation free energies from molecular simulations. II. The hydration free energy of the sodium cation

The raw ionic solvation free energies computed from atomistic (explicit-solvent) simulations are extremely sensitive to the boundary conditions (finite or periodic system, system shape, and size) and treatment of electrostatic interactions (Coulombic, lattice sum, or cutoff based) used during these simulations. In the present article, it is shown that correction terms can be derived for the effect of (A) an incorrect solvent polarization around the ion due to the use of an approximate (not strictly Coulombic) electrostatic scheme; (B) the finite size or artificial periodicity of the simulated system; (C) an improper summation scheme to evaluate the potential at the ion site and the possible presence of a liquid-vacuum interface in the simulated system. Taking the hydration free energy of the sodium cation as a test case, it is shown that the raw solvation free energies obtained using seven different types of boundary conditions and electrostatic schemes commonly used in explicit-solvent simulations (for a total of 72 simulations differing in the corresponding simulation parameters) can be corrected so as to obtain a consistent value for this quantity.     

Homogeneous nucleation of droplets from supersaturated vapor in a closed system

Kinetic equations describing homogeneous nucleation kinetics within standard model are solved numerically under the condition of a constant number of molecules in the considered system. It has consequences to decrease the supersaturation of the supersaturated vapor during the process of the formation of small droplets of a new phase. The decrease of supersaturation occurs in a short time and reaches some value which remains unchanged for a relatively long time (quasistationary regime), especially at lower initial supersaturations. This time interval decreases with increasing value of the initial supersaturation. In the quasistationary regime the nucleation rate reaches its stationary value. At higher initial supersaturation, the rate of formation of nuclei goes to some maximum value corresponding to the stationary nucleation rate and then decreases with time due to the decrease of supersaturation.     

Computation of surface tensions using expanded ensemble simulations

A method for the direct simulation of the surface tension is examined. The technique is based on the thermodynamic route to the interfacial tension and makes use of the expanded ensemble simulation method for the calculation of the free energy difference between two inhomogeneous systems with the same number of particles, temperature, and volume, but different interfacial area. The method is completely general and suitable for systems with either continuous or discontinuous interactions. The adequacy of the expanded ensemble method is assessed by computing the interfacial tension of the planar vapor-liquid interface of Lennard-Jones, Lennard-Jones dimers, Gay-Berne, and square-well model fluids; in the latter, the interactions are discontinuous and the present method does not exhibit the asymmetry of other related methods, such as the test area. The expanded ensemble simulation results are compared with simulation data obtained from other techniques (mechanical and test area) with overall good agreement.     

Hiking down the energy landscape: progress toward the Kauzmann temperature via vapor deposition

Physical vapor deposition was employed to prepare amorphous samples of indomethacin and 1,3,5-(tris)naphthylbenzene. By depositing onto substrates held somewhat below the glass transition temperature and varying the deposition rate from 15 to 0.2 nm/s, glasses with low enthalpies and exceptional kinetic stability were prepared. Glasses with fictive temperatures that are as much as 40 K lower than those prepared by cooling the liquid can be made by vapor deposition. As compared to an ordinary glass, the most stable vapor-deposited samples moved about 40% toward the bottom of the potential energy landscape for amorphous materials. These results support the hypothesis that enhanced surface mobility allows stable glass formation by vapor deposition. A comparison of the enthalpy content of vapor-deposited glasses with aged glasses was used to evaluate the difference between bulk and surface dynamics for indomethacin; the dynamics in the top few nanometers of the glass are about 7 orders of magnitude faster than those in the bulk at Tg - 20 K.     

Molecular analysis of small drops in a metastable vapor

The molecular theory of curved vapor-liquid interfaces within the lattice gas model is applied to analyze supersaturated vapor states in dependence on the new phase size and system temperature. The molecular interaction is considered in the quasi-chemical approximation which describes effects of direct correlations of the nearest molecules. Two methods for determining the surface tension are discussed: equimolecular and according to the surface tension minimum in the intermediate region, i.e., on the tension surface. It is shown that a tension surface exists for metastable drops in supersaturated vapor over the temperature range, but its use leads to multivaluedness of solutions for its position; the minimum range of the existence of metastable drops is close to the previously determined lower limit for equilibrium.