Search

Determining collective variables (CVs) for conformational transitions is crucial to understanding their dynamics and targeting them in enhanced sampling simulations. Often, CVs are proposed based on intuition or prior knowledge of a system. However, the problem of systematically determining a proper reaction coordinate (RC) for a specific process in terms of a set of putative CVs can be achieved using committor analysis (CA). Identifying essential degrees of freedom that govern such transitions using CA remains elusive because of the high dimensionality of the conformational space. Various schemes exist to leverage the power of machine learning (ML) to extract an RC from CA. Here, we extend these studies and compare the ability of 17 different ML schemes to identify accurate RCs associated with conformational transitions. We tested these methods on an alanine dipeptide in vacuum and on a sarcosine dipeptoid in an implicit solvent. Our comparison revealed that the light gradient boosting machine method outperforms other methods. In order to extract key features from the models, we employed Shapley Additive exPlanations analysis and compared its interpretation with the "feature importance" approach. For the alanine dipeptide, our methodology identifies ϕ and θ dihedrals as essential degrees of freedom in the C7ax to C7eq transition. For the sarcosine dipeptoid system, the dihedrals ψ and ω are the most important for the cisαD to transαD transition. We further argue that analysis of the full dynamical pathway, and not just endpoint states, is essential for identifying key degrees of freedom governing transitions. [ABSTRACT FROM AUTHOR]

Quantum time correlation functions (TCFs) involving two states are important for describing nonadiabatic dynamical processes such as charge transfer (CT). Based on a previous single-state method, we propose an imaginary-time open-chain path-integral (OCPI) approach for evaluating the two-state symmetrized TCFs. Expressing the forward and backward propagation on different electronic potential energy surfaces as a complex-time path integral, we then transform the path variables to average and difference variables such that the integration over the difference variables up to the second order can be performed analytically. The resulting expression for the symmetrized TCF is equivalent to sampling the open-chain configurations in an effective potential that corresponds to the average surface. Using importance sampling over the extended OCPI space via open path-integral molecular dynamics, we tested the resulting path-integral approximation by calculating the Fermi's golden rule CT rate constant within a widely used spin-boson model. Comparing with the real-time linearized semiclassical method and analytical result, we show that the imaginary-time OCPI provides an accurate two-state symmetrized TCF and rate constant in the typical turnover region. It is shown that the first bead of the open chain corresponds to physical zero-time and that the endpoint bead corresponds to final time t; oscillations of the end-to-end distance perfectly match the nuclear mode frequency. The two-state OCPI scheme is seen to capture the tested model's electronic quantum coherence and nuclear quantum effects accurately. [ABSTRACT FROM AUTHOR]

A method for calculating the free energy difference between two structurally defined conformational states of a chemical system is developed. A path is defined using a previously reported collective variable that interpolates between two or more conformations, and a restraint is introduced in order to keep the system close to the path. The evolution of the system along the path, which typically presents a high free energy barrier, is generated using enhanced sampling schemes. Although the formulation of the method in terms of a path is quite general, an important advance in this work is the demonstration that prior knowledge of the path is, in fact, not needed and that the free energy difference can be obtained using a simplified definition of the path collective variable that only involves the endpoints. We first validate this method on cyclohexane isomerization. The method is then tested for an extensive conformational change in a realistic molecular system by calculating the free energy difference between the α-helix and β-hairpin conformations of deca-alanine in solution. Finally, the method is applied to a biologically relevant system to calculate the free energy difference of an observed and a hypothetical conformation of an antigenic peptide bound to a major histocompatibility complex. [ABSTRACT FROM AUTHOR]

Enhanced sampling techniques that target a set of collective variables and that use molecular dynamics as the driving engine have seen widespread application in the computational molecular sciences as a means to explore the free-energy landscapes of complex systems. The use of molecular dynamics as the fundamental driver of the sampling requires the introduction of a time step whose magnitude is limited by the fastest motions in a system. While standard multiple time-stepping methods allow larger time steps to be employed for the slower and computationally more expensive forces, the maximum achievable increase in time step is limited by resonance phenomena, which inextricably couple fast and slow motions. Recently, we introduced deterministic and stochastic resonance-free multiple time step algorithms for molecular dynamics that solve this resonance problem and allow ten- to twenty-fold gains in the large time step compared to standard multiple time step algorithms [P. Minary et al., Phys. Rev. Lett. 93, 150201 (2004); B. Leimkuhler et al., Mol. Phys. 111, 3579-3594 (2013)]. These methods are based on the imposition of isokinetic constraints that couple the physical system to Nosé-Hoover chains or Nos'e-Hoover Langevin schemes. In this paper, we show how to adapt these methods for collective variable-based enhanced sampling techniques, specifically adiabatic free-energy dynamics/temperature-accelerated molecular dynamics, unified free-energy dynamics, and by extension, metadynamics, thus allowing simulations employing these methods to employ similarly very large time steps. The combination of resonance-free multiple time step integrators with free-energy-based enhanced sampling significantly improves the efficiency of conformational exploration. [ABSTRACT FROM AUTHOR]

Path integral-based methodologies play a crucial role for the investigation of nuclear quantum effects by means of computer simulations. However, these techniques are significantly more demanding than corresponding classical simulations. To reduce this numerical effort, we recently proposed a method, based on a rigorous Hamiltonian formulation, which restricts the quantum modeling to a small but relevant spatial region within a larger reservoir where particles are treated classically. In this work, we extend this idea and show how it can be implemented along with state-of-the-art path integral simulation techniques, including path-integral molecular dynamics, which allows for the calculation of quantum statistical properties, and ring-polymer and centroid molecular dynamics, which allow the calculation of approximate quantum dynamical properties. To this end, we derive a new integration algorithm that also makes use of multiple time-stepping. The scheme is validated via adaptive classical-path-integral simulations of liquid water. Potential applications of the proposed multiresolution method are diverse and include efficient quantum simulations of interfaces as well as complex biomolecular systems such as membranes and proteins. [ABSTRACT FROM AUTHOR]

Higher order factorization schemes are developed for path integral molecular dynamics in order to improve the convergence of estimators for physical observables as a function of the Trotter number. The methods are based on the Takahashi-Imada and Susuki decompositions of the Boltzmann operator. The methods introduced improve the averages of the estimators by using the classical forces needed to carry out the dynamics to construct a posteriori weighting factors for standard path integral molecular dynamics. The new approaches are straightforward to implement in existing path integral codes and carry no significant overhead. The Suzuki higher order factorization was also used to improve the end-to-end distance estimator in open path integral molecular dynamics. The new schemes are tested in various model systems, including an ab initio path integral molecular dynamics calculation on the hydrogen molecule and a quantum water model. The proposed algorithms have potential utility for reducing the cost of path integral molecular dynamics calculations of bulk systems. [ABSTRACT FROM AUTHOR]

Extraction of relaxation times, lifetimes, and rates associated with the transport of topological charge defects in hydrogen-bonded networks from molecular dynamics simulations is a challenge because proton transfer reactions continually change the identity of the defect core. In this paper, we present a statistical mechanical theory that allows these quantities to be computed in an unbiased manner. The theory employs a set of suitably defined indicator or population functions for locating a defect structure and their associated correlation functions. These functions are then used to develop a chemical master equation framework from which the rates and lifetimes can be determined. Furthermore, we develop an integral equation formalism for connecting various types of population correlation functions and derive an iterative solution to the equation, which is given a graphical interpretation. The chemical master equation framework is applied to the problems of both hydronium and hydroxide transport in bulk water. For each case it is shown that the theory establishes direct links between the defect's dominant solvation structures, the kinetics of charge transfer, and the mechanism of structural diffusion. A detailed analysis is presented for aqueous hydroxide, examining both reorientational time scales and relaxation of the rotational anisotropy, which is correlated with recent experimental results for these quantities. Finally, for OH-(aq) it is demonstrated that the 'dynamical hypercoordination mechanism' is consistent with available experimental data while other mechanistic proposals are shown to fail. As a means of going beyond the linear rate theory valid from short up to intermediate time scales, a fractional kinetic model is introduced in the Appendix in order to describe the nonexponential long-time behavior of time-correlation functions. Within the mathematical framework of fractional calculus the power law decay ∼t-σ, where σ is a parameter of the model and depends on the dimensionality of the system, is obtained from Mittag-Leffler functions due to their long-time asymptotics, whereas (stretched) exponential behavior is found for short times. [ABSTRACT FROM AUTHOR]

The problems of ergodicity and internal consistency in the centroid and ring-polymer molecular dynamics methods are addressed in the context of a comparative study of the two methods. Enhanced sampling in ring-polymer molecular dynamics (RPMD) is achieved by first performing an equilibrium path integral calculation and then launching RPMD trajectories from selected, stochastically independent equilibrium configurations. It is shown that this approach converges more rapidly than periodic resampling of velocities from a single long RPMD run. Dynamical quantities obtained from RPMD and centroid molecular dynamics (CMD) are compared to exact results for a variety of model systems. Fully converged results for correlations functions are presented for several one dimensional systems and para-hydrogen near its triple point using an improved sampling technique. Our results indicate that CMD shows very similar performance to RPMD. The quality of each method is further assessed via a new χ2 descriptor constructed by transforming approximate real-time correlation functions from CMD and RPMD trajectories to imaginary time and comparing these to numerically exact imaginary time correlation functions. For para-hydrogen near its triple point, it is found that adiabatic CMD and RPMD both have similar χ2 error. [ABSTRACT FROM AUTHOR]

An efficient computational approach to perform Car–Parrinello ab initio molecular dynamics (CPAIMD) simulations under cluster (free) boundary conditions is presented. The general approach builds upon a recent real-space CPAIMD formalism using discrete variable representation (DVR) basis sets [Y. Liu et al., Phys. Rev. B 12, 125110 (2003); H.-S. Lee and M. E. Tuckerman, J. Phys. Chem. A 110, 5549 (2006)]. In order to satisfy cluster boundary conditions, a DVR based on sinc functions is utilized to expand the Kohn–Sham orbitals and electron density. Poisson’s equation is solved in order to calculate the Hartree potential via an integral representation of the 1/r singularity. Excellent convergence properties are achieved with respect to the number of grid points (or DVR functions) and the size of the simulation cell. A straightforward implementation of the present approach leads to near linear scaling [O(N4/3)] of the computational cost with respect to the system size (N) for the solution of Poisson’s equation. The accuracy and stability of CPAIMD simulations based on sinc DVR are tested for a model problem as well as for N2 and a water dimer. [ABSTRACT FROM AUTHOR]

Dynamical properties of liquid water were studied using Car-Parrinello [Phys. Rev. Lett. 55, 2471 (1985)] ab initio molecular dynamics (AIMD) simulations within the Kohn-Sham (KS) density functional theory employing the Becke-Lee-Yang-Parr exchange-correlation functional for the electronic structure. The KS orbitals were expanded in a discrete variable representation basis set, wherein the complete basis set limit can be easily reached and which, therefore, provides complete convergence of ionic forces. In order to minimize possible nonergodic behavior of the simulated water system in a constant energy (NVE) ensemble, a long equilibration run (30 ps) preceded a 60 ps long production run. The temperature drift during the entire 60 ps trajectory was found to be minimal. The diffusion coefficient [0.055 Å2/ps] obtained from the present work for 32 D2O molecules is a factor of 4 smaller than the most up to date experimental value, but significantly larger than those of other recent AIMD studies. Adjusting the experimental result so as to match the finite-sized system used in the present study brings the comparison between theory and experiment to within a factor of 3. More importantly, the system is not observed to become “glassy” as has been reported in previous AIMD studies. The computed infrared spectrum is in good agreement with experimental data, especially in the low frequency regime where the translational and librational motions of water are manifested. The long simulation length also made it possible to perform detailed studies of hydrogen bond dynamics. The relaxation dynamics of hydrogen bonds observed in the present AIMD simulation is slower than those of popular force fields, such as the TIP4P potential, but comparable to that of the TIP5P potential. [ABSTRACT FROM AUTHOR]

We present a rigorous description of chemical space within a molecular grand-canonical ensemble multi-component density functional theory framework. A total energy density functional for chemical compounds in contact with an electron and a proton bath is introduced using Lagrange multipliers which correspond to the energetic response to changes of the elementary particle densities. From a generalized Gibbs-Duhem equation analog, reactivity indices such as the nuclear hardness and a molecular Fukui function, which couples the grand-canonical electronic and nuclear degrees of freedom, are obtained. Maxwell relations between composition particles, ionic displacements, and the external potential are discussed. Numerical results for the molecular Fukui function are presented as well as finite temperature estimates for the oxidation of ammonia. [ABSTRACT FROM AUTHOR]

Structural properties of liquid water at ambient temperature were studied using Car-Parrinello [Phys. Rev. Lett. 55, 2471 (1985)] ab initio molecular dynamics (CPAIMD) simulations combined with the Kohn-Sham (KS) density functional theory and the BLYP exchange-correlation functional for the electronic structure. Unlike other recent work on the same subject, where plane-wave (PW) or hybrid Gaussian/plane-wave basis sets were employed, in the present paper, a discrete variable representation (DVR) basis set is used to expand the KS orbitals, so that with the real-space grid adapted in the present work, the properties of liquid water could be obtained very near the complete basis set limit. Structural properties of liquid water were extracted from a 30 ps CPAIMD-BLYP/DVR trajectory at 300 K. The radial distribution functions (RDFs), spatial distribution functions, and hydrogen bond geometry obtained from the CPAIMD-BLYP/DVR simulation are generally in good agreement with the most up to date experimental measurements. Compared to recent ab initio MD simulations based on PW basis sets, less significant overstructuring was found in the RDFs and the distributions of hydrogen bond angles, suggesting that previous plane-wave and Gaussian basis set calculations have exaggerated the tendency toward overstructuring. [ABSTRACT FROM AUTHOR]

A new molecular dynamics method for calculating free energies associated with transformations of the thermodynamic state or chemical composition of a system (also known as alchemical transformations) is presented. The new method extends the adiabatic dynamics approach recently introduced by Rosso et al. [J. Chem. Phys. 116, 4389 (2002)] and is based on the use of an additional degree of freedom, λ, that is used as a switching parameter between the potential energy functions that characterize the two states. In the new method, the coupling parameter λ is introduced as a fictitious dynamical variable in the Hamiltonian, and a system of switching functions is employed that leads to a barrier in the λ free energy profile between the relevant thermodynamic end points. The presence of such a barrier, therefore, enhances sampling in the end point (λ=0 and λ=1) regions which are most important for computing relevant free energy differences. In order to ensure efficient barrier crossing, a high temperature Tλ is assigned to λ and a fictitious mass mλ is introduced as a means of creating an adiabatic separation between λ and the rest of the system. Under these conditions, it is shown that the λ free energy profile can be directly computed from the adiabatic probability distribution function of λ without any postprocessing or unbiasing of the output data. The new method is illustrated on two model problems and in the calculation of the solvation free energy of amino acid side-chain analogs in TIP3P water. Comparisons to previous work using thermodynamic integration and free energy perturbation show that the new λ adiabatic free energy dynamics method results in very precise free energy calculations using significantly shorter trajectories. [ABSTRACT FROM AUTHOR]

A multistate empirical valence bond model for proton transport in water, which explicitly includes solvent polarization, is presented. Polarization is included for each valence-bond state via induced point dipoles, and the model is parametrized to be used with an effective path integral derived potential surface, so as to include quantum effects of the transferring proton. The new model is shown to reproduce ab initio geometries and energetics for small protonated clusters. It is also shown that the new model gives a diffusion constant for the excess proton in water, which is in good agreement with experiment, and that the qualitative features of ab initio path integral simulations [D. Marx, M. E. Tuckerman, J. Hutter, and M. Parrinello, Nature (London) 397, 601 (1999)] are well reproduced. [ABSTRACT FROM AUTHOR]

The theoretical principles underpinning the calculation of infrared spectra for condensed-phase systems in the context of ab initio molecular dynamics have been recently developed in literature. At present, most ab initio molecular dynamics calculations are restricted to relatively small systems and short simulation times. In this paper we devise a method that allows well-converged results for infrared spectra from ab initio molecular dynamics simulations using small systems and short trajectories characteristic of simulations typically performed in practice. We demonstrate the utility of our approach by computing the imaginary part of the dielectric constant [variant_greek_epsilon]″(ω) for H2O and D2O in solid and liquid phases and show that it compares well with experimental data. We further demonstrate that maximally localized Wannier orbitals can be used to separate the individual contributions of different molecular species to the linear spectrum of complex systems. The new spectral decomposition method is shown to be useful in present-day ab initio molecular dynamics calculations to compute the magnitude of the “continuous absorption” generated by excess protons in aqueous solutions with good accuracy even when other species present in the solutions absorb strongly in the same frequency window. [ABSTRACT FROM AUTHOR]

There are many atomic scale systems in materials, chemistry, and biology that can be effectively modeled as finite in two of the physical spatial dimensions and periodically replicated in the third including nanoscale metallic and semiconducting wires, carbon nanotubes, and DNA. However, it is difficult to design techniques to treat long range forces in these systems without truncation or recourse to slowly convergent supercells or computationally inefficient Poisson solvers. In this paper, a rigorous reciprocal space based formalism which permits long range forces on wires to be evaluated simply and easily via a small modification of existing methods for three dimensional periodicity is derived. The formalism is applied to determine long range interactions both between point particles using an Ewald-like approach and the continuous charge distributions that appear in electronic structure calculations. In this way, both empirical force field calculations and, for example, plane-wave based density functional theory computations on wires can be performed easily. The methodology is tested on model and realistic systems including a lithium doped carbon nanotube. © 2004 American Institute of Physics. [ABSTRACT FROM AUTHOR]

The ab initio molecular-dynamics formalism of Car and Parrinello is extended to preserve the locality of the orbitals. The supplementary term in the Lagrangian does not affect the nuclear dynamics, but ensures “on the fly” localization of the electronic orbitals within a periodic supercell in the Γ-point approximation. The relationship between the resulting equations of motion and the formation of a gauge-invariant Lagrangian combined with a gauge-fixing procedure is briefly discussed. The equations of motion can be used to generate a very stable and easy to implement numerical integration algorithm. It is demonstrated that this algorithm can be used to compute the trajectory of the maximally localized orbitals, known as Wannier orbitals, in ab initio molecular dynamics with only a modest increase in the overall computer time. In the present paper, the new method is implemented within the generalized gradient approximation to Kohn–Sham density-functional theory employing plane wave basis sets and atomic pseudopotentials. In the course of the presentation, we briefly discuss how the present approach can be combined with localized basis sets to design fast linear scaling ab initio molecular-dynamics methods. © 2004 American Institute of Physics. [ABSTRACT FROM AUTHOR]

In this paper (Paper I) and a companion paper (Paper II), novel new algorithms and applications of the isokinetic ensemble as generated by Gauss' principle of least constraint, pioneered for use with molecular dynamics 20 years ago, are presented for biophysical, path integral, and Car-Parrinello based ab initio molecular dynamics. In Paper I, a new "extended system" version of the isokinetic equations of motion that overcomes the ergodicity problems inherent in the standard approach, is developed using a new theory of non-Hamiltonian phase space analysis [M. E. Tuckerman et al., Europhys. Lett. 45, 149 (1999); J. Chem. Phys. 115, 1678 (2001)]. Reversible multiple time step integrations schemes for the isokinetic methods, first presented by Zhang [J. Chem. Phys. 106, 6102 (1997)] are reviewed. Next, holonomic constraints are incorporated into the isokinetic methodology for use in fast efficient biomolecular simulation studies. Model and realistic examples are presented in order to evaluate, critically, the performance of the new isokinetic molecular dynamic schemes. Comparisons are made to the, now standard, canonical dynamics method, Nosé-Hoover chain dynamics [G. J. Martyna et al., J. Chem. Phys. 97, 2635 (1992)]. The new isokinetic techniques are found to yield more efficient sampling than the Nosé-Hoover chain method in both path integral molecular dynamics and biophysical molecular dynamics calculations. In Paper II, the use of isokinetic methods in Car-Parrinello based ab initio molecular dynamics calculations is presented. [ABSTRACT FROM AUTHOR]

In this paper (Paper II), the isokinetic dynamics scheme described in Paper I is combined with the plane-wave based Car-Parrinello (CP) ab initio molecular dynamics (MD) method [R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985)] to enable the efficient study of chemical reactions and metallic systems. The Car-Parrinello approach employs "on the fly" electronic structure calculations as a means of generating accurate internuclear forces for use in a molecular dynamics simulation. This is accomplished by the introduction of an extended Lagrangian that contains the electronic orbitals as fictitious dynamical variables (often expressed directly in terms of the expansion coefficients of the orbitals in a particular basis set). Thus, rather than quench the expansion coefficients to obtain the ground state energy and nuclear forces at every time step, the orbitals are "propagated" under conditions that allow them to fluctuate rapidly around their global minimum and, hence, generate an accurate approximation to the nuclear forces as the simulation proceeds. Indeed, the CP technique requires the dynamics of the orbitals to be both fast compared to the nuclear degrees of freedom while keeping the fictitious kinetic energy that allows them to be propagated dynamically as small as possible. While these conditions can be easy to achieve in many types of systems, in metals and highly exothermic chemical reactions difficulties arise. (Note, the CP dynamics of metals is incorrect because the nuclear motion does not occur on the ground state electronic surface but it can, nonetheless, provide useful information.) In order to alleviate these difficulties the isokinetic methods of Paper I are applied to derive isokinetic CP equations of motion. The efficacy of the new isokinetic CPMD method is demonstrated on model and realistic systems. The latter include, metallic systems, liquid aluminum, a small silicon sample, the 2 × 1 reconstruction of the silicon 100... [ABSTRACT FROM AUTHOR]

The transport of protons through aqueous, partially aqueous, or nonaqueous hydrogen-bonded media is a fundamental process in many biologically and technologically important systems. Liquid methanol is an example of a hydrogen-bonded system that, like water, supports anomalously fast proton transport. Using the methodology of ab initio molecular dynamics, in which internuclear forces are computed directly from electronic structure calculations as the simulation proceeds, we have investigated the microscopic mechanism of the proton transport process in liquid methanol at 300 K. It is found that the defect structure associated with an excess proton in liquid methanol is a hydrogen-bonded cationic chain whose length generally exceeds the average chain length in pure liquid methanol. Hydrogen bonds in the first and second solvation shells of the excess proton are considerably shorter and stronger than ordinary methanol-methanol hydrogen bonds. Along this chain, proton transfer reactions occur in an essentially random manner described by Poisson statistics. Structural diffusion of the defect structure is possible if the proton migrates toward an end of the defect chain, which causes a weakening of the hydrogen bonds at the opposite end. The latter can, therefore, be easily ruptured by ordinary thermal fluctuations. At the end of the chain where the proton resides, new hydrogen bonds are likely to form due to the strong associative nature of the excess proton. It is through this "snake-like" mechanism that the defect structure is able to diffuse through the hydrogen-bond network of the liquid. The estimated activation enthalpy of this proposed mechanism is found to be in reasonable agreement with the experimentally determined activation enthalpy. [ABSTRACT FROM AUTHOR]

A new formalism designed to treat long range interactions on surfaces, systems which are infinitely replicated in two spatial directions but have finite extent in the third, is developed. The new formalism is based in reciprocal space and, thus, permits the facile extension of standard plane-wave based density functional theory, Ewald summation, and smooth particle-mesh Ewald methods to handle surfaces efficiently. The method is tested on both model (body centered cubic lattices) and realistic problems (an ice surface with a defect and the 2 x 1 surface reconstruction of silicon) and found to be accurate, efficient, and a marked improvement on existing formulations in speed, accuracy, and utility. [ABSTRACT FROM AUTHOR]

A new molecular dynamics method for calculating free energy profiles for rare events is presented. The new method is based on the creation of an adiabatic separation between a reaction coordinate subspace and the remaining degrees of freedom within a molecular dynamics run. This is achieved by associating with the reaction coordinate(s) a high temperature and large mass, thereby allowing the activated process to occur while permitting the remaining degrees of freedom to respond adiabatically. In this limit, by applying a formal multiple time scale Liouville operator factorization, it can be rigorously shown that the free energy profile is obtained directly from the probability distribution of the reaction coordinate subspace and, therefore, no unbiasing of the configuration space or postprocessing of the output data is required. The new method is applied to a variety of model problems including a two-dimensional free energy surface and its performance tested against free energy calculations using the “blue moon ensemble” approach. The comparison shows that free energy profiles can be calculated with greater ease and efficiency using the new method. © 2002 American Institute of Physics. [ABSTRACT FROM AUTHOR]

Presents a reciprocal space based formalism for treating long range forces in clusters. Incorporation of the formalism into plane-wave based density function theory calculations, Ewald summation calculations and smooth particle-mesh Ewald calculations to yield accurate and numerically efficient descriptions of long range interactions in cluster systems.

Investigates the double proton transfer in the formic acid dimer using ab initio path integral molecular dynamics (PIMD). Description of the electronic structure by Kohn-Sham formulation of the density functional theory; Details on the PIMD approach of Tuckerman; Results of the thermal and quantum fluctuations around the equilibrium structure.

Examines the empirical valence model for proton transfer in aqueous solutions. Properties of the gas-phase H5O2+ complex; Use of gradient-corrected density functional theory in the study; Simulations of the solvated hydronium.

Focuses on the nonequilibrium molecular dynamics generated by the SLLOD algorithm for fluid flow. Existence of a conserved energy for the equations of motion; View that coupling each degree of freedom in the system to a Nose-Hoover chain thermostat as a means of efficiently generating the phase space equilibrium.

In path integral molecular dynamics, efficient sampling of the phase space is not guaranteed due to the stiff harmonic part of the action arising from the quantum kinetic energy. This problem has been eliminated by incorporating a sufficient number of thermostats into the dynamical scheme and by introducing a transformation of the path ‘‘bead’’ variables. In this paper, an efficient Car–Parrinello path integral molecular dynamics algorithm, sufficiently general to include the use of ultrasoft pseudopotentials is introduced. Difficulties encountered when combining thermostats and transformations of the Cartesian ‘‘bead’’ coordinates with the generalized orthonormality condition are circumvented by employing a constrained nonorthogonal orbital method. © 1996 American Institute of Physics. [ABSTRACT FROM AUTHOR]

A new fourth order predictor–corrector integration scheme is presented. The unique feature of the new algorithm and what distinguishes it from a Gear predictor–corrector is that the method is derived from the Trotter decomposition of a specially formulated evolution operator and as such, is both symplectic and reversible. In addition, the method retains the useful property of Gear methods that only one force evaluation per time step is required. The new integrator is tested on a harmonic plus quartic oscillator and the Henon–Heiles system. Comparisons are made to the second order velocity Verlet integrator, the true fourth order Yoshida/Suzuki schemes and fourth order Gear. In all cases, the new method works well, giving energy conservation and trajectories of much better quality than velocity Verlet and of comparable quality to the results of the true fourth order schemes for the same computational cost as velocity Verlet. © 1995 American Institute of Physics. [ABSTRACT FROM AUTHOR]

In this paper, new techniques for integrating the Car–Parrinello equations based on multiple time scale methodology are presented. The formalism of multiple time scale methodology based on operator factorizations of the classical Liouville propagator is reviewed. It is shown how the techniques are applied to Car–Parrinello for use with the velocity Verlet and Gaussian dynamics schemes presented in the preceding paper [M. E. Tuckerman and M. Parrinello, J. Chem. Phys. 101, 1302 (1994)], and a detailed discussion is presented of how a reference system for Car–Parrinello simulations may be chosen. It is shown that the use of such techniques can save up to a factor of 5–10 in cpu time over the standard methods. [ABSTRACT FROM AUTHOR]

In this paper and in a companion paper [M. E. Tuckerman and M. Parrinello, J. Chem. Phys. 101, 1316 (1994)] the problem of integrating the equations of motion in Car–Parrinello simulations is addressed. In this paper, new techniques for treating the constraint problem based on the velocity Verlet integrator and the Gaussian dynamics are presented. Questions of adiabaticity and temperature control are discussed, and it is shown how to combine the new techniques with the recently developed Nosé–Hoover chain thermostat method. All new techniques are described using the formalism of operator factorizations applied to the classical Liouville propagator. In the companion paper, the formalism and application of multiple time scale methodology in Car–Parrinello simulations are discussed. [ABSTRACT FROM AUTHOR]

The reference system propagator algorithm of Tuckerman et al. is applied to a simple model of flexible CO2. The forces are broken up according to three different time scales corresponding to the internal motion, short-range impulsive collisions, and long-range interactions, and three time steps are introduced to integrate the equations of motion. The use of this scheme allows a factor of 8 to 12 to be saved in CPU time over standard integration schemes. [ABSTRACT FROM AUTHOR]

New path integral molecular dynamics (PIMD) and path integral hybrid Monte Carlo (PIHMC) algorithms are developed. It is shown that the use of a simple noncanonical change of variables that naturally divides the quadratic part of the action into long and short wavelength modes and multiple time scale integration techniques results in very efficient algorithms. The PIMD method also employs a constant temperature MD technique that has been shown to give canonical averages even for stiff systems. The new methods are applied to the simple quantum mechanical harmonic oscillator and to electron solvation in fluid helium and xenon. Comparisons are made with PIMC and the more basic PIMD and PIHMC methods. [ABSTRACT FROM AUTHOR]

Rate constants evaluated from (1) the energy-loss turnover theory of Pollak, Grabert, and Hänggi (PGH), (2) the Grote–Hynes extension of Kramers theory (GH), and (3) the microcanonical variational transition state theory for dissipative systems of Tucker and Pollak (μVTST) are compared with rate constants determined from direct computer simulations of generalized Langevin dynamics. The comparisons are made for a cubic oscillator under the influence of a slow bath characterized by a Gaussian friction kernel. In the μVTST calculations, which are based on an effective two degree of freedom Hamiltonian, barrier crossing due to energy transfer from the bath to the effective Hamiltonian is neglected. This neglect is significant only at very strong coupling, where it causes the μVTST results to drop below the simulation results. Both GH and μVTST theories fail (as expected) in the energy diffusion regime, while PGH theory is only moderately successful. The μVTST results agree extremely well with the simulation results in the spatial diffusion regime, providing a significant improvement over the GH results at intermediate coupling strengths and over the PGH results at strong coupling strengths. This improvement is a result of nonlinear effects which are included in the μVTST approach but neglected in the PGH and GH theories. [ABSTRACT FROM AUTHOR]

The generalized Langevin equation (GLE) has been used to model a wide variety of systems in which a subset of the degrees of freedom move on a potential of mean force surface subject to fluctuating forces and dynamic friction. When there is a wide separation in the time scales for motion on the potential surface and for relaxation of the friction kernel, direct integration of the GLE is very costly in CPU time. In this paper we introduce an integrator based on our previous work using numerical analytical propagator algorithm (NAPA) and reference system propagator algorithm (RESPA) that greatly accelerates such simulations. We also discuss sampling methods for the random force. Accuracy of this algorithm is assessed by comparisons with an analytically solvable example. Introducing dynamic friction kernels determined from full molecular dynamics (MD) simulations allows us to compare the accuracy of the GLE simulations with full scale molecular dynamics simulations. [ABSTRACT FROM AUTHOR]

A frequently encountered problem in molecular dynamics is how to treat the long times that are required to simulate condensed systems consisting of particles interacting through long range forces. Standard methods require the calculation of the forces at every time step. Because each particle interacts with all particles within the interaction range of the potential the longer the range of the potential the larger the number forces that must be calculated at each time step. In this note we present a variant of the RESPA (reference system propagator algorithm), which we developed for handling systems with multiple time scales like disparate mass mixtures. This version of RESPA greatly reduces the number of forces that must be computed at each time step and thereby leads to a dramatic acceleration of such simulations. The RESPA method uses ideas similar to NAPA, an algorithm we invented to treat high frequency oscillators interacting with low frequency bath. The method is based on a choice of a reference system in which the particles interact through short range forces. The reference system is numerically integrated for n time steps δt and the error incurred by using short range forces is corrected by solving a rigorous set of equations once every Δt=nδt. This method reduces the cpu time dramatically. It is shown that this approach and suitable generalizations should be very useful for future simulations of quantum and classical condensed matter systems. [ABSTRACT FROM AUTHOR]

A frequently encountered problem in molecular dynamics is how to treat the long times that are required to simulate condensed systems consisting of mixtures of light and heavy particles. Standard methods require the choice of time step sufficiently small to guarantee stable solution for the low mass component with the consequence that these simulations require a very large number of central processing unit cycles to treat the relaxation of the heavier component. In this note, we present a new method that allows one to use a time step appropriate for the heavy particles. This method uses a similar idea to numerical analytical propogator algorithm, an algorithm we invented to treat high frequency oscillators interacting with low frequency baths and is based on a choice of a reference system for the light particle motions. The method is applied to the case of a liquid containing 864 Lennard-Jones spheres, 824 of these particles having a mass, M=100 and 40 spheres picked at random have a mass m=1. It is shown that molecular dynamics using the new algorithm runs seven to ten times faster than standard methods and this approach as well as suitable generalizations should be very useful for future simulations of quantum and classical condensed matter systems. [ABSTRACT FROM AUTHOR]

A frequently encountered problem in molecular dynamics simulations is the long runs required to study condensed systems consisting of both high frequency and low frequency degrees of freedom. Standard integrators require the choice of time step sufficiently small to guarantee stable solution of the highest frequency motion with the consequence that simulations require a very large number of central processing unit (CPU) cycles. In this note we present a new integrator that allows one to use a time step appropriate for the low frequency degrees of freedom without making any approximations related to the separation of time scales. This method is based on a choice of an analytically solvable reference system for the high frequency motion. We show how the analytical solution can be incorporated into a numerical integrator. The method is applied to two cases which are paradigms for this problem and it is shown that this approach and suitable generalizations should be very useful for future simulations of quantum and classical condensed matter systems. [ABSTRACT FROM AUTHOR]

Comments on the scientific paper 'Simple Reversible Molecular Dynamics Algorithms for Nose-Hoover Chain Dynamics,' by S. Jang and J.A. Voth. Criticism on the proposed leap-frog algorithms of Jang and Voth; Failure to consider energy conservation in the Nose-Hoover chain dynamics.

When there are high and low frequency motions in systems with long and short range forces a judicious choice of reference system leads to very large accelerations in molecular dynamic (MD) simulations. Building on our previous work where we have developed reference system methods for systems with high frequency oscillators, disparate masses, or long range forces, we present a double reference system method which leads to acceleration of as much as 20 in systems consisting of 864 molecules with Lennard-Jones (12-6) forces. Much larger savings should be achieved when this method is applied to longer range forces and larger systems. [ABSTRACT FROM AUTHOR]

We introduce a scheme for approximating quantum time correlation functions numerically within the Feynman path integral formulation. Starting with the symmetrized version of the correlation function expressed as a discretized path integral, we introduce a change of integration variables often used in the derivation of trajectory-based semiclassical methods. In particular, we transform to sum and difference variables between forward and backward complex-time propagation paths. Once the transformation is performed, the potential energy is expanded in powers of the difference variables, which allows us to perform the integrals over these variables analytically. The manner in which this procedure is carried out results in an open-chain path integral (in the remaining sum variables) with a modified potential that is evaluated using imaginary-time path-integral sampling rather than requiring the generation of a large ensemble of trajectories. Consequently, any number of path integral sampling schemes can be employed to compute the remaining path integral, including Monte Carlo, path-integral molecular dynamics, or enhanced path-integral molecular dynamics. We believe that this approach constitutes a different perspective in semiclassical-type approximations to quantum time correlation functions. Importantly, we argue that our approximation can be systematically improved within a cumulant expansion formalism. We test this approximation on a set of one-dimensional problems that are commonly used to benchmark approximate quantum dynamical schemes. We show that the method is at least as accurate as the popular ring-polymer molecular dynamics technique and linearized semiclassical initial value representation for correlation functions of linear operators in most of these examples and improves the accuracy of correlation functions of nonlinear operators. [ABSTRACT FROM AUTHOR]

Responds to the criticism made on an article about the nonequilibrium molecular dynamics of fluid flow. Basis of the nonequilibrium steady state of the Liouville equation; Precision temperature control in chain thermostats; Stabilization of linear velocity profile in the absence of boundary conditions.

We correct an error in our original paper. Calculations originally reported for the alanine tripeptide were performed using the earlier driven adiabatic free energy dynamics (d-AFED) approach rather than the new unified free energy dynamics (UFED) scheme. When UFED is applied to this system, the rate of convergence of the free energy surface improves significantly over that originally reported for d-AFED. [ABSTRACT FROM AUTHOR]

It is generally believed that studies of liquid water using the generalized gradient approximation to density functional theory require dispersion corrections in order to obtain reasonably accurate structural and dynamical properties. Here, we report on an ab initio molecular dynamics study of water in the isothermal-isobaric ensemble using a converged discrete variable representation basis set and an empirical dispersion correction due to Grimme [J. Comp. Chem. 27, 1787 (2006)]. At 300 K and an applied pressure of 1 bar, the density obtained without dispersion corrections is approximately 0.92 g/cm3 while that obtained with dispersion corrections is 1.07 g/cm3, indicating that the empirical dispersion correction overestimates the density by almost as much as it is underestimated without the correction for this converged basis. Radial distribution functions exhibit a loss of structure in the second solvation shell. Comparison of our results with other studies using the same empirical correction suggests the cause of the discrepancy: the Grimme dispersion correction is parameterized for use with a particular basis set; this parameterization is sensitive to this choice and, therefore, is not transferable to other basis sets. [ABSTRACT FROM AUTHOR]

We propose a general framework for the efficient sampling of conformational equilibria in complex systems and the generation of associated free energy hypersurfaces in terms of a set of collective variables. The method is a strategic synthesis of the adiabatic free energy dynamics approach, previously introduced by us and others, and existing schemes using Gaussian-based adaptive bias potentials to disfavor previously visited regions. In addition, we suggest sampling the thermodynamic force instead of the probability density to reconstruct the free energy hypersurface. All these elements are combined into a robust extended phase-space formalism that can be easily incorporated into existing molecular dynamics packages. The unified scheme is shown to outperform both metadynamics and adiabatic free energy dynamics in generating two-dimensional free energy surfaces for several example cases including the alanine dipeptide in the gas and aqueous phases and the met-enkephalin oligopeptide. In addition, the method can efficiently generate higher dimensional free energy landscapes, which we demonstrate by calculating a four-dimensional surface in the Ramachandran angles of the gas-phase alanine tripeptide. [ABSTRACT FROM AUTHOR]

Mixed ab initio/empirical force-field simulation studies, calculations in which one part of the system is treated using a fully ab initio description and another part is treated using an empirical description, are becoming increasingly popular. Here, the ability of the commonly used, plane wave-based generalized gradient approximation to density functional theory is extended to model systems in which the electrons are assumed to be localized in a single small region of space, that is, itself, embedded within a large chemically inert bath. This is accomplished by introducing two length scales, so that the rapidly varying, short range, electron–electron and electron–atom interactions, arising from the region where the electrons are localized, can be treated using an appropriately large plane wave basis, while the corresponding, slowly varying, long range interactions of the electrons with the full system or bath, can be treated using a small basis. Briefly, a novel Cardinal B-spline based formalism is employed to derive a smooth, differentiable, and rapidly convergent (with respect to the small basis) expression for the total electronic energy, which explicitly contains the two length scales. The method allows reciprocal space based techniques designed to treat clusters, wires, surfaces and solids/liquids (open, and 1-D and 2-D periodic boundary conditions, respectively) to be utilized. Other plane wave-based “mixed” methods are restricted to clusters. The new methodology, which scales as N log N at fixed size of the chemically active region, has been implemented for parallel computing platforms and tested through applications to both model and realistic problems including an enzyme, human carbonic anhydrase II solvated in an explicit bath of water molecules. © 2001 American Institute of Physics. [ABSTRACT FROM AUTHOR]

The use of non-Hamiltonian dynamical systems to perform molecular dynamics simulation studies is becoming standard. However, the lack of a sound statistical mechanical foundation for non-Hamiltonian systems has caused numerous misconceptions about the phase space distribution functions generated by these systems to appear in the literature. Recently, a rigorous classical statistical mechanical theory of non-Hamiltonian systems has been derived, [M. E. Tuckerman, et al., Europhys. Lett. 45, 149 (1999)]. In this paper, the new theoretical formulation is employed to develop the non-Hamiltonian generalization of the usual Hamiltonian based statistical mechanical phase space principles. In particular, it is shown how the invariant phase space measure and the complete sets of conservation laws of the dynamical system can be combined with the generalized Liouville equation for non-Hamiltonian systems to produce a well defined expression for the phase space distribution function. The generalization provides a systematic, controlled procedure for designing non-Hamiltonian molecular dynamics algorithms which can be used to generate nonmicrocanonical ensembles, stationary nonequilibrium flows, and/or the dynamics of constrained systems. In light of this new general analysis, molecular dynamics algorithms for the canonical and isothermal–isobaric ensembles are examined, potential difficulties are illuminated, and the limitations of previous theoretical treatments are elucidated. © 2001 American Institute of Physics. [ABSTRACT FROM AUTHOR]

A new method for generating the canonical ensemble via continuous dynamics is presented. The new method is based on controlling the fluctuations of an arbitrary number of moments of the multidimensional Gaussian momentum distribution function. The equations of motion are non-Hamiltonian, and hence have a nonvanishing phase space compressibility. By applying the statistical mechanical theory of non-Hamiltonian systems recently introduced by the authors [M. E. Tuckerman, C. J. Mundy, and G. J. Martyna, Europhys. Lett. 45, 149 (1999)], the equations are shown to produce the correct canonical phase space distribution function. Reversible integrators for the new equations of motion are derived based on a Trotter-type factorization of the classical Liouville propagator. The new method is applied to a variety of simple one-dimensional example problems and is shown to generate ergodic trajectories and correct canonical distribution functions of both position and momentum. The new method is further shown to lead to rapid convergence in molecular dynamics based calculations of path integrals. The performance of the new method in these examples is compared to that of another canonical dynamics method, the Nosé-Hoover chain method [G. J. Martyna, M. L. Klein, and M. E. Tuckerman, J. Chem. Phys. 97, 2635 (1992)]. The comparison demonstrates the improvements afforded by the new method as a molecular dynamics tool. Finally, when employed in molecular dynamics simulations of biological macromolecules, the new method is shown to provide better energy equipartitioning and temperature control and to lead to improved spatial sampling over the Nosé-Hoover chain method in a realistic application. © 2000 American Institute of Physics. [ABSTRACT FROM AUTHOR]

Presents corrections to the article 'Molecular Dynamics Algorithm for Multiple Time Scales: Systems with Disparate Masses,' by Mark E. Tuckerman, Bruce J. Berne and Angelo Rossi from the 1991 issue of the 'Journal of Chemical Physics.'

Dynamical properties of liquid water were studied using Car-Parrinello [Phys. Rev. Lett. 55, 2471 (1985)] ab initio molecular dynamics (AIMD) simulations within the Kohn-Sham (KS) density functional theory employing the Becke-Lee-Yang-Parr exchange-correlation functional for the electronic structure. The KS orbitals were expanded in a discrete variable representation basis set, wherein the complete basis set limit can be easily reached and which, therefore, provides complete convergence of ionic forces. In order to minimize possible nonergodic behavior of the simulated water system in a constant energy (NVE) ensemble, a long equilibration run (30 ps) preceded a 60 ps long production run. The temperature drift during the entire 60 ps trajectory was found to be minimal. The diffusion coefficient [0.055 A2/ps] obtained from the present work for 32 D2O molecules is a factor of 4 smaller than the most up to date experimental value, but significantly larger than those of other recent AIMD studies. Adjusting the experimental result so as to match the finite-sized system used in the present study brings the comparison between theory and experiment to within a factor of 3. More importantly, the system is not observed to become "glassy" as has been reported in previous AIMD studies. The computed infrared spectrum is in good agreement with experimental data, especially in the low frequency regime where the translational and librational motions of water are manifested. The long simulation length also made it possible to perform detailed studies of hydrogen bond dynamics. The relaxation dynamics of hydrogen bonds observed in the present AIMD simulation is slower than those of popular force fields, such as the TIP4P potential, but comparable to that of the TIP5P potential.