Understanding Molecular Dynamics with Stochastic Processes via Real or Virtual Dynamics

Molecular dynamics with the stochastic process provides a convenient way to compute structural and thermodynamic properties of chemical, biological, and materials systems. It is demonstrated that the virtual dynamics case that we proposed for the Langevin equation[J. Chem. Phys. 147 , 184104 (2017)] in principle exists in other types of stochastic thermostats as well. The recommended "middle" scheme[J. Chem. Phys. 147 , 034109 (2017)] of the Andersen thermostat is investigated as an example. As shown by both analytic and numerical results, while the real and virtual dynamics cases approach the same plateau of the characteristic correlation time in the high collision frequency limit, the accuracy and efficiency of sampling are relatively insensitive to the value of the collision frequency in a broad range. After we compare the behaviors of the Andersen thermostat to those of Langevin dynamics, a heuristic schematic representation is proposed for understanding efficient stochastic thermostatting processes with molecular dynamics.     

A path-integral Langevin equation treatment of low-temperature doped helium clusters

We present an implementation of path integral molecular dynamics for sampling low temperature properties of doped helium clusters using Langevin dynamics. The robustness of the path integral Langevin equation and white-noise Langevin equation [M. Ceriotti, M. Parrinello, T. E. Markland, and D. E. Manolopoulos, J. Chem. Phys. 133, 124104 (2010)] sampling methods are considered for those weakly bound systems with comparison to path integral Monte Carlo (PIMC) in terms of efficiency and accuracy. Using these techniques, convergence studies are performed to confirm the systematic error reduction introduced by increasing the number of discretization steps of the path integral. We comment on the structural and energetic evolution of He(N)-CO(2) clusters from N = 1 to 20. To quantify the importance of both rotations and exchange in our simulations, we present a chemical potential and calculated band origin shifts as a function of cluster size utilizing PIMC sampling that includes these effects. This work also serves to showcase the implementation of path integral simulation techniques within the molecular modelling toolkit [K. Hinsen, J. Comp. Chem. 21, 79 (2000)], an open-source molecular simulation package.     

An idealized model for nonequilibrium dynamics in molecular systems

The nonequilibrium dynamics of highly nonlinear and multidimensional systems can give rise to emergent chemical behavior which can often be tracked using low-dimensional order parameters such as a reaction path. Such behavior cannot be readily surmised by stationary projected stochastic representations such as those described by the Langevin equation or the generalized Langevin equation (GLE). The irreversible generalized Langevin equation (iGLE) contains a nonstationary friction kernel that in certain limits reduces to the GLE with space-dependent friction. For more general forms of the friction kernel, the iGLE was previously shown to be the projection of a mechanical system with a time-dependent Hamiltonian [R. Hernandez, J. Chem. Phys. 110, 7701 (1999)]. In the present work, the corresponding open Hamiltonian system is shown to be amenable to numerical integration despite the presence of a nonlocal term. Simulations of this mechanical system further confirm that the time dependence of the observed total energy and the correlations of the solvent force are in precise agreement with the projected iGLE. This extended nonstationary Hamiltonian is thus amenable to the study of nonequilibrium bounds and fluctuation theorems.     

Implementations of Nosé-Hoover and Nosé-Poincaré thermostats in mesoscopic dynamic simulations with the united-residue model of a polypeptide chain

Molecular dynamics (MD) simulations generate a canonical ensemble only when integration of the equations of motion is coupled to a thermostat. Three extended phase space thermostats, one version of Nose-Hoover and two versions of Nose-Poincare, are compared with each other and with the Berendsen thermostat and Langevin stochastic dynamics. Implementation of extended phase space thermostats was first tested on a model Lennard-Jones fluid system; subsequently, they were implemented with our physics-based protein united-residue (UNRES) force field MD. The thermostats were also implemented and tested for the multiple-time-step reversible reference system propagator (RESPA). The velocity and temperature distributions were analyzed to confirm that the proper canonical distribution is generated by each simulation. The value of the artificial mass constant, Q, of the thermostat has a large influence on the distribution of the temperatures sampled during UNRES simulations (the velocity distributions were affected only slightly). The numerical stabilities of all three algorithms were compared with each other and with that of microcanonical MD. Both Nose-Poincare thermostats, which are symplectic, were not very stable for both the Lennard-Jones fluid and UNRES MD simulations started from nonequilibrated structures which implies major changes of the potential energy throughout a trajectory. Even though the Nose-Hoover thermostat does not have a canonical symplectic structure, it is the most stable algorithm for UNRES MD simulations. For UNRES with RESPA, the "extended system inside-reference system propagator algorithm" of the RESPA implementation of the Nose-Hoover thermostat was the only stable algorithm, and enabled us to increase the integration time step.     

Activation processes with memory

We propose a mathematical treatment of the activated processes governed by stochastic Langevin dynamics with a colored random force, corresponding to a noise generated by an Ornstein-Uhlenbeck process. Such non-Markovian dynamics take place in a variety of chemical and biological systems. Using the path integral approach, we constructed the conditional probability for passing between two stationary states in configurational space. Our relations can be used for Monte Carlo sampling of evolution trajectories for systems with many degrees of freedom as well as for determining the reaction coordinate used in transition state theory. On the basis of our relation for a conditional probability, we generalize the method of determining the most probable path to the case of colored random force. Using the simple three-hole potential, we examine numerically the effect of nonzero correlation time (memory) on the evolution of the most probable path for a finite temperature.     

Accelerating the convergence of path integral dynamics with a generalized Langevin equation

The quantum nature of nuclei plays an important role in the accurate modelling of light atoms such as hydrogen, but it is often neglected in simulations due to the high computational overhead involved. It has recently been shown that zero-point energy effects can be included comparatively cheaply in simulations of harmonic and quasiharmonic systems by augmenting classical molecular dynamics with a generalized Langevin equation (GLE). Here we describe how a similar approach can be used to accelerate the convergence of path integral (PI) molecular dynamics to the exact quantum mechanical result in more strongly anharmonic systems exhibiting both zero point energy and tunnelling effects. The resulting PI-GLE method is illustrated with applications to a double-well tunnelling problem and to liquid water.     

New observations regarding deterministic, time-reversible thermostats and Gauss's principle of least constraint

Deterministic thermostats are frequently employed in nonequilibrium molecular dynamics simulations in order to remove the heat produced irreversibly over the course of such simulations. The simplest thermostat is the Gaussian thermostat, which satisfies Gauss's principle of least constraint and fixes the peculiar kinetic energy. There are of course infinitely many ways to thermostat systems, e.g., by fixing sigma(i)/p(i)/mu+l. In the present paper we provide, for the first time, convincing arguments as to why the conventional Gaussian isokinetic thermostat (mu = 1) is unique in this class. We show that this thermostat minimizes the phase space compression and is the only thermostat for which the conjugate pairing rule holds. Moreover, it is shown that for finite sized systems in the absence of an applied dissipative field, all other thermostats (mu not = 1) perform work on the system in the same manner as a dissipative field while simultaneously removing the dissipative heat so generated. All other thermostats (mu not = 1) are thus autodissipative. Among all mu, thermostats, only the mu = 1 Gaussian thermostat permits an equilibrium state.     

Conserving the linear momentum in stochastic dynamics: Dissipative particle dynamics as a general strategy to achieve local thermostatization in molecular dynamics simulations

Stochastic dynamics is a widely employed strategy to achieve local thermostatization in molecular dynamics simulation studies; however, it suffers from an inherent violation of momentum conservation. Although this short‐coming has little impact on structural and short‐time dynamic properties, it can be shown that dynamics in the long‐time limit such as diffusion is strongly dependent on the respective thermostat setting. Application of the methodically similar dissipative particle dynamics (DPD) provides a simple, effective strategy to ensure the advantages of local, stochastic thermostatization while at the same time the linear momentum of the system remains conserved. In this work, the key parameters to employ the DPD thermostats in the framework of periodic boundary conditions are investigated, in particular the dependence of the system properties on the size of the DPD‐region as well as the treatment of forces near the cutoff. Structural and dynamical data for light and heavy water as well as a Lennard–Jones fluid have been compared to simulations executed via stochastic dynamics as well as via use of the widely employed Nose–Hoover chain and Berendsen thermostats. It is demonstrated that a small size of the DPD region is sufficient to achieve local thermalization, while at the same time artifacts in the self‐diffusion characteristic for stochastic dynamics are eliminated. © 2016 Wiley Periodicals, Inc.     

Quantum generalized Langevin equation approach to gas/solid collisions

The classical generalized Langevin equation (GLE) approach to gas/solid collisions is generalized to quantum scattering. Using Feynman's method of partial path integration, the full gas/solid propagator is reduced to a form in which only the dynamics of the incident atom and the surface oscillator(s) directly struck appear explicitly. Solving this effective dynamical problem in the semiclassical limit yields a stationary phase equation of motion identical in form to the classical GLE. The noise, however, is distributed according to quantum rather than classical statistics. From the GLE a quantum phase can be constructed and an S-matrix computed. The resulting theory is capable of describing inelastic-diffractive scattering which has been seen experimentally by Williams.     

Synchronization of trajectories in canonical molecular-dynamics simulations: observation, explanation, and exploitation

For two methods commonly used to achieve canonical-ensemble sampling in a molecular-dynamics simulation, the Langevin thermostat and the Andersen [H. C. Andersen, J. Chem. Phys. 72, 2384 (1980)] thermostat, we observe, as have others, synchronization of initially independent trajectories in the same potential basin when the same random number sequence is employed. For the first time, we derive the time dependence of this synchronization for a harmonic well and show that the rate of synchronization is proportional to the thermostat coupling strength at weak coupling and inversely proportional at strong coupling with a peak in between. Explanations for the synchronization and the coupling dependence are given for both thermostats. Observation of the effect for a realistic 97-atom system indicates that this phenomenon is quite general. We discuss some of the implications of this effect and propose that it can be exploited to develop new simulation techniques. We give three examples: efficient thermalization (a concept which was also noted by Fahy and Hamann [S. Fahy and D. R. Hamann, Phys. Rev. Lett. 69, 761 (1992)]), time-parallelization of a trajectory in an infrequent-event system, and detecting transitions in an infrequent-event system.     

Why and how do systems react in thermally fluctuating environments?

Many chemical reactions, including those of biological importance, take place in thermally fluctuating environments. Compared to isolated systems, there arise markedly different features due to the effects of energy dissipation through friction and stochastic driving by random forces reflecting the fluctuation of the environment. Investigation of how robustly the system reacts under the influence of thermal fluctuation, and elucidating the role of thermal fluctuation in the reaction are significant subjects in the study of chemical reactions. In this article, we start with overviewing the generalized Langevin equation (GLE), which has long been used and continues to be a powerful tool to describe a system surrounded by a thermal environment. It has been also generalized further to treat a nonstationary environment, in which the conventional fluctuation-dissipation theorem no longer holds. Then, within the framework of the Langevin equation we present a method recently developed to extract a new reaction coordinate that is decoupled from all the other coordinates in the region of a rank-one saddle linking the reactant and the product. The reaction coordinate is buried in nonlinear couplings among the original coordinates under the influence of stochastic random force. It was ensured that the sign of this new reaction coordinate (= a nonlinear functional of the original coordinates, velocities, friction, and random force) at any instant is sufficient to determine in which region, the reactant or the product, the system finally arrives. We also discuss how one can extend the method to extract such a coordinate from the GLE framework in stationary and nonstationary environments, where memory effects exist in dynamics of the reaction.     

Direct assessment of quantum nuclear effects on hydrogen bond strength by constrained-centroid ab initio path integral molecular dynamics

The impact of quantum nuclear effects on hydrogen (H-) bond strength has been inferred in earlier work from bond lengths obtained from path integral molecular dynamics (PIMD) simulations. To obtain a direct quantitative assessment of such effects, we use constrained-centroid PIMD simulations to calculate the free energy changes upon breaking the H-bonds in dimers of HF and water. Comparing ab initio simulations performed using PIMD and classical nucleus molecular dynamics (MD), we find smaller dissociation free energies with the PIMD method. Specifically, at 50 K, the H-bond in (HF)(2) is about 30% weaker when quantum nuclear effects are included, while that in (H(2)O)(2) is about 15% weaker. In a complementary set of simulations, we compare unconstrained PIMD and classical nucleus MD simulations to assess the influence of quantum nuclei on the structures of these systems. We find increased heavy atom distances, indicating weakening of the H-bond consistent with that observed by direct calculation of the free energies of dissociation.     

Approximate first passage time distribution for barrier crossing in a double well under fractional Gaussian noise

The distribution of waiting times, f(t), between successive turnovers in the catalytic action of single molecules of the enzyme beta-galactosidase has recently been determined in closed form by Chaudhury and Cherayil [J. Chem. Phys. 125, 024904 (2006)] using a one-dimensional generalized Langevin equation (GLE) formalism in combination with Kramers' flux-over-population approach to barrier crossing dynamics. The present paper provides an alternative derivation of f(t) that eschews this approach, which is strictly applicable only under conditions of local equilibrium. In this alternative derivation, a double well potential is incorporated into the GLE, along with a colored noise term representing protein conformational fluctuations, and the resulting equation transformed approximately to a Smoluchowski-type equation. f(t) is identified with the first passage time distribution for a particle to reach the barrier top starting from an equilibrium distribution of initial points, and is determined from the solution of the above equation using local boundary conditions. The use of such boundary conditions is necessitated by the absence of definite information about the precise nature of the boundary conditions applicable to stochastic processes governed by non-Markovian dynamics. f(t) calculated in this way is found to have the same analytic structure as the distribution calculated by the flux-over-population method.     

Effective thermostat induced by coarse graining of simple point charge water

We investigate how the transport properties of a united atom fluid with a dissipative particle dynamics thermostat depend on the functional form and magnitude of both the conservative and the stochastic interactions. We demonstrate how the thermostat strongly affects the hydrodynamics, especially diffusion, viscosity, and local escape times. As model system we use simple point charge (SPC) water, from which projected trajectories are used to determine the effective interactions in the united atom model. The simulation results support our argument that the thermostat should be viewed as an integral part of the coarse-grained dynamics rather than a tool for approaching thermal equilibrium. As our main result we show that the united atom model with the adjusted effective interactions approximately reproduces the diffusion constant and the viscosity of the underlying detailed SPC water model.     

Basis set study of classical rotor lattice dynamics

The reorientational relaxation of molecular systems is important in many phenomenon and applications. In this paper, we explore the reorientational relaxation of a model Brownian rotor lattice system with short range interactions in both the high and low temperature regimes. In this study, we use a basis set expansion to capture collective motions of the system. The single particle basis set is used in the high temperature regime, while the spin wave basis is used in the low temperature regime. The equations of motion derived in this approach are analogous to the generalized Langevin equation, but the equations render flexibility by allowing nonequilibrium initial conditions. This calculation shows that the choice of projection operators in the generalized Langevin equation (GLE) approach corresponds to defining a specific inner-product space, and this inner-product space should be chosen to reveal the important physics of the problem. The basis set approach corresponds to an inner-product and projection operator that maintain the orthogonality of the spherical harmonics and provide a convenient platform for analyzing GLE expansions. The results compare favorably with numerical simulations, and the formalism is easily extended to more complex systems.     

On the short-time limit of ring polymer molecular dynamics

We examine the short-time accuracy of a class of approximate quantum dynamical techniques that includes the centroid molecular dynamics (CMD) and ring polymer molecular dynamics (RPMD) methods. Both of these methods are based on the path integral molecular dynamics (PIMD) technique for calculating the exact static equilibrium properties of quantum mechanical systems. For Kubo-transformed real-time correlation functions involving operators that are linear functions of positions or momenta, the RPMD and (adiabatic) CMD approximations differ only in the choice of the artificial mass matrix of the system of ring polymer beads that is employed in PIMD. The obvious ansatz for a general method of this type is therefore to regard the elements of the PIMD (or Parrinello-Rahman) mass matrix as an adjustable set of parameters that can be chosen to improve the accuracy of the resulting approximation. We show here that this ansatz leads uniquely to the RPMD approximation when the criterion that is used to select the mass matrix is the short-time accuracy of the Kubo-transformed correlation function. In particular, we show that the leading error in the RPMD position autocorrelation function is O(t(8)) and the error in the velocity autocorrelation function is O(t(6)), for a general anharmonic potential. The corresponding errors in the CMD approximation are O(t(6)) and O(t(4)), respectively.     

Improved configuration space sampling: Langevin dynamics with alternative mobility

We present a new and efficient method for determining optimal configurations of a large number (N) of interacting particles. We use a coarse-grained stochastic Langevin equation in the overdamped limit to describe the dynamics of this system and replace the standard mobility by an effective space dependent inverse Hessian correlation matrix. Due to the analogy of the drift term in the Langevin equation and the update scheme in Newton's method, we expect accelerated dynamics or improved convergence in the convex part of the potential energy surface Phi. The stochastic noise term, however, is not only essential for proper thermodynamic sampling but also allows the system to access transition states in the concave parts of Phi. We employ a Broyden-Fletcher-Goldfarb-Shannon method for updating the local mobility matrix. Quantitative analysis for one and two dimensional systems shows that the new method is indeed more efficient than standard methods with constant effective friction. Due to the construction, our effective mobility adapts high values/low friction in configurations which are less optimal and low values/high friction in configurations that are more optimal.     

The Langevin Hull: Constant pressure and temperature dynamics for non-periodic systems

We have developed a new isobaric-isothermal (NPT) algorithm which applies an external pressure to the facets comprising the convex hull surrounding the system. A Langevin thermostat is also applied to the facets to mimic contact with an external heat bath. This new method, the "Langevin Hull", can handle heterogeneous mixtures of materials with different compressibilities. These systems are problematic for traditional affine transform methods. The Langevin Hull does not suffer from the edge effects of boundary potential methods, and allows realistic treatment of both external pressure and thermal conductivity due to the presence of an implicit solvent. We apply this method to several different systems including bare metal nanoparticles, nanoparticles in an explicit solvent, as well as clusters of liquid water. The predicted mechanical properties of these systems are in good agreement with experimental data and previous simulation work.