Free energy calculations on systems of rigid molecules: An application to the TIP4P model of H2O

Enthalpies and entropies of hydration from Monte Carlo simulations

The changes in free energy, enthalpy, and entropy for transfer of a solute from the gas phase into solution are the fundamental thermodynamic quantities that characterize the solvation process. Owing to the development of methods based on free-energy perturbation theory, computation of free energies of solvation has become routine in conjunction with Monte Carlo (MC) statistical mechanics and molecular dynamics (MD) simulations. Computation of the enthalpy change and by inference the entropy change is more challenging. Two methods are considered in this work corresponding to direct averaging for the solvent and solution and to computing the temperature derivative of the free energy in the van't Hoff approach. The application is for neutral organic solutes in TIP4P water using long MC simulations to improve precision. Definitive results are also provided for pure TIP4P water. While the uncertainty in computed free energies of hydration is ca. 0.05 kcal mol-1, it is ca. 0.4 kcal mol-1 for the enthalpy changes from either van't Hoff plots or the direct method with sampling for 5 billion MC configurations. Partial molar volumes of hydration are also computed by the direct method; they agree well with experimental data with an average deviation of 3 cm3 mol-1. In addition, the results permit breakdown of the errors in the free energy changes from the OPLS-AA force field into their enthalpic and entropic components. The excess hydrophobicity of organic solutes is enthalpic in origin.

Melting Points of OPC and OPC3 Water Models

A recently introduced family of globally optimal water models, OPC, has shown promise in a variety of biomolecular simulations, but properties of these water models outside of the liquid phase remain mostly unexplored. Here, we contribute to filling the gap by reporting melting temperatures of ice I _h of OPC and OPC3 water models. Through the direct coexistence method, which we make available in the AMBER package, the melting points of OPC and OPC3 are estimated as 242 and 210 K, similar to TIP4P-Ew and SPC/E models, respectively, and appreciably below the experimental value of 273.15 K under 1 bar pressure. Water models of the OPC family were optimized to best reproduce water properties in the liquid phase where these models offer noteworthy accuracy advantages over many models of previous generations. It is not surprising that the accuracy of OPC models in describing the phase transition to the solid state does not appear to offer similar improvements. The new anisotropic barostat option implemented in AMBER may benefit system preparation and simulation outside of the direct coexistence applications, such as modeling of membranes or very long DNA strands.

Free Energy of Classical Molecular Crystals by Thermodynamic Integration from a Harmonic Reference

We develop an algorithm for calculating the normal modes of vibration of mechanical systems with constraints, particularly of molecules with rigid bonds and models of rigid molecules, and use it to obtain the harmonic free energy of a crystal. The anharmonic correction is then calculated by the conventional thermodynamic integration over temperature in the NVT ensemble. Attention is paid to finite-size errors, tail corrections, thermostat choice, ergodicity, and other sources of inaccuracies. The calculated free energy of ice XIV modeled by the TIP4P/2005 potential agrees with the previously reported value and is by one order more accurate.

Self-assembly and reconfigurability of shape-shifting particles

Reconfigurability of two-dimensional colloidal crystal structures assembled by anisometric particles capable of changing their shape were studied by molecular dynamics computer simulation. We show that when particles change shape on cue, the assembled structures reconfigure into different ordered structures, structures with improved order, or more densely packed disordered structures, on faster time scales than can be achieved via self-assembly from an initially disordered arrangement. These results suggest that reconfigurable building blocks can be used to assemble reconfigurable materials, as well as to assemble structures not possible otherwise, and that shape shifting could be a promising mechanism to engineer assembly pathways to ordered and disordered structures.

Computing the free energy of molecular solids by the Einstein molecule approach: ices XIII and XIV, hard-dumbbells and a patchy model of proteins

The recently proposed Einstein molecule approach is extended to compute the free energy of molecular solids. This method is a variant of the Einstein crystal method of Frenkel and Ladd [J. Chem. Phys. 81, 3188 (1984)]. In order to show its applicability, we have computed the free energy of a hard-dumbbell solid, of two recently discovered solid phases of water, namely, ice XIII and ice XIV, where the interactions between water molecules are described by the rigid nonpolarizable TIP4P/2005 model potential, and of several solid phases that are thermodynamically stable for an anisotropic patchy model with octahedral symmetry which mimics proteins. Our calculations show that both the Einstein crystal method and the Einstein molecule approach yield the same results within statistical uncertainty. In addition, we have studied in detail some subtle issues concerning the calculation of the free energy of molecular solids. First, for solids with noncubic symmetry, we have studied the effect of the shape of the simulation box on the free energy. Our results show that the equilibrium shape of the simulation box must be used to compute the free energy in order to avoid the appearance of artificial stress in the system that will result in an increase in the free energy. In complex solids, such as the solid phases of water, another difficulty is related to the choice of the reference structure. As in some cases there is no obvious orientation of the molecules; it is not clear how to generate the reference structure. Our results will show that, as long as the structure is not too far from the equilibrium structure, the calculated free energy is invariant to the reference structure used in the free energy calculations. Finally, the strong size dependence of the free energy of solids is also studied.

Vapor-liquid equilibria from the triple point up to the critical point for the new generation of TIP4P-like models: TIP4P/Ew, TIP4P/2005, and TIP4P/ice

The vapor-liquid equilibria of three recently proposed water models have been computed using Gibbs-Duhem simulations. These models are TIP4P/Ew, TIP4P/2005, and TIP4P/ice and can be considered as modified versions of the TIP4P model. By design TIP4P reproduces the vaporization enthalpy of water at room temperature, whereas TIP4P/Ew and TIP4P/2005 match the temperature of maximum density and TIP4P/ice the melting temperature of water. Recently, the melting point for each of these models has been computed, making it possible for the first time to compute the complete vapor-liquid equilibria curve from the triple point to the critical point. From the coexistence results at high temperature, it is possible to estimate the critical properties of these models. None of them is capable of reproducing accurately the critical pressure or the vapor pressures and densities. Additionally, in the cases of TIP4P and TIP4P/ice the critical temperatures are too low and too high, respectively, compared to the experimental value. However, models accounting for the density maximum of water, such as TIP4P/Ew and TIP4P/2005 provide a better estimate of the critical temperature. In particular, TIP4P/2005 provides a critical temperature just 7 K below the experimental result as well as an extraordinarily good description of the liquid densities from the triple point to the critical point. All TIP4P-like models present a ratio of the triple point temperature to the critical point temperature of about 0.39, compared with the experimental value of 0.42. As is the case for any effective potential neglecting many body forces, TIP4P/2005 fails in describing simultaneously the vapor and the liquid phases of water. However, it can be considered as one of the best effective potentials of water for describing condensed phases, both liquid and solid. In fact, it provides a completely coherent view of the phase diagram of water including fluid-solid, solid-solid, and vapor-liquid equilibria.

The melting point of ice Ih for common water models calculated from direct coexistence of the solid-liquid interface

In this work we present an implementation for the calculation of the melting point of ice I(h) from direct coexistence of the solid-liquid interface. We use molecular dynamics simulations of boxes containing liquid water and ice in contact. The implementation is based on the analysis of the evolution of the total energy along NpT simulations at different temperatures. We report the calculation of the melting point of ice I(h) at 1 bar for seven water models: SPC/E, TIP4P, TIP4P-Ew, TIP4P/ice, TIP4P/2005, TIP5P, and TIP5P-E. The results for the melting temperature from the direct coexistence simulations of this work are in agreement (within the statistical uncertainty) with those obtained previously by us from free energy calculations. By taking into account the results of this work and those of our free energy calculations, recommended values of the melting point of ice I(h) at 1 bar for the above mentioned water models are provided.

Free energy of the solid C60 fullerene orientational order-disorder transition

The free energies of the orientationally ordered crystal phase of C60 at low temperatures and the disordered crystal phase at high temperatures are calculated to an accuracy of +/-0.05 kJ/mol using the expanded ensemble Monte Carlo method with the potential model of Sprik et al. [J. Phys. Chem. 96, 2027 (1992)]. The order-disorder transition temperature at zero pressure is determined directly from these free energies, and is found to be consistent with the abrupt changes in configurational energy and unit cell size also found in simulation. A modification of the potential results in predictions of the transition temperature of 257 K and the entropy change of 18.1 J/mol K at this transition, which are in good agreement with the experimental values of 260 K and 19 J/mol K, respectively. The orientational distinguishability in the ordered phase and the indistinguishability in the disordered phase lead to a contribution to the entropy difference of k ln 60, with 60 being the symmetry number of C60. This quantum mechanical correction is important for the accurate prediction of the phase transition properties of the C60 crystals.

Growth Mechanism of a Gas Clathrate Hydrate from a Dilute Aqueous Gas Solution: A Molecular Dynamics Simulation of a Three-Phase System

A molecular dynamics simulation of a three-phase system including a gas clathrate, liquid water, and a gas was carried out at 298 K and high pressure in order to investigate the growth mechanism of the clathrate from a dilute aqueous gas solution. The simulation indicated that the clathrate grew on interfaces between the clathrate and the liquid water, after transfer of the gas molecules from the gas phase to the interfaces. The results suggest a two-step process for growth: first, gas molecules are arranged at cage sites, and second, H(2)O molecules are ordered near the gas molecules. The results also suggest that only the H(2)O molecules, which are surrounded or sandwiched by the gas molecules, form the stable polygons that constitute the cages of the clathrate. In addition, the growth of the clathrate from a concentrated aqueous gas solution was also simulated, and the results suggested a growth mechanism in which many H(2)O and gas molecules correctively form the structure of the clathrate. The clathrate grown from the concentrated solution contained some empty cages, whereas the formation of empty cages was not observed during the growth from the dilute solution. The results obtained by both simulations are compared with the results of an experimental study, and the growth mechanism of the clathrate in a real system is discussed.

Entropy-based measure of structural order in water

We analyze the nature of the structural order established in liquid TIP4P water in the framework provided by the multiparticle correlation expansion of the statistical entropy. Different regimes are mapped onto the phase diagram of the model upon resolving the pair entropy into its translational and orientational components. These parameters are used to quantify the relative amounts of positional and angular order in a given thermodynamic state, thus allowing a structurally unbiased definition of low-density and high-density water. As a result, the structurally anomalous region--within which both types of order are simultaneously disrupted by an increase of pressure at constant temperature--is clearly identified through extensive molecular-dynamics simulations.

Hybrid integral equation/simulation model for enhancing free energy computations

Integral equation theory is used for extrapolating free energy data from molecular simulations of a reference state with respect to a modification of the interaction potential. The methodology is applied to the correction of artefacts arising from potential shifting and truncation. Corrective contributions for the hydration free energy with respect to the full potential are analysed for the case that both the solute-solvent as well as the solvent-solvent potentials are truncated and modified by a shifted-force term, reaching beyond the range of the dielectric continuum approximation and simple long-range correction expressions. The model systems argon in water and pure water are used as examples for apolar and polar solutes, revealing significant correction contributions even for the short-ranged dispersive interactions and the magnitude of solute-solvent and solvent-solvent components. In comparison with simulation-based extrapolation techniques the integral equation method is shown to be capable of quantitatively predicting truncation artefacts at negligible computational overhead.

Reproducible Polypeptide Folding and Structure Prediction using Molecular Dynamics Simulations

The folding of a polypeptide from an extended state to a well-defined conformation is studied using microsecond classical molecular dynamics (MD) simulations and replica exchange molecular dynamics (REMD) simulations in explicit solvent and in vacuo. It is shown that the solvated peptide folds many times in the REMD simulations but only a few times in the conventional simulations. From the folding events in the classical simulations we estimate an approximate folding time of 1-2 micros. The REMD simulations allow enough sampling to deduce a detailed Gibbs free energy landscape in three dimensions. The global minimum of the energy landscape corresponds to the native state of the peptide as determined previously by nuclear magnetic resonance (NMR) experiments. Starting from an extended state it takes about 50 ns before the native structure appears in the REMD simulations, about an order of magnitude faster than conventional MD. The calculated melting curve is in good qualitative agreement with experiment. In vacuo, the peptide collapses rapidly to a conformation that is substantially different from the native state in solvent.

Relation between the melting temperature and the temperature of maximum density for the most common models of water

Water exhibits a maximum in density at normal pressure at 4 degrees above its melting point. The reproduction of this maximum is a stringent test for potential models used commonly in simulations of water. The relation between the melting temperature and the temperature of maximum density for these potential models is unknown mainly due to our ignorance about the melting temperature of these models. Recently we have determined the melting temperature of ice I(h) for several commonly used models of water (SPC, SPC/E, TIP3P, TIP4P, TIP4P/Ew, and TIP5P). In this work we locate the temperature of maximum density for these models. In this way the relative location of the temperature of maximum density with respect to the melting temperature is established. For SPC, SPC/E, TIP3P, TIP4P, and TIP4P/Ew the maximum in density occurs at about 21-37 K above the melting temperature. In all these models the negative charge is located either on the oxygen itself or on a point along the H-O-H bisector. For the TIP5P and TIP5P-E models the maximum in density occurs at about 11 K above the melting temperature. The location of the negative charge appears as a geometrical crucial factor to the relative position of the temperature of maximum density with respect to the melting temperature.

Computer simulation of the 13 crystalline phases of ice

As a reference for follow-up studies toward more accurate model parametrizations, we performed molecular-dynamics and Monte Carlo simulations for all known crystalline phases of ice, as described by the simple point-charge/extended and TIP4P water models. We started from the measured structures, densities, and temperatures, and carried out classical canonical simulations for all these arrangements. All simulated samples were cooled down close to 0 K to facilitate the comparison with theoretical estimates. We determined configurational internal energies as well as pressures, and monitored how accurately the measured configurations were preserved during the simulations. While these two models predicted very similar thermophysical and structural properties for water at ambient conditions, the predicted features for the corresponding ice polymorphs may differ significantly.

A potential model for the study of ices and amorphous water: TIP4P/Ice

The ability of several water models to predict the properties of ices is discussed. The emphasis is put on the results for the densities and the coexistence curves between the different ice forms. It is concluded that none of the most commonly used rigid models is satisfactory. A new model specifically designed to cope with solid-phase properties is proposed. The parameters have been obtained by fitting the equation of state and selected points of the melting lines and of the coexistence lines involving different ice forms. The phase diagram is then calculated for the new potential. The predicted melting temperature of hexagonal ice (Ih) at 1 bar is 272.2 K. This excellent value does not imply a deterioration of the rest of the properties. In fact, the predictions for both the densities and the coexistence curves are better than for TIP4P, which previously yielded the best estimations of the ice properties.

The melting temperature of the most common models of water

The melting temperature of ice I(h) for several commonly used models of water (SPC, SPC/E,TIP3P,TIP4P, TIP4P/Ew, and TIP5P) is obtained from computer simulations at p = 1 bar. Since the melting temperature of ice I(h) for the TIP4P model is now known [E. Sanz, C. Vega, J. L. F. Abascal, and L. G. MacDowell, Phys. Rev. Lett. 92, 255701 (2004)], it is possible to use the Gibbs-Duhem methodology [D. Kofke, J. Chem. Phys. 98, 4149 (1993)] to evaluate the melting temperature of ice I(h) for other potential models of water. We have found that the melting temperatures of ice I(h) for SPC, SPC/E, TIP3P, TIP4P, TIP4P/Ew, and TIP5P models are T = 190 K, 215 K, 146 K, 232 K, 245 K, and 274 K, respectively. The relative stability of ice I(h) with respect to ice II for these models has also been considered. It turns out that for SPC, SPC/E, TIP3P, and TIP5P the stable phase at the normal melting point is ice II (so that ice I(h) is not a thermodynamically stable phase for these models). For TIP4P and TIP4P/Ew, ice I(h) is the stable solid phase at the standard melting point. The location of the negative charge along the H-O-H bisector appears as a critical factor in the determination of the relative stability between the I(h) and II ice forms. The methodology proposed in this paper can be used to investigate the effect upon a coexistence line due to a change in the potential parameters.

Simulations of proton order and disorder in ice Ih

Computer simulations of ice Ih with different proton orientations are presented. Simulations of proton disordered ice are carried out using a Monte Carlo method which samples over proton degree of freedom, allowing for the calculation of the dielectric constant and for the examination of the degree of proton disorder. Simulations are also presented for two proton ordered structures of ice Ih, the ferroelectric Cmc2(1) structure or ice XI and the antiferroelectric Pna2(1) structure. These simulations indicate that a transition to a proton ordered phase occurs at low temperatures (below 80 K). The symmetry of the ordered phase is found to be dependent on the water potential. The stability of the two proton ordered structures is due to a balance of short-ranged interactions which tend to stabilize the Pna2(1) structure and longer-range interactions which stabilize the Cmc2(1) structure.

Combinatorial entropy and phase diagram of partially ordered ice phases

A close analytical estimate for the combinatorial entropy of partially ordered ice phases is presented. The expression obtained is very general, as it can be used for any ice phase obeying the Bernal-Fowler rules. The only input required is a number of crystallographic parameters, and the experimentally observed proton site occupancies. For fully disordered phases such as hexagonal ice, it recovers the result deduced by Pauling, while for fully ordered ice it is found to vanish. Although the space groups determined for ice I, VI, and VII require random proton site occupancies, it is found that such random allocation of protons does not necessarily imply random orientational disorder. The theoretical estimate for the combinatorial entropy is employed together with free energy calculations in order to obtain the phase diagram of ice from 0 to 10 GPa. Overall qualitative agreement with experiment is found for the TIP4P model of water. An accurate estimate of the combinatorial entropy is found to play an important role in determining the stability of partially ordered ice phases, such as ice III and ice V.

Temperature dependence of the hydrophobic hydration and interaction of simple solutes: An examination of five popular water models

We examine five different popular rigid water models (SPC, SPCE, TIP3P, TIP4P, and TIP5P) using molecular dynamics simulations in order to investigate the hydrophobic hydration and interaction of apolar Lennard-Jones solutes as a function of temperature in the range between 275 and 375 K along the 0.1 MPa isobar. For all investigated models and state points we calculate the excess chemical potential for the noble gases and methane employing the Widom particle insertion technique. All water models exhibit too small hydration entropies, but show a clear hierarchy. TIP3P shows poorest agreement with experiment, whereas TIP5P is closest to the experimental data at lower temperatures and SPCE is closest at higher temperatures. As a first approximation, this behavior can be rationalized as a temperature shift with respect to the solvation behavior found in real water. A rescaling procedure inspired by the information theory model of Hummer et al. [Chem. Phys. 258, 349 (2000)] suggests that the different solubility curves for the different models and real water can be largely explained on the basis of the different density curves at constant pressure. In addition, the models that give a good representation of the water structure at ambient conditions (TIP5P, SPCE, and TIP4P) show considerably better agreement with the experimental data than the ones which exhibit less structured O-O correlation functions (SPC and TIP3P). In the second part of the paper we calculate the hydrophobic interaction between xenon particles directly from a series of 60 ns simulation runs. We find that the temperature dependence of the association is to a large extent related to the strength of the solvation entropy. Nevertheless, differences between the models seem to require a more detailed molecular picture. The TIP5P model shows by far the strongest temperature dependence. The suggested density rescaling is also applied to the chemical potential in the xenon-xenon contact-pair configuration, indicating the presence of a temperature where the hydrophobic interaction turns into purely repulsive. The predicted association for xenon in real water suggests the presence of a strong variation with temperature, comparable to the behavior found for TIP5P water. Comparing different water models and experimental data we conclude that a proper description of density effects is an important requirement for a water model to account correctly for the correct description of the hydrophobic effects. A water model exhibiting a density maximum at the correct temperature is desirable.

Sticking of CO to crystalline and amorphous ice surfaces

We present results of classical trajectory calculations on the sticking of hyperthermal CO to the basal plane (0001) face of crystalline ice Ih and to the surface of amorphous ice Ia. The calculations were performed for normal incidence at a surface temperature Ts = 90 K for ice Ia, and at Ts = 90 and 150 K for ice Ih. For both surfaces, the sticking probability can be fitted to a simple exponentially decaying function of the incidence energy, Ei: Ps = 1.0e(-Ei(kJ/mol)/90(kJ/mol)) at Ts = 90 K. The energy transfer from the impinging molecule to the crystalline and the amorphous surface is found to be quite efficient, in agreement with the results of molecular beam experiments on the scattering of the similar molecule, N2, from crystalline and amorphous ice. However, the energy transfer is less efficient for amorphous than for crystalline ice. Our calculations predict that the sticking probability decreases with Ts for CO scattering from crystalline ice, as the energy transfer from the impinging molecule to the warmer surfaces becomes less efficient. At high Ei (up to 193 kJ/mol), no surface penetration occurs in the case of crystalline ice. However, for CO colliding with the amorphous surface, a penetrating trajectory was observed to occur into a large water pore. The molecular dynamics calculations predict that the average potential energy of CO adsorbed to ice Ih is -10.1 +/- 0.2 and -8.4 +/- 0.2 kJ/mol for CO adsorbed to ice Ia. These values are in agreement with previous experimental and theoretical data. The distribution of the potential energy of CO adsorbed to ice Ia was found to be wider (with a standard deviation sigma of 2.4 kJ/mol) than that of CO interacting with ice Ih (sigma = 2.0 kJ/mol). In collisions with ice Ia, the CO molecules scatter at larger angles and over a wider distribution of angles than in collisions with ice Ih.

Water in nanopores. I. Coexistence curves from Gibbs ensemble Monte Carlo simulations

Coexistence curves of water in cylindrical and slitlike nanopores of different size and water-substrate interaction strength were simulated in the Gibbs ensemble. The two-phase coexistence regions cover a wide range of pore filling level and temperature, including ambient temperature. Five different kinds of two-phase coexistence are observed. A single liquid-vapor coexistence is observed in hydrophobic and moderately hydrophilic pores. Surface transitions split from the main liquid-vapor coexistence region, when the water-substrate interaction becomes comparable or stronger than the water-water pair interaction. In this case prewetting, one and two layering transitions were observed. The critical temperature of the first layering transition decreases with strengthening water-substrate interaction towards the critical temperature expected for two-dimensional systems and is not sensitive to the variation of pore size and shape. Liquid-vapor phase transition in a pore with a wall which is already covered with two water layers is most typical for hydrophilic pores. The critical temperature of this transition is very sensitive to the pore size, in contrast to the liquid-vapor critical temperature in hydrophobic pores. The observed rich phase behavior of water in pores evidences that the knowledge of coexistence curves is of crucial importance for the analysis of experimental results and a prerequiste of meaningful simulations.