On the direct calculation of the free energy of quantization for molecular systems in the condensed phase

Using the path integral formalism or the Feynman-Hibbs approach, various expressions for the free energy of quantization for a molecular system in the condensed phase can be derived. These lead to alternative methods to directly compute quantization free energies from molecular dynamics computer simulations, which were investigated with an eye to their practical use. For a test system of liquid neon, two methods are shown to be most efficient for a direct evaluation of the excess free energy of quantization. One of them makes use of path integral simulations in combination with a single-step free energy perturbation approach and was previously reported in the literature. The other method employs a Feynman-Hibbs effective Hamiltonian together with the thermodynamic integration formalism. However, both methods are found to give less accurate results for the excess free energy of quantization than the estimate obtained from explicit path integral calculations on the excess free energy of the neon liquid in the classical and quantum mechanical limit. Suggestions are made to make both methods more accurate.     

Molecular dynamics simulations of atmospheric oxidants at the air-water interface: solvation and accommodation of OH and O3

A comparative study of OH, O3, and H2O equilibrium aqueous solvation and gas-phase accommodation on liquid water at 300 K is performed using a combination of ab initio calculations and molecular dynamics simulations. Polarizable force fields are developed for the interaction potential of OH and O3 with water. The free energy profiles for transfer of OH and O3 from the gas phase to the bulk liquid exhibit a pronounced minimum at the surface, but no barrier to solvation in the bulk liquid. The calculated surface excess of each oxidant is comparable to calculated and experimental values for short chain, aliphatic alcohols. Driving forces for the surface activity are discussed in terms of the radial distribution functions and dipole orientation distributions for each molecule in the bulk liquid and at the surface. Simulations of OH, O3, and H2O impinging on liquid water with a thermal impact velocity are used to calculate thermal accommodation (S) and mass accommodation (alpha) coefficients. The values of S for OH, O3, and H2O are 0.95, 0.90, and 0.99, respectively. The approaching molecules are accelerated toward the liquid surface when they are approximately 5 angstroms above it. The molecules that reach thermal equilibrium with the surface do so within 2 ps of striking the surface, while those that do not scatter into the gas phase with excess translational kinetic energy in the direction perpendicular to the surface. The time constants for absorption and desorption range from approximately 35 to 140 ps, and the values of alpha for OH, O3, and H2O are 0.83, 0.047, and 0.99, respectively. The results are consistent with previous formulations of gas-phase accommodation from simulations, in which the process occurs by rapid thermal and structural equilibration followed by diffusion on the free energy profile. The implications of these results with respect to atmospheric chemistry are discussed.     

Pressure and temperature dependence of hydrophobic hydration: volumetric, compressibility, and thermodynamic signatures

The combined effect of pressure and temperature on hydrophobic hydration of a nonpolar methanelike solute is investigated by extensive simulations in the TIP4P model of water. Using test-particle insertion techniques, free energies of hydration under a range of pressures from 1 to 3000 atm are computed at eight temperatures ranging from 278.15 to 368.15 K. Corresponding enthalpy, entropy, and heat capacity accompanying the hydration process are estimated from the temperature dependence of the free energies. Partial molar and excess volumes calculated using pressure derivatives of the simulated free energies are consistent with those determined by direct volume simulations; but direct volume determination offers more reliable estimates for compressibility. At 298.15 K, partial molar and excess isothermal compressibilities of methane are negative at 1 atm. Partial molar and excess adiabatic (isentropic) compressibilities are estimated to be also negative under the same conditions. But partial molar and excess isothermal compressibilities are positive at high pressures, with a crossover from negative to positive compressibility at approximately 100-1000 atm. This trend is consistent with experiments on aliphatic amino acids and pressure-unfolded states of proteins. For the range of pressures simulated, hydration heat capacity exhibits little pressure dependence, also in apparent agreement with experiment. When pressure is raised at constant room temperature, hydration free energy increases while its entropic component remains essentially constant. Thus, the increasing unfavorability of hydration under raised pressure is seen as largely an enthalpic effect. Ramifications of the findings of the authors for biopolymer conformational transitions are discussed.     

Elastic properties of the (001) face of xenon crystals

In this work we investigated the elastic properties of the (001) face of xenon crystal. The slabs (twodimensional crystals) defined by (001) planes are generated, their structures are optimized and the slabs thermodynamic functions in excess to the crystal bulk calculated. The calculations are based on the Lennard-Jones 6?12 force field, classical elasticity theory and surface thermodynamics. In this work, the number of planes undergoing relaxation is not a priori constrained but it follows from the minimization of the free energy of the slabs and of the bulk, in respect to atomic positions. The value of the surface free energy is calculated as a function of the homogeneous strain of the 2D (001) cell measured relatively to the cell of the stable 3D crystal. At 0 K, when strain is not applied, the specific surface free energy is about 0.064 Jm^-2 and decreases by about 6% at 50 K. The surface stress is positive amounting to 0.010 Jm^-2 at 0 K, and it decreases by about 50% at 50 K. We find that the surface stress can be released by a reorganization of the interatomic distances at the crystal surfaces. The surface excess mean value of the slab elastic constants at 0 K is small (0.012 GPa) and it decreases by about 35% at 50 K. The method proposed can be alternative to molecular dynamics simulations in order to assess the excess surface properties of materials having a complex structure.     

Isobaric vapour–liquid phase diagram and excess properties for the binary system 1,4-dioxane + water at 298.15 K, 318.15 K and 338.15 K

Isobaric vapour–liquid equilibrium (VLE) data have been measured for the binary system (1,4-dioxane + water) at atmospheric pressure. The compositions of the two phases in equilibrium were taken from the refractive indices curve of the system. The equilibrium data show a positive deviation from ideality and reveal an azeotrope. VLE data have been used to calculate the activity coefficients of the two components at 298.15 K, 318.15 K and 338.15 K. The excess properties have been deduced at these temperatures. Molar Gibbs free energy at liquid state was expressed as a function of composition and temperature. At the liquid state, excess entropies of the present system are negative over the entire composition range. Molecular-scale segregation in this system was proved.     

Free energy of solvation of simple ions: molecular-dynamics study of solvation of Cl- and Na+ in the ice/water interface

Molecular-dynamics simulations of Cl(-) and Na(+) ions are performed to calculate ionic solvation free energies in both bulk simple point-charge/extended water and ice 1 h at several different temperatures, and at the basal ice 1 h/water interface. For the interface we calculate the free energy of "transfer" of the ions across the ice/water interface. For the ions in bulk water in the NPT ensemble at 298 K and 1 atm, results are found to be in good agreement with experiments, and with other simulation results. Simulations performed in the NVT ensemble are shown to give equivalent solvation free energies, and this ensemble is used for the interfacial simulations. Solvation free energies of Cl(-) and Na(+) ions in ice at 150 K are found to be approximately 30 and approximately 20 kcal mol(-1), respectively, less favorable than for water at room temperature. Near the melting point of the model the solvation of the ions in water is the same (within statistical error) as that measured at room temperature, and in the ice is equivalent and approximately 10 kcal mol(-1) less favorable than the liquid. The free energy of transfer for each ion across ice/water interface is calculated and is in good agreement with the bulk observations for the Cl(-) ion. However, for the model of Na(+) the long-range electrostatic contribution to the free energy was more negative in the ice than the liquid, in contrast with the results observed in the bulk calculations.     

Application of the constrained fluid lambda-integration path to the calculation of high temperature Au(110) surface free energies

Recently a method termed constrained fluid lambda-integration was proposed for calculating the free energy difference between bulk solid and liquid reference states via the construction of a reversible thermodynamic integration path; coupling the two states in question. The present work shows how the application of the constrained fluid lambda-integration concept to solid/liquid slab simulation cells makes possible a generally applicable computer simulation methodology for calculating the free energy of any surface and/or surface defect structure, including surfaces requiring variations in surface atom or density number, such as the (1 x 5) Au(100) or (1 x 2) missing row Au(110) reconstructed surfaces or excess adatom/vacancy/step populated surfaces. We evaluate the methodology by calculating the free energy of various disordered high temperature Au(110) embedded atom method surfaces constrained to differing excess surface atom numbers [including those corresponding to the (1 x 2) missing row reconstructed surface] and obtained the interesting result that at 1000 K (as distinct from lower temperatures) the free energy difference between these surfaces is reduced to zero; a result which is consistent with an expected order-disorder phase transition for the Au(110) surface at such high temperatures.     

Transitions between disordered phases in supercooled liquid silicon

We have investigated the transitions between disordered phases in supercooled liquid silicon using computer simulations. The thermodynamic properties were directly obtained from the free energy, which was computed using the recently proposed reversible scaling method. The calculated free energies of the crystalline and liquid phases of silicon at zero pressure, obtained using the environment dependent interatomic potential, are in excellent agreement with the available experimental data. The results show that, at zero pressure, a weak first-order liquid-liquid transition occurs at 1135 K and a continuous liquid-amorphous transition takes place at 843 K. These results are consistent with the existence of a second critical point for the liquid-liquid transition at a negative pressure.     

Hydroxyl radical at the air-water interface

Interaction of the hydroxyl radical with the liquid water surface was studied using classical molecular dynamics computer simulations. From a series of scattering trajectories, the thermal and mass accommodation coefficients of OH on liquid water at 300 K were determined to be 0.95 and 0.83, respectively. The calculated free energy profile for transfer of OH across the air-water interface at 300 K exhibits a minimum in the interfacial region, with the free energy of adsorbtion (DeltaGa) being about 1 kcal/mol more negative than the hydration free energy (DeltaGs). The propensity of the hydroxyl radical for the air-water interface manifests itself in partitioning of OH radicals between the bulk water and the surface. The enhancement of the surface concentration of OH relative to its concentration in the aqueous phase suggests that important OH chemistry may be occurring in the interfacial layer of water droplets, aqueous aerosol particles, and thin water films adsorbed on solid surfaces. This has profound consequences for modeling heterogeneous atmospheric chemical processes.     

Free energy landscapes for homogeneous nucleation of ice for a monatomic water model

We simulate the homogeneous nucleation of ice from supercooled liquid water at 220 K in the isobaric-isothermal ensemble using the MW monatomic water potential. Monte Carlo simulations using umbrella sampling are performed in order to determine the nucleation free energy barrier. We find the Gibbs energy profile to be relatively consistent with that predicted by classical nucleation theory; the free energy barrier to nucleation was determined to be ~18 k(B)T and the critical nucleus comprised ~85 ice particles. Growth from the supercooled liquid gives clusters that are predominantly cubic, whilst starting with a pre-formed subcritical nucleus of cubic or hexagonal ice results in the growth of predominantly that phase of ice only.     

Modeling the cavitation free energy

A new expression to compute the cavitation free energy has been derived by integrating a new model to fit its derivative with respect to the cavity radius. The derivatives were obtained from Monte Carlo simulations data of the contact values of distribution functions for hard-sphere solutes in TIP4P water at 298 K and 1 atm. The new expression, formulated in the framework of the thermodynamics of surfaces and unlike the classical simple models, gives good results also for very small cavities with a radius of approximately 1 A. The contribution to the free energy of a term, which depends on the excess number of molecules at the dividing surface, has been taken into account and discussed for the assumed dependence on r of the surface tension. The asymptotic behavior of the derivative has thus been considered, and a function t(r), which is 0 at r = 0 and 1 at infinity, has been introduced to describe the transition from small to large length regimes. The value of the surface tension obtained by fitting is in very good agreement with that obtained from a simulation of the liquid/vapor interface by using the TIP4P model.     

Adsorption and solvation of ethanol at the water liquid-vapor interface: a molecular dynamics study

The free energy profiles of methanol and ethanol at the water liquid-vapor interface at 310K were calculated using molecular dynamics computer simulations. Both alcohols exhibit a pronounced free energy minimum at the interface and, therefore, have positive adsorption at this interface. The surface excess was computed from the Gibbs adsorption isotherm and was found to be in good agreement with experimental results. Neither compound exhibits a free energy barrier between the bulk and the surface adsorbed state. Scattering calculations of ethanol molecules from a gas phase thermal distribution indicate that the mass accommodation coefficient is 0.98, and the molecules become thermalized within 10 ps of striking the interface. It was determined that the formation of the solvation structure around the ethanol molecule at the interface is not the rate-determining step in its uptake into water droplets. The motion of an ethanol molecule in a water lamella was followed for 30 ns. The time evolution of the probability distribution of finding an ethanol molecule that was initially located at the interface is very well described by the diffusion equation on the free energy surface.     

Thermodynamics of liquid mixtures of krypton+ethane

The total vapour pressure, the excess Gibbs energy Q^E (at 115.77 K), the excess enthalpy H^E (at 117.0 K) and the excess molar volume V^E (at 115.77 K) are reported for liquid mixtures of krypton and ethane. The results are interpreted in the light of some recent statistical theories of liquid mixtures.     

Free energy calculations using flexible-constrained, hard-constrained and non-constrained molecular dynamics simulations.

A comparison of different treatments of bond-stretching interactions in molecular dynamics simulation is presented. Relative free energies from simulations using rigid bonds maintained with the SHAKE algorithm, using partially rigid bonds maintained with a recently introduced flexible constraints algorithm, and using fully flexible bonds are compared in a multi-configurational thermodynamic integration calculation of changing liquid water into liquid methanol. The formula for the free energy change due to a changing flexible constraint in a flexible constraint simulation is derived. To allow for a more direct comparison between these three methods, three different pairs of models for water and methanol were used: a flexible model (simulated without constraints and with flexible constraints), a rigid model (simulated with standard hard constraints), and an alternative flexible model (simulated with flexible constraints and standard hard constraints) in which the ideal or constrained bond lengths correspond to the average bond lengths obtained from a short simulation of the unconstrained flexible model. The particular treatment of the bonds induces differences of up to 2 % in the liquid densities, whereas (excess) free energy differences of up to 5.7 (4.3) kJ mol(-1) are observed. These values are smaller than the differences observed between the three different pairs of methanol/water models: up to 5 % in density and up to 8.5 kJ mol(-1) in (excess) free energy.     

Crystallization kinetics and excess free energy of H2O and D2O nanoscale films of amorphous solid water

Temperature-programmed desorption (TPD) and reflection absorption infrared spectroscopy (RAIRS) are used to investigate the crystallization kinetics and measure the excess free energy of metastable amorphous solid water films (ASW) of H(2)O and D(2)O grown using molecular beams. The desorption rates from the amorphous and crystalline phases of ASW are distinct, and as such, crystallization manifests can be observed in the TPD spectrum. The crystallization kinetics were studied by varying the TPD heating rate from 0.001 to 3 K/s. A coupled desorption-crystallization kinetic model accurately simulates the desorption spectra and accurately predicts the observed temperature shifts in the crystallization. Isothermal crystallization studies using RAIRS are in agreement with the TPD results. Furthermore, highly sensitive measurements of the desorption rates were used to determine the excess free energy of ASW near 150 K. The excess entropy obtained from these data is consistent with there being a thermodynamic continuity between ASW and supercooled liquid water.     

Simulation of adsorption of liquid mixtures of N2 and O2 in a model faujasite cavity at 77.5 K

Grand canonical Monte Carlo simulations of adsorption of N₂ and O₂ and their mixtures in a model zeolitic cavity 14 Å in diameter were performed at 77.5 K for pressures ranging from zero up to saturation, where the adsorbed phase is in equilibrium with coexisting vapor and liquid phases. The same intermolecular potential functions were employed for gas-gas interactions in the vapor, liquid, and adsorbed phases. The gas-solid interaction potential includes dispersion-repulsion energy, induced electrostatic energy, and an ion-quadrupole term to model the interaction of the electric field in zeolites like NaX with polar molecules like N₂. The simulation of the coexisting vapor and liquid phases reproduces the saturation properties of pure liquid oxygen and nitrogen at 77.5 K. Activity coefficients in the adsorbed phase derived from simulations as a function of cavity filling and composition show negative deviations from Raoult's law, even though the non-idealities in the bulk liquid phase have the opposite sign. The simulation of the surface excess isotherm for adsorption from liquid mixtures exhibits preferential adsorption of N₂ and has the commonly-observed quadratic shape skewed toward the more strongly adsorbed component. Micropore condensation is observed for oxygen but not for nitrogen. The condensation of oxygen is similar to a first order phase transition but because of the small number of molecules that can fit into a micropore, coexistence of the two phases is replaced by oscillations between gas- and liquid-like densities.     

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.     

Hydrated excess proton at water-hydrophobic interfaces

The behavior of the hydrated excess proton is investigated at the water-vapor, water-hydrophobic wall, and water-carbon tetrachloride interfaces through molecular dynamics simulations using the third-generation multistate empirical valence bond model (MS-EVB3). The MS-EVB3 simulations show a surface preference of the excess proton at the water-vapor interface, consistent with the discovery of this effect using an earlier version of the MS-EVB model (Petersen et al. J. Phys. Chem. B 2004, 108, 14804) and with the experimental results. The preference of the hydrated excess proton for other water-hydrophobic interfaces is also analyzed for the first time. The hydrated proton structures and charge defect delocalization effects at these interfaces are discussed in detail. By decomposing the free energy profiles into the internal energy and entropic contributions, the thermodynamic (free energy) driving forces for the surface preference of the excess proton are also elaborated. These results indicate that the "rigid" hydrated proton structures at all the interfaces are energetically (as opposed to entropically) stabilized due to the "amphiphilic" nature of the hydrated excess proton, resulting in its overall interfacial concentration enhancement. The effects of acid concentration and nuclear quantization are also explored.     

A first principles molecular dynamics study of excess electron and lithium atom solvation in water-ammonia mixed clusters: structural, spectral, and dynamical behaviors of [(H2O)5NH3]- and Li(H2O)5NH3 at finite temperature

First principles molecular dynamics simulations are carried out to investigate the solvation of an excess electron and a lithium atom in mixed water-ammonia cluster (H(2)O)(5)NH(3) at a finite temperature of 150 K. Both [(H(2)O)(5)NH(3)](-) and Li(H(2)O)(5)NH(3) clusters are seen to display substantial hydrogen bond dynamics due to thermal motion leading to many different isomeric structures. Also, the structures of these two clusters are found to be very different from each other and also very different from the corresponding neutral cluster without any excess electron or the metal atom. Spontaneous ionization of Li atom occurs in the case of Li(H(2)O)(5)NH(3). The spatial distribution of the singly occupied molecular orbital shows where and how the excess (or free) electron is primarily localized in these clusters. The populations of single acceptor (A), double acceptor (AA), and free (NIL) type water and ammonia molecules are found to be significantly high. The dangling hydrogens of these type of water or ammonia molecules are found to primarily capture the free electron. It is also found that the free electron binding motifs evolve with time due to thermal fluctuations and the vertical detachment energy of [(H(2)O)(5)NH(3)](-) and vertical ionization energy of Li(H(2)O)(5)NH(3) also change with time along the simulation trajectories. Assignments of the observed peaks in the vibrational power spectra are done and we found a one to one correlation between the time-averaged populations of water and ammonia molecules at different H-bonding sites with the various peaks of power spectra. The frequency-time correlation functions of OH stretch vibrational frequencies of these clusters are also calculated and their decay profiles are analyzed in terms of the dynamics of hydrogen bonded and dangling OH modes. It is found that the hydrogen bond lifetimes in these clusters are almost five to six times longer than that of pure liquid water at room temperature.