Molecular Dynamics Simulations: The Limits and Beyond* Herman J.C.Berendsen BIOSON Research Institute and Dept of Biophysical Chemistry,University of Groningen,Nijenborgh 4,9747 AG Groningen,the Netherlands Abstract.This article reviews the present state of Molecular Dynamics (MD) simulations and tries to give an outlook into future developments.First an overview is given of methods,algorithms and force fields.After considering the limitations of the standard present-day techniques,developments that reach beyond the present limitations are considered.These concern three major directions:(a)inclusion of quantum dynamics,(b)reduction of complexity by reducing the number of degrees of freedom and averaging over interactions with less important degrees of freedom, (c)reduction to mesoscopic dynamics by considering particle densities rather than positions.It is concluded that MD is a mature technique for classical simulations of all-atom systems in the nanosecond time range,but is still in its infancy in reaching reliably into longer time scales. 1 Introduction and a Bit of History This conference on Algorithms for Molecular Simulation is a good occasion to pause for a moment and consider where we are now and where we go in this field.The book by Allen and Tildesley [2]describes most of the techniques that are still in use today. Molecular Dynamics (MD)simulations were first carried out in 1957 by Alder and Wainwright on a hard-sphere fluid 2];the first fluid with soft interactions was simulated by Rahman in 1964 [3]and the first complex fluid (water)was simulated by Rahman and Stillinger in 1971 [4].The first MD simulation of a protein was carried out in 1976 by Andrew McCammon [2],then postdoc in Martin Karplus'group,during a CECAM workshop in Orsay,France.That workshop [6],which brought about 20 physicists (one was the 'father'of MD,Anees Rahman)and protein specialists (one was Jan Hermans)together for two months,has been a seminal event in the development of simulation of biological macromolecules.Since then methods and force fields have improved and computers have become a thousandfold more powerful.Routine simulations now comprise fully hydrated systems with tens of thousands of atoms and extend over nanoseconds.Simulated systems include DNA,liquid crystals,polymers and lipid membranes. *A contribution from the Groningen Biomolecular Sciences and Biotechnology Institute
4 Berendsen But the methods have not really changed.The Verlet algorithm to solve Newton's equations,introduced by Verlet in 1967 [7],and it's variants are still the most popular algorithms today,possibly because they are time-reversible and symplectic,but surely because they are simple.The force field descrip- tion was then,and still is,a combination of Lennard-Jones and Coulombic terms,with(mostly)harmonic bonds and periodic dihedrals.Modern exten- sions have added many more parameters but only modestly more reliability. The now almost universal use of constraints for bonds (and sometimes bond angles)was already introduced in 1977 [8].That polarisability would be nec- essary was realized then [9],but it is still not routinely implemented today. Long-range interactions are still troublesome,but the methods that now be- come popular date back to Ewald in 1921 10 and Hockney and Eastwood in198111. What has been developed within the last 20 years is the computation of thermodynamic properties including free energy and entropy [12,13,14].But the ground work for free energy perturbation was done by Valleau and Torrie in 1977 [15],for particle insertion by Widom in 1963 and 1982 [16,17]and for umbrella sampling by Torrie and Valleau in 1974 and 1977 [18,19.These methods were primarily developed for use with Monte Carlo simulations; continuous thermodynamic integration in MD was first described in 1986 [201. Another topic that received increasing attention is the incorporation of quantum methods into dynamic simulations.True quantum dynamics for hundreds of particles is beyond any foreseeable computational capability,and only approximations are viable.We should distinguish: ()The application of quantum corrections to classical MD.An early exam- ple is the application of quantum corrections to water based on classical frequency distributions by Behrens et al.[21]. (i)The use of quantum methods to derive potentials for the heavy particles in the Born-Oppenheimer approximation during the MD simulation.This is now a very active field,with important applications for the study of chemical reactions in the condensed phase.Pioneering work using semi- empirical QM was done by Warshel [22],and in the groups of Jorgensen [23],Kollman [24],and Karplus [25,26].A methodological breakthrough was the introduction in 1985 of DFT(quantum density functional theory) to solve instantaneous forces on the heavy particles in a method that is now called ab initio molecular dymamics,by Car and Parrinello [27. This method solves the dynamical evolution of the ground-state electronic wave function described as a linear combination of plane waves.Selloni et al.28 solved the dynamical evolution of the ground state wave function of a solvated electron on a 3-D grid.These methods are adiabatic in the sense that they do not incorporate excitations of the quantum particles. (i)The use of quantum methods to obtain correct statistical static (but not dynamic)averages for heavy'quantum particles.In this category path-integral methods were developed on the basis of Feynman's path
MD:The Limits and Beyond 5 integral formulation [29]mostly by Chandler and Wolynes [30].For a lucid description see [31]. (iv)The incorporation of quantum-dynamical evolution for selected particles or degrees of freedom in a classical environment.With the inclusion of non-adiabatic transfers to excited states this is a rather new field,with im- portant applications to proton transfer processes in the condensed phase. We shall return to these methods in section 3.2. The important enquiry into long time scales has also been a subject of interest over many years,but the progress has been slow.Computer capabili- ties have increased so rapidly that it has often been worthwhile to wait a few years to obtain the required increase in speed with standard methods rather than invent marginal improvements by faster algorithms or by using reduced systems.Many attempts to replace the time-consuming solvent molecules by potentials of mean force(see for example [32])or to construct an appropriate outer boundary without severe boundary effects [43,34]have been made, but none of these are fully satisfactory.Really slow events cannot be mod- eled by such simplifications:a drastic reduction in the number of degrees of freedom is needed.When events are slow because an identifiable barrier must be crossed,good results can be obtained by calculating the free energy at the barrier in one or a few degrees of freedom.However,when events are slow because a very large multidimensional configurational space must be explored (as in protein folding or macromolecular aggregation),the appro- priate methods are still lacking.We shall return to this important topic in Section 3.3. 2 Where Are We Now? With the danger of severe oversimplification,which unavoidably leads to improper underevaluation of important recent developments,I shall try to indicate where traditional,classical MD has brought us today,or will bring us tomorrow.This concerns the techniques rather than the applications,which cannot be reviewed in the present context.The main aspects to consider concern algorithms and force fields. 2.1 Algorithms As remarked in the introduction,the reversible Verlet algorithm or any of its disguised forms as velocity-Verlet or leap-frog,has remained strong and sturdy.Non-reversible higher-order algorithms of the predictor-corrector type,such as Gear's algorithms,may be useful if very high accuracy is re- quired,but offer little advantage in cases where the evaluation of forces is accompanied by noisy errors [35]. An elegant derivation of the Verlet-type algorithms has been given by Tuckerman et al.[36]and is useful in multiple timestep implementations
6 Berendsen called RESPA (REference System Propagator Algorithms)[37,20].Because of their elegance I cannot resist to quote the principle of these algorithms. They are based on the Trotter-Suzuki expansion of exponential operators [39] which is much used in quantum simulations: ei(A+B)teiAt/2eiBteiAt/2 where A and B are non-commuting operators.This expansion is obviously time reversible and gives an error of third order in t.If applied to the classical Liouville operator acting on phase space,and separating the components act- ing on coordinates and on momenta,equations for momenta and coordinate updates per time step are obtained.Let us make it simple:In cartesian coor- dinates a point in phase space is represented by a vector of positions x and velocities v.Evolving the vector(x,v)T over a time(step)t means applying both a force propagator (）=(+e/m） and a velocity propagator ()-() Each of these operators is unitary:U(-t)=U-1(t).Updating a time step with the propagator Uf(At)U(At)Uf(At)yields the velocity-Verlet algo- rithm.Concatenating the force operator for successive steps yields the leap- frog algorithm: )v-+F(At/m 1 z(t+At)=(t)+(+t)At A double-timestep algorithm with short-and long-range forces is obtained by applying the propagator [36]: w.((tU.((), where Us and U are the propagators for the short-range,resp.the long-range force,and At =not.These algorithms are not only time reversible,but they are sympletic and conserve volume in phase space [40]. In practice modifications are made to incorporate thermostats or barostats that may destroy the time-reversible and symplectic properties. While extended-system algorithms such as Nose dynamics [41]can be de- signed on the principles of the reversible operators,methods that use pro- portional velocity or coordinate scaling [42]cannot.Such methods are very
MD:The Limits and Beyond 7 convenient and practical,but -unless the time constant used for coupling with an external bath are long with respect to the intrinsic time constant for kinetic energy dispersal -they give undefined ensembles with unreliable fluctuations and may produce spurious transfer of kinetic energy to uncou- pled degrees of freedom,such as the overall translation and rotation of the system.Standard programs correct for such effects. The incorporation of holonomic constraints for covalent bond lengths(and sometimes bond angles)saves roughly a factor of 4 in the allowed time step for molecular systems and has been common practice for many years.Conserving constraints involves the solution of a set of nonlinear equations,which can be solved iteratively,either by solving a matrix equation after linearization, or by iteratively solving successive equations for each constraint.The latter method is employed in the widely used SHAKE program [8].Recently a linear constraint solver LINCS has been introduced [43 which is much faster and more stable than SHAKE and is better suited for implementation in programs for parallel computers.It is built in our MD package GROMACS![44]. The question whether constraints for covalent bonds give better physics than harmonic oscillators is not really resolved.A mathematical argument can be given that specific motions which occur in a frequency range clearly separated from all other motions can be considered uncoupled and can then be treated as constraints without affecting the overall motion.For molecular systems such a separation is valid for bond length constraints,but not gener- ally for bond angle constraints;the latter should therefore not be constrained in large molecules [3].One should in principle take care of corrections due to the Jacobian of the transformation when using constraints,related to the con- figuration dependence of the extent of phase space on a curved constrained surface [46,47],but this effect is negligible if only bond lengths are con- strained.A physical argument for using bond constraints is that real covalent bonds correspond to quantum-mechanical harmonic oscillators with frequen- cies well above keT/h.They are thus permanently in the stationary ground state and do not take up any additional kinetic or potential energy.Treat- ing them as classical oscillators would provoke unphysical energy exchange and requires the application of quantum corrections to the energy.While this is true,the treatment of bonds as constraints denies the generation of low- frequency modes from coupled vibrations,and also prevents the oscillators to relax under the influence of an external force.The latter effect,if adiabat- ically imposed,absorbs an energy of/(where F is the external force and k is the harmonic force constant)from the environment,both for classi- cal and quantum oscillators.For an OH stretch in water,with external forces around 10-9 N/m [48]this energy is about 0.5 kJ/mol.Finally,both clas- sical vibrations and constraints neglect the configuration dependence of the zero-point energy of quantum oscillators;this effect amounts in water (with a 300 cm-1 shift in vibrational frequency)to 1.8 kJ/mol per OH oscillator. 1 See http://rugmdo.chem.rug.nl/gromacs/