Computation of the free energy of solids

We present a procedure to compute the absolute free energy of solid phases by Monte Carlo simulation. The method is based on the so-called "Einstein-crystal" method of Frenkel and Ladd [J. Chem. Phys. 81, 3188 (1984)]. The new technique is more general and simplifies the calculation for systems with hard core interactions. In addition, the reference Einstein crystal is built up to fulfill translational invariance, which seems to reduce the system size dependence of the results.     

Revisiting the Frenkel-Ladd method to compute the free energy of solids: the Einstein molecule approach

In this paper a new method to evaluate the free energy of solids is proposed. The method can be regarded as a variant of the method proposed by Frenkel and Ladd [J. Chem. Phys. 81, 3188 (1984)]. The main equations of the method can be derived in a simple way. The method can be easily implemented within a Monte Carlo program. We have applied the method to determine the free energy of hard spheres in the solid phase for several system sizes. The obtained free energies agree within the numerical uncertainty with those obtained by Polson et al. [J. Chem. Phys. 112, 5339 (2000)]. The fluid-solid equilibria has been determined for several system sizes and compared to the values published previously by Wilding and Bruce [Phys. Rev. Lett. 85, 5138 (2000)] using the phase switch methodology. It is shown that both the free energies and the coexistence pressures present a strong size dependence and that the results obtained from free energy calculations agree with those obtained using the phase switch method, which constitutes a cross-check of both methodologies. From the results of this work we estimate the coexistence pressure of the fluid-solid transition of hard spheres in the thermodynamic limit to be p*=11.54(4), which is slightly lower than the classical value of Hoover and Ree (p*=11.70) [J. Chem. Phys. 49, 3609 (1968)]. Taking into account the strong size dependence of the free energy of the solid phase, we propose to introduce finite size corrections, which allow us to estimate approximately the free energy of the solid phase in the thermodynamic limit from the known value of the free energy of the solid phase with N molecules. We have also determined the free energy of a Lennard-Jones solid by using both the methodology of this work and the finite size correction. It is shown how a relatively good estimate of the free energy of the system in the thermodynamic limit is obtained even from the free energy of a relatively small system.     

Efficient calculation of temperature dependence of solid-phase free energies by overlap sampling coupled with harmonically targeted perturbation

We examine a method for computing the change in free energy with temperature of a crystalline solid. In the method, the free-energy difference between nearby temperatures is calculated via overlap-sampling free-energy perturbation with Bennett's optimization. Coupled to this is a harmonically targeted perturbation that displaces the atoms in a manner consistent with the temperature change, such that for a harmonic system, the free-energy difference would be recovered with no error. A series of such perturbations can be assembled to bridge larger gaps in temperature. We test this harmonically targeted temperature perturbation (HTTP) method through the application to the inverse-power soft potential, u(r)=ε(σ/r)(n), over a range of temperatures up to the melting condition. Three exponent values (n=12, 9, and 6) for the potential are studied with different crystal structures, specifically face-centered cubic (fcc), body-centered cubic (bcc), and hexagonal close packing. Absolute free energies (classical only) for each system are obtained by implementing the series to near-zero temperature, where the harmonic model becomes very accurate. The HTTP method is shown to provide very precise results, with errors in the free energy smaller than two parts in 10(5). An analysis of the thermodynamic stability of the various structures in the infinite-system limit confirms previous findings. In particular, for n=12 and 9, the fcc structure is stable for all temperatures up to melting, and for n=6, the bcc crystal becomes stable relative to fcc for temperatures above kT/ε=0.802±0.001. The effects of vacancies and other defects are not considered in the analysis.     

Melting points and thermal expansivities of proton-disordered hexagonal ice with several model potentials

A method of free energy calculation is proposed, which enables to cover a wide range of pressure and temperature. The free energies of proton-disordered hexagonal ice (ice Ih) and liquid water are calculated for the TIP4P [J. Chem. Phys. 79, 926 (1983)] model and the TIP5P [J. Chem. Phys. 112, 8910 (2000)] model. From the calculated free energy curves, we determine the melting point of the proton-disordered hexagonal ice at 0.1 MPa (atmospheric pressure), 50 MPa, 100 MPa, and 200 MPa. The melting temperatures at atmospheric pressure for the TIP4P ice and the TIP5P ice are found to be about T(m)=229 K and T(m)=268 K, respectively. The melting temperatures decrease as the pressure is increased, a feature consistent with the pressure dependence of the melting point for realistic proton-disordered hexagonal ice. We also calculate the thermal expansivity of the model ices. Negative thermal expansivity is observed at the low temperature region for the TIP4P ice, but not for the TIP5P ice at the ambient pressure.     

Calculation of solid-liquid interfacial free energy: a classical nucleation theory based approach

We present a simple approach to calculate the solid-liquid interfacial free energy. This new method is based on the classical nucleation theory. Using the molecular dynamics simulation, we employ spherical crystal nuclei embedded in the supercooled liquids to create an ideal model of a homogeneous nucleation. The interfacial free energy is extracted by fitting the relation between the critical nucleus size and the reciprocal of the critical undercooling temperature. The orientationally averaged interfacial free energy is found to be 0.302+/-0.002 (in standard LJ unit). The temperature dependence of the interfacial free energy is also obtained in this work. We find that the interfacial free energy increases slightly with increasing temperature. The positive temperature coefficient of the interfacial free energy is in qualitative agreement with Spaepen's analysis [Solid State Phys. 47, FS181 (1994)] and Turnbull's empirical estimation [J. Appl. Phys. 21, 1022 (1950)].     

Molecular dynamics simulations of surface-initiated melting of nitromethane

The melting of nitromethane initiated at solid-vacuum interfaces has been investigated using molecular dynamics nvt simulations with a realistic force field [D. C. Sorescu et al., J. Phys. Chem. B 104, 8406 (2000)]. The calculated melting point (251+/-5 K) is in good agreement with experiment (244.73 K) and values obtained previously (approximately 255.5 and 266.5+/-8 K) using other simulation methods [P. M. Agrawal et al., J. Chem. Phys. 119, 9617 (2003)]. Analyses of the molecular orientations and diffusion during the simulations as functions of the distance from the exposed surfaces show that the melting is a direct crystal-to-liquid transition, in which the molecules first gain rotational freedom, then mobility. There is a slight dependence of the melting temperature on the exposed crystallographic face.     

Comparison of selected polarizable and nonpolarizable water models in molecular dynamics simulations of ice I(h)

We present a molecular dynamics simulation study in which we determined the melting point of ice I(h) for the polarizable SWM4-NDP water model (Lamoureux et al., Chem. Phys. Lett., 2006, 418, 245-249) and compared the performance of several popular water force fields, both polarizable and nonpolarizable, in terms of melting temperature, stability and orientational structuring of ice. The simulations yield the melting temperature of SWM4-NDP ice as low as T(m) = 185 ± 10 K, despite the quadrupole moment of a SWM4-NDP water molecule being close to the experimental gas phase value. The results thus show that the dependence of T(m) on the molecular quadrupole, observed for the three- and four-site water models, is generally lost if polarization is explicitly included. The study also shows that adding polarizability to a planar three-charge water model increases orientational disorder in hexagonal ice. In addition, analysis of the tetrahedral order in bulk ice reveals a correlation between the pre-existing degree of orientational disorder in ice simulated using different polarizable and nonpolarizable models and the melting temperature of the models. Our findings thus suggest some new considerations regarding the role of polarization forces in a crystalline solid that may guide future development of reliable polarizable water models for ice.     

Molecular dynamics simulations of hexahydro-1,3,5-trinitro-1,3,5-s-triazine (RDX) using a combined Sorescu-Rice-Thompson AMBER force field

We present the results of molecular dynamics simulations of crystalline hexahydro-1,3,5-trinitro-1,3,5-s-triazine (RDX) using the SRT-AMBER force field (P. M. Agrawal et al., J. Phys. Chem. B 2006, 110, 5721), which combines the rigid-molecule force field developed by Sorescu-Rice-Thompson (D. C. Sorescu, B. M. Rice, and D. L. Thompson, J. Phys. Chem. B 1997, 101, 798) with the intramolecular interactions obtained from the Generalized AMBER Force Field (Wang et al., J. Comput. Chem. 2004, 25, 1157). The calculated crystal density at room conditions is about 10% lower than the measured value, while the lattice parameters and thermodynamic melting point are within about 5% at ambient pressure. The chair and inverted chair conformation, bond lengths, and bond angles of the RDX molecule are accurately predicted; however, there are some inaccuracies in the calculated orientations of the NO2 groups. The SRT-AMBER force field predicts overall reasonable results, but modifications, probably in the torsional parameters, are needed for a more accurate force field.     

Molecular dynamics simulations of melting and the glass transition of nitromethane

Molecular dynamics simulations have been used to investigate the thermodynamic melting point of the crystalline nitromethane, the melting mechanism of superheated crystalline nitromethane, and the physical properties of crystalline and glassy nitromethane. The maximum superheating and glass transition temperatures of nitromethane are calculated to be 316 and 160 K, respectively, for heating and cooling rates of 8.9 x 10(9) Ks. Using the hysteresis method [Luo et al., J. Chem. Phys. 120, 11640 (2004)] and by taking the glass transition temperature as the supercooling temperature, we calculate a value of 251.1 K for the thermodynamic melting point, which is in excellent agreement with the two-phase result [Agrawal et al., J. Chem. Phys. 119, 9617 (2003)] of 255.5 K and measured value of 244.73 K. In the melting process, the nitromethane molecules begin to rotate about their lattice positions in the crystal, followed by translational freedom of the molecules. A nucleation mechanism for the melting is illustrated by the distribution of the local translational order parameter. The critical values of the Lindemann index for the C and N atoms immediately prior to melting (the Lindemann criterion) are found to be around 0.155 at 1 atm. The intramolecular motions and molecular structure of nitromethane undergo no abrupt changes upon melting, indicating that the intramolecular degrees of freedom have little effect on the melting. The thermal expansion coefficient and bulk modulus are predicted to be about two or three times larger in crystalline nitromethane than in glassy nitromethane. The vibrational density of states is almost identical in both phases.     

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.     

Direct calculation of the crystal-melt interfacial free energy via molecular dynamics computer simulation

We review our recent work on the direct calculation of the interfacial free energy, gamma, of the crystal-melt interface via molecular dynamics computer simulation for a number of model systems. The value of gamma as a function of crystal orientation is determined using a thermodynamic integration technique employing moving cleaving walls [Phys. Rev. Lett. 2000, 85, 4751]. The calculation is sufficiently precise to resolve the small anisotropy in gamma, which is crucial in determining the kinetics and morphology of dendritic growth. We report values of gamma for the hard-sphere and Lennard-Jones systems, as well as recent results on the series of inverse-power potentials. For the inverse sixth-, seventh-, and eighth-power systems, we determine gamma for both fcc and bcc crystal structures. For these systems, the bcc-melt gamma is lower than that for fcc crystals by about 25%, consistent with recent experiments and computer simulations on fcc-forming systems that show preferential formation of bcc nuclei in the initial stages of crystallization. In addition, we show that our results give a molecular interpretation of Turnbull's rule, which is the empirical relationship between gamma and the enthalpy of fusion.     

Path integral calculation of free energies: quantum effects on the melting temperature of neon

The path integral formulation has been combined with several methods to determine free energies of quantum many-body systems, such as adiabatic switching and reversible scaling. These techniques are alternatives to the standard thermodynamic integration method. A quantum Einstein crystal is used as a model to demonstrate the accuracy and reliability of these free energy methods in quantum simulations. Our main interest focuses on the calculation of the melting temperature of Ne at ambient pressure, taking into account quantum effects in the atomic dynamics. The free energy of the solid was calculated by considering a quantum Einstein crystal as reference state, while for the liquid, the reference state was defined by the classical limit of the fluid. Our findings indicate that, while quantum effects in the melting temperature of this system are small, they still amount to about 6% of the melting temperature, and are therefore not negligible. The particle density as well as the melting enthalpy and entropy of the solid and liquid phases at coexistence is compared to results obtained in the classical limit and also to available experimental data.     

Liquid-solid transition in fully ionized hydrogen at ultra-high pressures

We study the phase diagram of an effective ion model of fully ionized hydrogen at ultra-high pressure. We assume that the protons interact with a screened Coulomb potential derived from a static linear response theory. This model accurately reproduces the physical properties of hydrogen for densities greater than g/ρ(m)=10 cm(3) corresponding to the range of the coupling parameter r(s) ? 0.6. The pressure range, P ? 20 TPa, is well beyond present experimental limitations. Assuming classical protons, we find that the zero temperature enthalpy of the perfect bcc crystal is slightly lower than for other structures at g/ρ(m)=12.47 cm(3) while the fcc structure gains stability at higher density. Using Monte Carlo calculations, we compute the free energy of various phases and locate the melting transition versus density. We find that on melting, bcc is energetically favored with respect to fcc over the entire range investigated. In the solid phase the system undergoes a structural transition from bcc at higher temperature to fcc at lower temperature. The free energy difference between these two structures is very small so that obtaining a quantitative estimate of this second transition line requires accuracy beyond that provided by our method. We estimate the effect of proton zero point motion on the bcc melting line for hydrogen, deuterium, and tritium by a path integral Monte Carlo method. Although zero point effects on hydrogen are large, since the two competing phases (bcc and liquid) have locally similar environments, the effect on the melting line is small; the melting temperature for hydrogen is lowered by about 10% with respect to the classical value.     

Vapor-to-droplet transition in a Lennard-Jones fluid: simulation study of nucleation barriers using the ghost field method

We report a comprehensive Monte Carlo (MC) simulation study of the vapor-to-droplet transition in Lennard-Jones fluid confined to a spherical container with repulsive walls, which is a case study system to investigate homogeneous nucleation. The focus is made on the application of a modified version of the ghost field method (Vishnyakov, A.; Neimark, A. V. J. Chem. Phys. 2003, 119, 9755) to calculate the nucleation barrier. This method allows one to build up a continuous trajectory of equilibrium states stabilized by the ghost field potential, which connects a reference droplet with a reference vapor state. Two computation schemes are employed for free energy calculations, direct thermodynamic integration along the constructed trajectory and umbrella sampling. The nucleation barriers and the size dependence of the surface tension are reported for droplets containing from 260 to 2000 molecules. The MC simulation study is complemented by a review of the simulation methods applied to computing the nucleation barriers and a detailed analysis of the vapor-to-droplet transition by means of the classical nucleation theory.     

Simulation study on the translocation of polymer chains through nanopores

The translocation of polymer chains through nanopores is simulated by dynamical Monte Carlo method. The free energy landscape for the translocation of polymer is calculated by scanning method. The dependence of the free energy barrier Fb and the chemical difference Deltamu on the concentration of chains can explain the behavior of polymer translocation at low and high concentration limits. The relationship between Deltamu and the escaping time tau(2) is in good agreement with the theoretical conclusions obtained by Muthukumar [J. Chem. Phys. 111, 10371 (1999)]. Our simulation results show that the relaxation time is mainly dominated by Fb, while the escaping time is mainly dominated by Deltamu.     

Optimal pairwise and non-pairwise alchemical pathways for free energy calculations of molecular transformation in solution phase

We estimate the global minimum variance path for computing the free energy insertion into or deletion of small molecules from a dense fluid. We perform this optimization over all pair potentials, irrespective of functional form, using functional optimization with a two-body approximation for the radial distribution function. Surprisingly, the optimal pairwise path obtained via this method is almost identical to the path obtained using a optimized generalized "soft core" potential reported by Pham and Shirts [J. Chem. Phys. 135, 034114 (2011)]. We also derive the lowest variance non-pairwise potential path for molecular insertion or deletion and compare its efficiency to the pairwise path. Under certain conditions, non-pairwise pathways can reduce the total variance by up to 60% compared to optimal pairwise pathways. However, optimal non-pairwise pathways do not appear generally feasible for practical free energy calculations because an accurate estimate of the free energy, the parameter that is itself is desired, is required for constructing this non-pairwise path. Additionally, simulations at most intermediate states of these non-pairwise paths have significantly longer correlation times, often exceeding standard simulation lengths for solvation of bulky molecules. The findings suggest that the previously obtained soft core pathway is the lowest variance pathway for molecular insertion or deletion in practice. The findings also demonstrate the utility of functional optimization for determining the efficiency of thermodynamic processes performed with molecular simulation.     

Free volume and water vapor transport properties of temperature‐sensitive polyurethanes

Segmented polyurethanes (PU) with crystalline soft segments were prepared with different crystalline polyols as soft segments. Morphology and microstructure of the PUs were investigated using Differential Scanning Calorimetry (DSC), Wide‐angle X‐ray Diffraction (WAXD), and Positron Annihilation Lifetime Spectra (PALS). Water vapor transport properties of the PU membranes were measured in the temperature range of crystal melting of their soft segments. Dependence of free volume of the PUs on temperature and the relationship between the free volume and water vapor permeability of the PU membranes were investigated. The results show that the mean free volume size and fractional free volume increase more rapidly in the temperature range of crystal melting than in other temperature intervals. In the specified temperature range, water vapor permeability of the polyester based PU membranes increases significantly, caused by the steep increase in free volume, due to crystal melting of the soft segments. Water vapor permeability of the polyester based PUs exhibits approximately direct correlation with the fractional free volume within the temperature range of crystal melting. © 2005 Wiley Periodicals, Inc. J Polym Sci Part B: Polym Phys 43: 1865–1872, 2005     

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.     

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.