Understanding free-energy perturbation calculations through a model of harmonic oscillators: theory and implications to improve the sampling efficiency by molecular simulation

Free-energy perturbation calculation is frequently used to calculate free-energy differences because it is easy to implement and the computation is fast. However, the calculation is subject to large inaccuracies in some circumstances due to the insufficient sampling of the relevant tails of the energy-difference distributions. Here we expand this knowledge of insufficient sampling into a two-dimensional (2D) energy space using a model of harmonic oscillators. We show analytically the relation between the energies of the sampling system and those of the desired target energy spaces, which provide the basis to understand the difficulties in free-energy perturbation calculations. We clarify the reasons of the inaccurate calculation in the different harmonic cases that stem from the spatial separations of the reference and the target energy pairs located in the two-dimensional energy space. The potential-energy space introduced into this 2D energy-space model provides additional clues to improve the sampling efficiency. Based on this understanding, we propose two ways to calculate the free-energy differences using the two schemes of the distribution method. We show that the distribution method implemented in the appropriate energy space--the energy-difference space and the potential-energy space, respectively--can improve the calculation of free energies in different circumstances. This analysis implies that the sampling can be improved if it is directed toward the appropriate region in the potential-energy space, which is easily implemented in various types of free-energy calculations. To test this, we calculate the free-energy surface of alanine dipeptide in gas phase and in aqueous phase, respectively. We demonstrate that the free-energy surface calculation is improved when the biased sampling of the potential energy is integrated into the sampling scheme.     

Rosenbluth-sampled nonequilibrium work method for calculation of free energies in molecular simulation

We present methods that introduce concepts from Rosenbluth sampling [M. N. Rosenbluth and A. W. Rosenbluth, J. Chem. Phys. 23, 356 (1955)] into the Jarzynski nonequilibrium work (NEW) free-energy calculation technique [C. Jarzynski, Phys. Rev. Lett. 78, 2690 (1997)]. The proposed hybrid modifies the way steps are taken in the NEW process. With it, each step is selected from a range of alternatives, with bias given to steps that contribute the least work. The definition of the work average is modified to account for the bias. We introduce two variants of this method, lambda-bias sampling and configuration-bias sampling, respectively; a combined lambda- and configuration-bias method is also considered. By reducing the likelihood that large nonequilibrated work values enter the ensemble average, the Rosenbluth sampling aids in remedying problems of inaccuracy of the calculation. We demonstrate the performance of the proposed methods through a model system of N independent harmonic oscillators. This model captures the difficulties involved in calculating free energies in real systems while retaining many tractable features that are helpful to the study. We examine four variants of this model that differ qualitatively in the nature of their phase-space overlap. Results indicate that the lambda-bias sampling method is most useful for systems with entropic sampling barriers, while the configuration-bias methods are best for systems with energetic sampling barriers. The Rosenbluth-sampling schemes yield much more accurate results than the unbiased nonequilibrium work method. Typically the accuracy can be improved by about an order of magnitude for a given amount of sampling; this improvement translates into two or more orders of magnitude less sampling required to obtain a given level of accuracy, owing to the generally slow convergence of the NEW calculation when the inaccuracy is large.     

Semiclassical calculation of nonadiabatic thermal rate constants: application to condensed phase reactions

The nonadiabatic transition state theory proposed recently by Zhao et al. [J. Chem. Phys. 121, 8854 (2004)] is extended to calculate rate constants of complex systems by using the Monte Carlo and umbrella sampling methods. Surface hopping molecular dynamics technique is incorporated to take into account the dynamic recrossing effect. A nontrivial benchmark model of the nonadiabatic reaction in the condensed phase is used for the numerical test. It is found that our semiclassical results agree well with those produced by the rigorous quantum mechanical method. Comparing with available analytical approaches, we find that the simple statistical theory proposed by Straub and Berne [J. Chem. Phys. 87, 6111 (1987)] is applicable for a wide friction region although their formula is obtained using Landau-Zener [Phys. Z. Sowjetunion 2, 46 (1932); Proc. R. Soc. London, Ser. A 137, 696 (1932)] nonadiabatic transition probability along a one-dimensional diffusive coordinate. We also investigate how the nuclear tunneling events affect the dependence of the rate constant on the friction.     

Folding kinetics of a lattice protein via a forward flux sampling approach

We implement a forward flux sampling approach [R. J. Allen et al., J. Chem. Phys. 124, 194111 (2006)] for calculating transition rate constants and for sampling paths of protein folding events. The algorithm generates trajectories for the transition between the unfolded and folded states as chains of partially connected paths, which can be used to obtain the transition-state ensemble and the properties that characterize these intermediates. We apply this approach to Monte Carlo simulations of a model lattice protein in open space and in confined spaces of varying dimensions. We study the effect of confinement on both protein thermodynamic stability and folding kinetics; the former by mapping free-energy landscapes and the latter by the determination of rate constants and mechanistic details of the folding pathway. Our results show that, for the range of temperatures where the native state is stable, confinement of a protein destabilizes the unfolded state by reducing its entropy, resulting in increased thermodynamic stability of the folded state. Relative to the folding in open space, we find that the kinetics can be accelerated at temperatures above the temperature at which the unconfined protein folds fastest and that the rate constant increases with the number of constrained dimensions. By examining the statistical properties of the transition-state ensemble, we detect signs of a classical nucleation folding mechanism for a core of native contacts formed at an early stage of the process. This nucleus acts as folding foci and is composed of those residues that have higher probability to form native contacts in the transition-state intermediates, which can vary depending on the confinement conditions of the system.     

Thermodynamical and structural properties of imidazolium based ionic liquids from molecular simulation

We have performed molecular dynamics simulations to determine the densities and heat of vaporization as well as structural information for the 1-alkyl-3-methyl-imidazolium based ionic liquids [amim][Cl] and [amim][BF(4)] in the temperature range from 298 to 363 K. In this simulation study, we used an united atom model of Liu et al. [Phys. Chem. Chem. Phys. 8, 1096 (2006)] for the [emim(+)] and [bmim(+)] cations, which we have extended for simulation in [hmim]-ILs and combined with parameters of Canongia Lopes et al. [J. Phys. Chem. B 108, 2038 (2004)] for the [Cl(-)] anion. Our simulation results prove that both the original united atoms approach by Liu et al. and our extension yield reasonable predictions for the ionic liquid with a considerably reduced computational expense than that required for all atoms models. Radial distribution functions and spatial distribution functions where employed to analyze the local structure of this ionic liquid, and in which way it is influenced by the type of the anion, the size of the cation, and the temperature. Our simulations give evidence for the occurrence of tail aggregations in these ionic liquids with increasing length of the side chain and also increasing temperature.     

Computer simulation of the role of torsional flexibility on mass and momentum transport for a series of linear alkanes

We present the results obtained from a systematic equilibrium molecular dynamics study of the effect of torsional flexibility on the diffusion and viscosity of a series of linear alkanes. To make unambiguous comparisons between molecules with torsional flexibility and those without, we use the frozen distribution sampling (FDS) method introduced by Travis et al. [J. Chem. Phys. 98, 1524 (1993); J. Chem. Phys. 102, 2174 (1995)] but modified and updated for increased efficiency. We first demonstrate comprehensively that FDS guarantees corresponding thermodynamic states. We then show that removal of torsional flexibility results in a significant lowering of the diffusion coefficient (and corresponding increase in shear viscosity) and furthermore that this effect increases with increasing chain length. The results are discussed in terms of the possible mechanism giving rise to this dynamic coupling phenomenon.     

Crystal nucleation of hard spheres using molecular dynamics, umbrella sampling, and forward flux sampling: a comparison of simulation techniques

Over the last number of years several simulation methods have been introduced to study rare events such as nucleation. In this paper we examine the crystal nucleation rate of hard spheres using three such numerical techniques: molecular dynamics, forward flux sampling, and a Bennett-Chandler-type theory where the nucleation barrier is determined using umbrella sampling simulations. The resulting nucleation rates are compared with the experimental rates of Harland and van Megen [Phys. Rev. E 55, 3054 (1997)], Sinn et al. [Prog. Colloid Polym. Sci. 118, 266 (2001)], Sch?tzel and Ackerson [Phys. Rev. E 48, 3766 (1993)], and the predicted rates for monodisperse and 5% polydisperse hard spheres of Auer and Frenkel [Nature 409, 1020 (2001)]. When the rates are examined in units of the long-time diffusion coefficient, we find agreement between all the theoretically predicted nucleation rates, however, the experimental results display a markedly different behavior for low supersaturation. Additionally, we examined the precritical nuclei arising in the molecular dynamics, forward flux sampling, and umbrella sampling simulations. The structure of the nuclei appears independent of the simulation method, and in all cases, the nuclei contains on average significantly more face-centered-cubic ordered particles than hexagonal-close-packed ordered particles.     

Ordering of limits in the Jarzynski equality

We consider the sampling problems encountered in computing free-energy differences using Jarzynski's nonequilibrium work relation [Phys. Rev. Lett. 56, 2690 (1997)]. This relation expresses the free-energy change of a system, on which finite-time work is done, as an average over all possible trajectories of the system. This average can then be expressed as a cumulant expansion of the work. We study the scaling of these cumulants with an appropriately defined measure of phase-space accessibility epsilon and particle number N for two simple changes in state. We find that the cumulant expansion is slowly convergent for predominantly entropic processes, those where epsilon is considerably altered during the course of the process. An accurate determination of the free-energy change requires some knowledge of the behavior of the tails of the work distribution associated with the process. Jarzynski's irreversible work relation is only valid with the correct ordering of the infinite limits of N and epsilon, clarifying the regime of its applicability.     

Melting temperature of ice Ih calculated from coexisting solid-liquid phases

We carried out molecular-dynamics simulations by using the two-phase coexistence method with the constant pressure, particle number, and enthalpy ensemble to compute the melting temperature of proton-disordered hexagonal ice I(h) at 1-bar pressure. Four models of water were considered, including the widely used TIP4P [W. L. Jorgensen, J. Chandrasekha, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem. Phys.79, 926 (1983)] and TIP5P [M. W. Mahoney and W. L. Jorgensen J. Chem. Phys.112, 8910 (2000)] models, as well as recently improved TIP4P and TIP5P models for use with Ewald techniques-the TIP4P-Ew [W. Horn, W. C. Swope, J. W. Pitera, J. C. Madura, T. J. Dick, G. L. Hura, and T. Head-Gordon, J. Chem. Phys.120, 9665 (2004)] and TIP5P-Ew [S. W. Rick, J. Chem. Phys.120, 6085 (2004)] models. The calculated melting temperature at 1 bar is T(m) = 229 +/- 1 K for the TIP4P and T(m) = 272.0 +/- 0.6 K for the TIP5P ice I(h), both are consistent with previous simulations based on free-energy methods. For the TIP4P-Ew and TIP5P-Ew models, the calculated melting temperature is T(m) = 257.0 +/- 1.1 K and T(m) = 253.9 +/- 1.1 K, respectively.     

Bond-order potentials with split-charge equilibration: application to C-, H-, and O-containing systems

A method for extending charge transfer to bond-order potentials, known as the bond-order potential/split-charge equilibration (BOP/SQE) method [P. T. Mikulski, M. T. Knippenberg, and J. A. Harrison, J. Chem. Phys. 131, 241105 (2009)], is integrated into a new bond-order potential for interactions between oxygen, carbon, and hydrogen. This reactive potential utilizes the formalism of the adaptive intermolecular reactive empirical bond-order potential [S. J. Stuart, A. B. Tutein, and J. A. Harrison, J. Chem. Phys. 112, 6472 (2000)] with additional terms for oxygen and charge interactions. This implementation of the reactive potential is able to model chemical reactions where partial charges change in gas- and condensed-phase systems containing oxygen, carbon, and hydrogen. The BOP/SQE method prevents the unrestricted growth of charges, often observed in charge equilibration methods, without adding significant computational time, because it makes use of a quantity which is calculated as part of the underlying covalent portion of the potential, namely, the bond order. The implementation of this method with the qAIREBO potential is designed to provide a tool that can be used to model dynamics in a wide range of systems without significant computational cost. To demonstrate the usefulness and flexibility of this potential, heats of formation for isolated molecules, radial distribution functions of liquids, and energies of oxygenated diamond surfaces are calculated.     

Computing accurate potentials of mean force in electrolyte solutions with the generalized gradient-augmented harmonic Fourier beads method

We establish the accuracy of the novel generalized gradient-augmented harmonic Fourier beads (ggaHFB) method in computing free-energy profiles or potentials of mean force (PMFs) through comparison with two independent conventional techniques. In particular, we employ umbrella sampling with one dimensional weighted histogram analysis method (WHAM) and free molecular dynamics simulation of radial distribution functions to compute the PMF for the Na(+)-Cl(-) ion-pair separation to 16 A in 1.0M NaCl solution in water. The corresponding ggaHFB free-energy profile in six dimensional Cartesian space is in excellent agreement with the conventional benchmarks. We then explore changes in the PMF in response to lowering the NaCl concentration to physiological 0.3 and 0.1M, and dilute 0.0M concentrations. Finally, to expand the scope of the ggaHFB method, we formally develop the free-energy gradient approximation in arbitrary nonlinear coordinates. This formal development underscores the importance of the logarithmic Jacobian correction to reconstruct true PMFs from umbrella sampling simulations with either WHAM or ggaHFB techniques when nonlinear coordinate restraints are used with Cartesian propagators. The ability to employ nonlinear coordinates and high accuracy of the computed free-energy profiles further advocate the use of the ggaHFB method in studies of rare events in complex systems.     

Long-range-corrected hybrids using a range-separated Perdew-Burke-Ernzerhof functional and random phase approximation correlation

We build on methods combining a short-range density functional approximation with a long-range random phase approximation [B. G. Janesko, T. M. Henderson, and G. E. Scuseria, J. Chem. Phys. 130, 081105 (2009)] or second-order screened exchange [J. Paier et al., J. Chem. Phys. 132, 094103 (2010)] by replacing the range-separated local density approximation functional with a range-separated generalized gradient approximation functional in the short range. We present benchmark results that show a marked improvement in the thermodynamic tests over the previous local density approximation-based methods while retaining those methods' excellent performance in van der Waals interactions.     

Rapid free energy calculation of peptide self-assembly by REMD umbrella sampling

We extend umbrella sampling with replica exchange steps to calculate free energies that are important in the self-assembly of peptides. This leads to a more than 10-fold speed up over conventional umbrella sampling, thereby providing a practical method to calculate these free-energy differences. This approach can also observe first-order phase transitions and pinpoint the location of the concomitant boundary. When conformational changes are involved, this method can handle peptides up to approximately 7 residues, providing a rapid and accurate assessment of the thermodynamic properties of model systems, and can thus be used to answer fundamental questions about peptide self-assembly. When no major conformational changes are involved, we expect the size limit to be equal to that of standard molecular dynamics.     

Rovibrational energy transfer in ortho-H2+para-H2 collisions

We present the results of a full-dimensional quantum mechanical study of the rovibrational energy transfer in the collision between ortho-H2 and para-H2 in the energy range of 0.1-1.0 eV. The multiconfiguration time-dependent Hartree algorithm has been used to propagate the wave packets on the global potential energy surface by Boothroyd et al. [J. Chem. Phys. 116, 666 (2002)] and on a modification of this surface where the short range anisotropy is reduced. State-to-state attributes such as probabilities or integral cross sections are obtained using the formalism of Tannor and Weeks [J. Chem. Phys. 98, 3884 (1993)] by Fourier transforming the correlation functions. The effect of initial rotation of the diatoms on the inelastic and de-excitation processes is investigated.     

Accurate and efficient method for the treatment of exchange in a plane-wave basis

We describe an accurate and efficient extension of Chawla and Voth's [J. Chem. Phys. 108, 4697 (1998)] plane-wave based algorithm for calculating exchange energies, exchange energy densities, and exchange energy gradients with respect to wave-function parameters in systems of electrons subject to periodic boundary conditions. The theory and numerical results show that the computational effort scales almost linearly with the number of plane waves and quadratically with the number of k vectors. To obtain high accuracy with relatively few k vectors, we use an adaptation of Gygi and Baldereschi's [Phys. Rev. B 34, 4405 (1986)] method for reducing Brillouin-zone integration errors.     

Theoretical scheme for the shear viscosity of Lennard-Jones fluids

We use the shear viscosity expression from the Enskog theory of dense gases in a perturbative scheme for the Lennard-Jones (LJ) fluid. This perturbative scheme is formulated by combining the analytic rational function approximation method of Bravo Yuste and Santos [Phys. Rev. A 43, 5418 (1991)] for the radial distribution function of hard-sphere fluids and the well known Mansoori-Canfield/Rasaiah-Stell perturbation theory to determine an effective diameter for the LJ fluid. The scheme is reliable on a wide range of temperatures and densities, and is very accurate around the critical point. Using this information, we build an accurate empirical formula for the shear viscosity in the liquid phase, which fits the recent data [K. Meier et al., J. Chem. Phys. 121, 3671 (2004)] in the whole simulation range.     

State-to-state rotational transitions in H2+H2 collisions at low temperatures

We present quantum mechanical close-coupling calculations of collisions between two hydrogen molecules over a wide range of energies, extending from the ultracold limit to the superthermal region. The two most recently published potential energy surfaces for the H(2)-H(2) complex, the so-called Diep-Johnson (DJ) [J. Chem. Phys. 112, 4465 (2000); 113, 3480 (2000)] and Boothroyd-Martin-Keogh-Peterson (BMKP) [J. Chem. Phys. 116, 666 (2002)] surfaces, are quantitatively evaluated and compared through the investigation of rotational transitions in H(2)+H(2) collisions within rigid rotor approximation. The BMKP surface is expected to be an improvement, approaching chemical accuracy, over all conformations of the potential energy surface compared to previous calculations of H(2)-H(2) interaction. We found significant differences in rotational excitation/deexcitation cross sections computed on the two surfaces in collisions between two para-H(2) molecules. The discrepancy persists over a large range of energies from the ultracold regime to thermal energies and occurs for several low-lying initial rotational levels. Good agreement is found with experiment B. Mate et al., [J. Chem. Phys. 122, 064313 (2005)] for the lowest rotational excitation process, but only with the use of the DJ potential. Rate coefficients computed with the BMKP potential are an order of magnitude smaller.     

Long-time and unitary properties of semiclassical initial value representations

We numerically compare the semiclassical "frozen Gaussian" Herman-Kluk propagator [Chem. Phys. 91, 27 (1984)] and the "thawed Gaussian" propagator put forward recently by Baranger et al. [J. Phys. A 34, 7227 (2001)] by studying the quantum dynamics in some nonlinear one-dimensional potentials. The reasons for the lack of long-time accuracy and norm conservation in the latter method are uncovered. We amend the thawed Gaussian propagator with a global harmonic approximation for the stability of the trajectories and demonstrate that this revised propagator is a true alternative to the Herman-Kluk propagator with similar accuracy.