8 Berendsen These numbers are not negligible.At present such effects are on the average compensated by other force field terms through empirical parametrization. We may conclude that the matter of optimal algorithms for integrating Newton's equations of motion is now nearly settled;however,their optimal and prudent use [28]has not been fully exploited yet by most programs and may still give us an improvement by a factor 3 to 5. 2.2 Force Fields While simulations reach into larger time spans,the inaccuracies of force fields become more apparent:on the one hand properties based on free energies, which were never used for parametrization,are computed more accurately and discrepancies show up;on the other hand longer simulations,particularly of proteins,show more subtle discrepancies that only appear after nanosec- onds.Thus force fields are under constant revision as far as their parameters are concerned,and this process will continue.Unfortunately the form of the potentials is hardly considered and the refinement leads to an increasing num- ber of distinct atom types with a proliferating number of parameters and a severe detoriation of transferability.The increased use of quantum mechanics to derive potentials will not really improve this situation:ab initio quantum mechanics is not reliable enough on the level of kT,and on-the-fly use of quantum methods to derive forces,as in the Car-Parrinello method,is not likely to be applicable to very large systems in the foreseeable future. This situation,despite the fact that reliability is increasing,is very unde- sirable.A considerable effort will be needed to revise the shape of the poten- tial functions such that transferability is greatly enhanced and the number of atom types can be reduced.After all,there is only one type of carbon; it has mass 12 and charge 6 and that is all that matters.What is obviously most needed is to incorporate essential many-body interactions in a proper way.In all present non-polarisable force fields many-body interactions are incorporated in an average way into pair-additive terms.In general,errors in one term are compensated by parameter adjustments in other terms,and the resulting force field is only valid for a limited range of environments. A useful force field should be accurate and simple.Therefore it is desirable that polarisability be incorporated by changing charges(positions or magni- tudes)rather than by incorporating induced dipole moments,which involve dipole field gradient tensors to be computed.The best candidate is the shell model,which represents electron clouds by charges on a spring;a detailed study of nitrogen in all its phases by Jordan [50]could serve as a guide. The task of generating a new general force field with proper many-body in- teractions comprises a complete overhaul of the present force fields and a completely new parametrization involving not only static data but also free energy evaluations.This non-rewarding task is not likely to be accomplished without an internationally organized concerted effort
MD:The Limits and Beyond 2.3 Long-Range Interactions The proper treatment of long-range interactions has not yet been settled in a quite satisfactory way.The use of a cut-off range for dispersion interactions with r-6 radial dependence does not present a severe problem,although continuum corrections for the range beyond the cut-off are often necessary, particularly for pressure calculations.But a simple cutoff for Coulombic terms can give disastreus effects,specially when ionic species are present or when dielectric properties are required.It has been observed(see e.g.[51])that in electrolyte solutions the radial distribution of like ions strongly peaks around the cutoff radius,while also the short-range structure is severely distorted In dipolar systems without explicit charges a cut-off is tolerable to compute structural and energetic properties,as long as dipoles are not broken in the cut-off treatment,but dielectric properties cannot be evaluated with any precision [52.Another effect of(sharp)cut-offs,among other artefacts,is the introduction of additional noise into the system.Sufficiently smooth cut-offs, on the other hand,severely deviate from the correct Coulomb interaction Several structural and dynamic artefacts have been described,see e.g.[45, 54,55,56,57,58,59,60,61,62].Therefore it is recommended that at least some method to incorporate the long-range part of electrostatic interactions be included. If the simulated system uses periodic boundary conditions,the logical long-range interaction includes a lattice sum over all particles with all their images.Apart from some obvious and resolvable corrections for self-energy and for image interaction between excluded pairs,the question has been raised if one really wishes to enhance the effect of the artificial boundary conditions by including lattice sums.The effect of the periodic conditions should at least be evaluated by simulation with different box sizes or by continuum corrections,if applicable(see below). A survey of available methods has,among others,been given by Smith and Van Gunsteren [51;see also Gilson [63.Here a short overview of methods, with some comment on their quality and properties,will be given. Cut-off Methods The use of an abrupt potential cut-off radius re for the evaluation of the (electrostatic)potential,while feasible for Monte Carlo simulations,implies a delta-function in the force at re.If a delta function is really implemented,it gives an unphysical force when particle distances hit the cutoff radius,and a dynamic behaviour that is very dependent on the cut-off range [64].The use of an abrupt force cut-off is compu- tationally more acceptable,although it introduces additional noise.For systems containing dipolar groups that are represented by charges it is mandatory that the cut-off is based on charge groups rather than on in- dividual charges [65].One should realize,however,that a force cut-off implies a shift in the potential,since the latter is continuous at re and zero beyond re.Therefore Monte Carlo and MD simulations with cut-offs are not expected to give the same equilibrium results.A useful extension
10 Berendsen is the use of a twin-range cut-off,where the forces from a second shell are evaluated every 10 to 20 steps simultaneously with the construction of a neighbour list for the first shell.Smooth cut-offs can be accomplished by switching functions applied to the potential.They generally cause a better behaviour of the integration algorithm,but also disguise the errors that are nevertheless made.Ding et al.[66]have shown that traditional switching functions cause large errors in structure and fluctuations when applied to a dendrimeric polypeptide;they suggest a large smoothing range between 0.4289re and re.Fun Lau et al.[61]found structural and dynamical influences of switching functions on water,and Perera et al. [62]found severe influences on the dynamic properties of ions in aque- ous solution.All cut-off methods suffer from severe distortions in systems containing full charges and deny the evaluation of dielectric properties The latter is due to the fact that fluctuations of the total dipole moment in a sphere are much reduced when the sphere is surrounded by vacuum [52. Reaction Field Methods In order to alleviate the quenching of dipole mo- ment fluctuations,a reaction field due to a polarizable environment be- yond the cut-off can be incorporated in conjunction with cut-off meth- ods.The reaction field [52,51]is proportional to the total dipole moment within the cut-off sphere and depends on the dielectric constant and ionic strength of the environment.On the same level,a reaction poten- tial (Born potential)should be applied,proportional to the total charge in the cut-off sphere.This is only applicable for homogeneous fluids,and therefore not generally useful.However,in the case the dielectric constant or ionic strength of the environment is taken to be infinite,conducting or tin-foil boundary conditions arise with simple expressions for the forces. Such reaction fields and potentials are of course in general also incorrect, but they produce well-behaved forces and allow better subsequent correc- tions based on continuum theory (see next item),especially in systems like macromolecules in aqueous solution.The expression for the electric field at particle ri is E(r)=(4ro)-1∑4r(r3-r。3方 (rij<re), which amounts to a well-behaved shifted force.By integration the total electrical interaction energy becomes: Vel 1 8πe0 ∑99 2re (r访<rc). This is close,but not equal to the tin-foil Born-corrected energy VBorn= 1 8πe0 i,≠
MD:The Limits and Beyond 11 The discrepancy is not large and the last term is zero for a system with- out net charge.Thus we see that the use of a shifted Coulomb force is equivalent to a tin-foil reaction field and almost equivalent to a tin-foil Born condition. Continuum Corrections If the geometry of the simulated system is not too complex,it is possible to make corrections for the 'incorrect'long- range treatment,based on continuum considerations.This has been con- vincingly shown by Wood (67]in a paper that has not received the at- tention it should.The idea is that the computational world'has its own physics (like cut-offs and periodic boundary conditions),and that the differences with the 'real world'are fairly smooth and therefore can be treated by continuum methods.Such corrections were made on earlier simulations by Straatsma (68]on the free energy of ionic hydration,using various cut-offs for the ion-water and the water-water interactions.While in the original paper the usual Born correction was made,a discrepancy remained due to the neglect of water-water interactions between pairs of molecules that are both correlated with the ion.Wood showed that all results fitted beautifully after correction based on the spatial distribution of solvent polarization.Such corrections could also be made had a tin-foil reaction field been used. The same idea was actually exploited by Neumann in several papers on dielectric properties (52,69,70].Using a tin-foil reaction field the relation between the (frequency-dependent)relative dielectric constant e(w)and the autocorrelation function of the total dipole moment M(t)becomes particularly simple: o)-1=T(盖MoMe e-iwt dt Hierarchical Methods Methods that group more particles together for in- creasing distances [32,72,73]scale roughly proportional to the number of particles N or to Nlog N,rather than to N2(as for straightforward summation over pairs).For large system sizes such linear hierarchical methods should win out over other methods.Hierarchical methods,rou- tinely applied in the simulation of star clusters and galaxies,have also been adopted for proteins [74]and implemented in simulation programs for large molecular clusters [75].These methods have been extensively compared to each other and to Ewald summation 76,with the result that they only surpass Ewald summation for particle numbers in the hundred thousands.Since it has been known for a long time that Ewald summation is considerably more expensive than a Poisson solver on a grid [11,77](see next item),I conclude that there is not much point in pursuing these methods for molecular simulation. Separation of Short-and Long-Range;Ewald and PPPM Methods If we split the total Coulomb interaction in a short-and a long-range contribution,chosen to be smooth functions of the distance,the two
12 Berendsen contributions can be more efficiently handled by different techniques [65].The short-range contribution is evaluated for each pair on the basis of a pair list.The long-range part can be recast in terms of a Poisson problem and then solved by an appropriate Poisson solver.The choice of Poisson solver depends on whether the system is periodic or not,on the preconditioning (i.e.,is an approximate solution available from the previous step?)and on the wish to implement the algorithm on a parallel computer.It is possible to update the long-range part less frequently in a multiple time-step algorithm.The popular Ewald summation 10,78,79 and its variants [80],the efficient Ewald-mesh technique [40],and the particle-particle particle-mesh (PPPM)method [11]are all special cases of this class of techniques. Consider2 a short-range interaction defined by the electric field at position T E(ri)=Anco The function f(r)is a force-switching function that goes smoothly from 1 ar r=0 to 0 at r =re.The long-range part of the field,i.e.,what remains from the complete Coulomb field: B=∑a-f暗 can be considered to be generated by a charge density p(r): p(r)=oV·E=∑99(r-r), with o(r)-gr) 4nr2 dr Thus the long-range field(and potential)is generated by a charge density which is the convolution of the charges in the system with a radial charge spread function g(r)dictated by the short-range force-switching function f(r).The task is to solve for the long-range potential (r)(the negative gradient of which is the long-range field)from the Poisson equation V26=-p/eo: If a gaussian function is chosen for the charge spread function,and the Poisson equation is solved by Fourier transformation (valid for periodic 2 The force functionf(rdiffers from that in ref.(65]by a factor,yielding simpler expressions.Some errors in that reference have been corrected