MD:The Limits and Beyond 13 boundary conditions),the Ewald summation method is recovered.If the Fourier solution is obtained on a grid,allowing the use of FFT,the Ewald- mesh method is obtained.However,the charge spread function need not be gaussian and can be chosen such that the short-range field and its derivative go exactly to zero at the cut-off radius.The shifted force ob- tained with the tin-foil reaction field(see above)corresponds to a charge spread function that spreads each charge homogeneously over a sphere with radius re,which is not an optimal choice to avoid truncation errors. Iterative Poisson solvers on a regular grid using Gauss-Seidel iteration with successive overrelaxation (SOR)[82]are less efficient than FFT methods,but are also applicable to non-periodic systems and are more easily parallellized.In an MD run,the previous step provides a near- solution and only a few iterations are needed per step 83. Summarizing,the most efficient way to handle electrostatic interactions cor- rectly seems to be the appropriate splitting into smooth short-range and long-range parts,and handling the latter by an efficient Poisson solver,using the knowledge available from the history of the trajectory,and exploiting the fact that the long-range part fluctuates on a longer time scale. 3 The Limits and Beyond 3.1 Limits to Traditional MD Despite all the shortcomings listed above,full particle classical MD can be considered mature [84].Even when all shortcomings will be overcome,we can now clearly delineate the limits for application.These are mainly in the size of the system and the length of the possible simulation.With the rapidly growing cheap computer memory shear size by itself is hardly a limitation: several tens of thousands of particles can be handled routinely (for example, we report a simulation of a porin trimer protein embedded in a phospholipid membrane in aqueous environment with almost 70,000 particles (85;see also the contribution of K.Schulten in this symposium)and a million particles could be handled should that be desired. The main limitation is in simulated time,which at present is in the order of nanoseconds for large systems.We may expect computational capabilities to increase by at least an order of magnitude every five years,but even at that rate it will take 15 years before we reach routinely into the microseconds. While adequate for many problems,such time scales are totally inadequate for many more other problems,particularly in molecular processes involving biological macromolecules.The obvious example of this is protein folding, which in favourable cases occurs in the real world on a time scale of seconds, to be reached in computero by the year 2040 if force fields are improved ac- cordingly.The latter condition will not automatically be fulfilled:an error in the total energy difference between native and unfolded state of only 8
14 Berendsen kJ/mol (on a total energy difference of some 250 kJ/mol,and with system energies measured in hundreds of MJ/mol)would shift a predicted 'melting' temperature by 10 C.It is instructive to realize that the Born correction in water for a univalent ion for a cut-off radius of 1 nm amounts to 70 kJ/mol, which clearly shows that force field improvements must be accompanied by very careful evaluation of long-range interactions.The incorporation of po- larisability is imperative if we realize that the polarization energy of a single carbon atom at a distance of 0.4 nm from a unit charge amounts to about 5 kJ/mol. The nanoseconds limit also indicates a limit in the configurational sam- pling that can be achieved by MD.Sufficient sampling of the configurational space accessible in an equilibrium condition is essential for the computation of thermodynamic properties that involve entropy,as the latter is a measure for the extent of accessible configurational space.Use of other equilibrium sampling techniques,like Monte Carlo simulation,does not really improve on the statistics.However,substantial improvements are obtained and still to be expected from multiconfigurational sampling,umbrella sampling,and other methods that bias sampling to include infrequently visited regions3, and from methods that circumvent barriers in configurational space such as the modification ('softening')of potential functions and the introduction of a fourth spatial dimension [86,14]. 3.2 Inclusion of Quantum Dynamics A limitation of classical force field-based MD is the restriction to covalent complexes,with exclusion of chemical reactions.The very important appli- cations to reactions in the condensed phase,including enzyme reactions and catalysis in general,need extension with dynamic behaviour of non-covalent intermediates.The latter must be described by quantum-mechanical meth- ods.Usually,except for fast electron transfer reactions,the Born-Oppenhei- mer approximation is valid for the electronic motion.Classical dynamics is generally sufficiently accurate for atoms heavier than hydrogen,but for pro- ton transfer reactions explicit quantum-dynamical treatment of the proton is required which fully includes tunneling as well as non-adiabatic involvement of excited states.We shall separately consider the two aspects:(i)computing the reactive Born-Oppenheimer surface in large condensed systems,and (ii) predicting reaction dynamics,including the quantum behaviour of protons. As ab initio MD for all valence electrons [27]is not feasible for very large systems,QM calculations of an embedded quantum subsystem are required. Since reviews of the various approaches that rely on the Born-Oppenheimer approximation and that are now in use or in development,are available (see Field [87],Merz [88],Aqvist and Warshel [89],and Bakowies and Thiel [90] and references therein),only some summarizing opinions will be given here. 3 see for example the contribution by Grubmuiller in this symposium
MD:The Limits and Beyond 15 The first point to remark is that methods that are to be incorporated in MD,and thus require frequent updates,must be both accurate and efficient. It is likely that only semi-empirical and density functional (DFT)methods are suitable for embedding.Semi-empirical methods include MO (molecular orbital)[90]and valence-bond methods [89],both being dependent on suit- able parametrizations that can be validated by high-level ab initio QM.The quality of DFT has improved recently by refinements of the exchange den- sity functional to such an extent that its accuracy rivals that of the best ab initio calculations [91].DFT is quite suitable for embedding into a classical environment [92].Therefore DFT is expected to have the best potential for future incorporation in embedded QM/MD. The second aspect,predicting reaction dynamics,including the quantum behaviour of protons,still has some way to go!There are really two sepa- rate problems:the simulation of a slow activated event,and the quantum- dynamical aspects of a reactive transition.Only fast reactions,occurring on the pico-to nanosecond time scale,can be probed by direct simulation;an in- teresting example is the simulation by ab initio MD of metallocene-catalysed cthylene polymerisation by Meier et al.93. Most reactions are too slow on a time scale of direct simulation,and the evaluation of reaction rates then requires the identification of a transition state (saddle point)in a reduced space of a few degrees of freedom(reaction coordinates),together with the assumption of equilibration among all other degrees of freedom.What is needed is the evaluation of the potential of mean force in this reduced space,using any of the available techniques to com- pute free energies.This defines the probability that the system resides in the transition-state region.Then the reactive flux in the direction of the unstable mode at the saddle point must be calculated.If friction for the motion over the barrier is neglected,a rate according to Eyrings transition-state theory is obtained.In general,the rate is smaller due to unsuccessful barrier crossing, as was first described by Kramers [94].The classical transition rate can be properly treated by the reactive flux method [95],see also the extensive re- view by Hanggi [96].The reactive flux can be obtained from MD simulations that start from the saddle point.An illustrative and careful application of the computational approach to classical barrier crossing,including a discussion of the effects due to the Jacobian of the transformation to reaction coordinates, has recently been described by den Otter and Briels [47]. While the classical approach to simulation of slow activated events,as described above,has received extensive attention in the literature and the methods are in general well established,the methods for quantum-dynamical simulation of reactive processes in complex systems in the condensed phase are still under development.We briefly consider electron and proton quantum dynamics. The proper quantumdynamical treatment of fast electronic transfer re- actions and reactions involving electronically excited states is very complex not only because the Born-Oppenheimer approximation brakes down but
16 Berendsen also because a reliable description of potential surfaces for excited states in complex molecules is very difficult,even within the Born-Oppenheimer approximation.With electron transfer processes in proteins one resorts to simplified descriptions,e.g.for electron transfer paths 97,98 or for the protein environment[99].In some cases,such as solvated electrons,the time- dependent electron wave function can be computed in a mixed QM/MD sim- ulation,given a suitably parametrized interaction function between electron and solvent atoms [28,100.The electron wave function is then usually defined on a grid,and its time evolution solved by a split-operator technique using fast Fourier transforms.In these adiabatic simulations the electron remains in its Born-Oppenheimer ground state,although it is possible to compute excited states (and thus spectra)as well.Non-adiabatic transitions between different states,driven by solvent fluctuations,are generally not important for such systems in equilibrium at normal temperatures. For the case of intramolecular energy transfer from excited vibrational states,a mixed quantum-classical treatment was given by Gerber et al.al- ready in 1982 [101].These authors used a time-dependent self-consistent field (TDSCF)approximation.In the classical limit of TDSCF averages over wave functions are replaced by averages over bundles of trajectories,each obtained by SCF methods. Proton transfer over hydrogen bonds,which is a rate-determining step in many enzyme reactions,is a case where both the quantum character of the particle (implying proton tunneling)and the non-adiabatic transitions between pure states are important [102].A quantum-dynamical treatment of the proton(s)requires embedding of one or more quantum degrees of freedom (d.o.f.)in a classical dynamic environment (yielding QD/MD)and there are two main points of discussion concerning the embedding procedure. One of the discussion points is how the quantum system reacts back on the classical d.o.f.,i.e.,how the forces on the classical system should be derived from the quantum system.One can use the gradient of the effective energy, i.e.,of the expectation value of the total energy 0 F:=RH(R)) where Ri are the coordinates of the i-th classical particle and is the wave function which is parametrically dependent on R.These forces are considered to be the 'exact'forces [103].Alternatively,the Hellmann-Feynman force,i.e., the expectation value of the gradient of the energy aH(BΦ. F:=-0R: can be used.These expressions are only equivalent if is a pure eigenstate. In other cases,extended Hellmann-Feynman forces [104 have been derived from the R-dependence of However,if represents a mired quanturn state
MD:The Limits and Beyond 17 as correct solution to the time-dependent Schrodinger equation (see below), the correct force is [105,106] (E1 F:=-R and not the gradient of the effective potential.This is imperative to produce conservation of momentum over the whole system as can be seen immediately from the Ehrenfest relation [107],which states that the rate of change of the expectation value of the momentum of a quantum particle is av(r). d(p)=r dt where r is the quantum coordinate.Since for a potential energy term between two particles of the form V(r-R)the gradients with respect to r and R are equal and opposite,the total momentum change will be zero.A proof that both energy and momentum are conserved has been given in [106].The correctness of this force has also been shown by Bornemann et al.108 who derived the equations of motion when one particle of two interacting quantum particles is systematically taken to the classical limit. The reason that non-adiabatic transitions must be included for protons is that fluctuations in the potential for the quantum degrees of freedom due to the environment (e.g.solvent)contain frequencies comparable to the transi- tion frequencies between protonic quantum states.In such cases pure quan- tum states do not persist. The second discussion point is how the actual quantum system is to be described:should one follow the time evolution of the time-dependent Schrodinger equation(TDSE)that allows mixed states to evolve,or should one insist on selecting a pure state,taking care of(sudden)transitions be- tween states by some additional action in order to satisfy the time evolution of probabilities of states as dictated by the TDSE?The former approach was fol- lowed,among others,by Bala et al.in wave packet dynamics applied to proton transfer in phospholipase A2 [109,110 and by us in the Density Matrix Evolu- tion(DME)method which describes the mixed time-dependent wave function on a simple,appropriately chosen,basis set.[105,111,112,113,114,106,115]. DME is obviously not capable of giving a correct response of the classical en- vironment to quantum transitions,but is perfectly able to describe initial rate processes or quantum systems that only weakly influence their environ- ment.In fact,DME is the common method used in the evolution of nuclear spin magnetization [116].The latter approach has led to the surface hop- ping method pioneered by Pechukas [117,with a modern formulation by Tully [118.The basic idea is that the dynamics of a pure quantum state is followed,simultaneous with the classical dynamics of the environment.At ev- ery step the probability of a transition to another quantum state is calculated and such transitions are realized on a stochastic basis.When a transition is made,velocities are scaled to conserve total energy.The method has been