dynamicsinsystemswithmultipletimescales:Systemswithstiffandsoftdegreesoffreedomandwithshortandlongrangeforces MarkE.Tuckerman”andBruceJ.BerneDepartmentofChemistry,ColumbiaUniversity,NewYork,NewYorkIO027 (Received9July1991;epted12August1991) Whentherearehighandlowfrequencymotionsinsystemswithlongandshortrangeforcesajudiciouschoiceofreferencesystemleadstoverylargeelerationsinmoleculardynamic(MD)simulations.Buildingonourpreviousworkwherewehavedevelopedreferencesystemmethodsforsystemswithhighfrequencyoscillators,disparatemasses,orlongrangeforces,wepresentadoublereferencesystemmethodwhichleadstoelerationofasmuchas20insystemsconsistingof864moleculeswithLeonard-Jones(12-6)forces.Muchlargersavingsshouldbeachievedwhenthismethodisappliedtolongerrangeforcesandlargersystems.
1.INTRODUCTION Inourrecentworkwehaveaddressedtheproblemofmultipletimescalesinmoleculardynamics(MD)simulations.1,53Themostobviouscasesinvolveeither“stiffoscillators”dissolvedin“softfluids”’(e.g.,N,inAr)orlowmasssolutesdissolvedinhighmasssolvents2(e.g.,HeinXe).Insuchcases,standardMDmethodsrequireveryshortintegrationtimesteps,St,toguaranteestableintegrationofthefastdegreesoffreedom.Thus,tofollowtherelaxationoftheslowdegreesoffreedomrequiresthegenerationofaverylargenumberofinteractions. Wehaveinventedreferencesystemmethodssuchasnumericalanalyticalpropagatoralgorithm(NAPA)’andreferencesystempropogatoralgorithm(RBSPA)273forintegratingsuchsystems.Whenthesystemsarediluteinthefastdegreesoffreedomsthesemethodsgreatlyeleratethesimulationslargelybecauseforcesfortheslowcoordinatesarerecalculatedmuchlessfrequentlythaninthestandardalgorithms.TheelerationfactorforRESPAoverstraightforwardvelocityVerletisdefinedtobe tq(&)srverlet(St) .(l.l) TRES(PSAt)’ whereTverl(eStt)isthecputimerequiredtocarryouta directsimulationusingthestraightforwardvelocityVerlet integrator4withatimestepStwhileTREsPA(St)isthecpu timerequiredbyRESPA(orNAPA)tosimulatethesame amountofrealtime,integratingthereferencesystemwitha timestepSt.Wehaveachieveduptoeightfoldelerations insystemswehavesimulatedtodate. intermolecularforcescanbesubdividedintoshortand longrangeparts.Theshortariesmuch morerapidlythanthelongponent,andthusde- terminestheintegrationtimestep.Toolargeatimestepwill leadtoanerrorinthenewpositionswithaitantly largeerrorintheshortrangeforcewhichsubsequentlyleads toamplifyingerrors.Eveninsystemsinwhichthereisno obviousseparationoftimescalesthisleadstoadefactomul- tipletimescaleproblem.Oneisforcedtousesmallintegra- tiontimestepsandtorecalculatethefullforceafterevery smalltimestep.Judiciouschoiceofareferencesystem(RE- a)Ph.D.studentintheDepartmentofPhysics,ColumbiaUniversity. SPA)hasallowedustoreducetheseputationsconsiderablyandtotherebyachievefactorsaslargeas6intheelerationofthesimulations3fortheuracyasmeasuredbytheenergyconservation,whichismeasuredby &L$Ni=I1E,0BE1)o(1.2) whereNisthetotalnumberofMDsteps,Elistheenergyatstepi,andE.istheinitialenergyofthesystem. Inthispaperweshow,forsystemsinwhichtherearemultipletimescales(stiffandsoftdegreesoffreedom)andinwhichtheforcescanbesubdividedintoshortandponents,thatadoubleapplicationofRESPAleadstoverylargeelerationsofthesimulationtimeforsimplesystems.Forsimplicitywetreatthetwosystemsalreadypresentedinourpreviouspapers:namely(a)onestiffdiatomic dissolvedin862Lennard-Jones(LJ)(12-6)atoms’and(b)40lightLJparticles(m=1)dissolvedinafluidconsistingof824heavyLJatoms(m=100)allinteractingwiththesameLJ(12-6)potential.2DoubleRESPAsolutionsyieldtwentyfoldelerationsinthecputimesrequiredforthesimulations. WeexpectthatforlargersystemsandsystemswithlongrangeforcesdoubleRBSPAcanyieldasmuchasafortyfoldeleration. II.SUBDIVISIONOFFORCES Thestartingpointbiningreferencesystemsintheseparationoftheinteratomicforcesintoshortandlongrangecomponentsaccordingto F(x)=
F,(x)$4(X). (2.1) TheseparationcanbeachievedbymeansofaWCAsubdivisionorbyusingaswitchingfunction.3WehaveshownthatbychoosingareferencesystembasedonF,(x),factorsofbetween2and4savingsincputimecanbeachievedforsimpleLJ(12-6)systems,andfactorsexceeding5canbereachedfor(12-1)systemswhenEwaldsummationisused.3TheRESPAalgorithmisimplementedbywritingthetrajectoryx(t)asasumofareferencetrajectoryX,(t)andacorrectionxl(t)whichsatisfytheequationsofmotion 8362
J.Chem.Phys.95(1l),1December19910021-9606/91/238362-03$03.00 @1991AmericanInstituteofPhysics
M.E.TuckermanandB.J.Berne:Moleculardynamicswithmultipletimescales 8363 mjz=,
F,(x,)+F,w, (2.2) Mr=K(xsf+)-
F,(x,)+
F,(x,+x,) -f’,(O), C-J.31 whereFl(0)isthevalueofthelongrangeforceatthebegin- ningofthetimestep.Theinitialconditionsaretakentobe x,(O)=x(O),
Z,(O)=wa, (2.4) Xl(O)=&(O)=o. (2.5) Equation(2.2)isintegratedfornlittletimestepsetsubject totheinitialconditions,Eq.(2.4)andthenthecorrectionxl putedusingabigtimestepAt=nStsubjecttothe initialconditions,Eq.(2.5).Thetruetrajectoryisthengiven byx(At)=x,(At)+xl(At).Theinitialconditionsarere- setsothat x,(O)=x(At),&(O)=k(At), (2.6) x,(O)=&(O)=
0, (2.7) andtheprocedureisrepeatedforeachstepdeterminingtheinitialconditionsforthenextstep. Ill.DISPARATEMASSSYSTEMS Indilutefluidmixturesconsistingoflightsoluteatomsandheavysolventatoms,thereisaseparationoftimescales:thelightatomsmovemuchmorequicklythanthesolventatoms.Asimpleexampleofthisisthatof864LJ(12-6)atomsconsistingof40lightsphereswithm=1and824sphereswithm=100allwiththesamediameter0andthesamewelldepthE.Inthissystem,thetimescalefortheheavyparticleis10timeslongerthanthatforthelightparticles.WehavealreadyshownhowthissystemcanbeintegratedusingRESPAresultinginasevenfoldspeedupofthesimulationoverstandardresultsusingthevelocityVerletintegrator.2HerewewishhiswithanotherversionofRESPAbasedonthesubdivisionoftheforcesintolongandshortponents. Considerasystemconsistingofamixtureof40lightLJparticleswithmL-1and824heavyLJparticleswithM=100.Denotingthesetoflightparticlecoordinatesasxandtheheavyparticlecoordinatesasy,theequationsofmotiontaketheform g=$%Y), j+&y). (3.1) Asbefore,theforcesandcoordinatesarebrokenupintoshortandlongponentsgivingequationsforthereferencesystemtrajectoriesandcorrections %=;[Fx.s(xs,ys)+Fx,(O)], (3.2) j;,=$[J--(&
Y,)+r;;,CO>], (3.3) X*l-=+ [
F,(x, ++
Y, SY,) -F.sh~s)] +; [F&s +x,,ys +~r) -MO)], (3.4) j;,=a[
F,(x,+w+
Y,)-F&,~
Y,)] +m[FJxf,x,,y+,
Y,-)-Fy#VJ.(3.5) Theinitialconditionsonx,,ys,xI,andy,aregivenbyEqs.(2.4)and(2.5).Tohandlethemassdisparity,wemakeafurthersubdivisionofx,andxzordingto x,=XZO+js,;XI=xi”’+SI, (3.6) whereXL”’andS,arechosentosatisfytheequationsofmotion m$0’_s--!
-@‘xs[x:“‘,Js] +
F,,(O)h (3.7) ;i‘s=i-@,s[xl’)+&,ys]-
F,,[x6”‘,A]
3, (3.8) wherejjXindicatesthattheheavyparticlereferencesystemisheldfixedwhilexl”’isintegrated.Asimplechoiceistofix,atitsinitialvalue.TheinitialconditionsaretakentobethesameasinEqs.(4.5)and(4.10).Similarly,XI(‘)andS,satisfytheequationsofmotion 2;”=f{
F,[xs+x:“,
Y,+PI]-F.s(x,,y>s3 +;{F&, +x;“,y,-I+Y-%(O)), (3.9)
8,=;@‘xs[xs+x:0)-I-
4,~s+YZ] -
F,[xs+x:“,
Y,+R]
3 ++{Fx[,xs+xl”’+
4,~s+YZ] -F&s+x~“,Y+,P1]
3. (3.10) Becauseoftheinitialconditionsonxl,xj”(O),2;‘)(O),
S,
(0),andb,(0)areall0.TheprocedureistointegrateEq. (3.7)forn,timestepsSt,andthentocorrectordingto Eq.(3.8)whilesimultaneouslyintegratingEqs.(3.9)and (3.3)withatimestepStz.Theinitialconditionsareresetandtheprocedureisrepeatedn2times.Finallythecorrec- tions6,andy,putedfromEqs.(3.10)and(3.5), respectivelyusingabigtimestepAt-n,St,=II,n,St,.WehavealreadyshownhowtoadoptthevelocityVerletalgo- rithmforusewithRESPA,3andtheextensiontodouble RESPAfollowsstraightforwardfromthistreatment.Letvibetheelerationfromthelightparticlereferencesystems andv2betheelerationfromtheshortrangereference system.Weexpecttheoverallelerationtobe17Iq2.Thetestiscarriedoutontwosystems.Oneisattem- perature0.67anddensity0.86whichcorrespo2dstothetri- plepoint.TheenergyconservationissetatAE=2Xlo-6whichrequirestheVerlettimesteptobechosenas St,=2x10-
3.ForRESPA,wechoosen,=10and
J.Chem.Phys.,Vol.95,No.11,lDecember1991 8364
M.E.TuckermanandB.J.Berne:Moleculardynamicswithmultipletimescales n,=
5.Thisgivesanoverallcpusavingfactorof77=13.Directmeasurementsof77,andv2yieldvaluesof7and2,respectivelysothatthepredictionisv1Q=14.Thesecondsystemisattemperature1.0anddensity0.86.KeepingthesyeVerlettimestepgivesanenergyconservationofAE=3X10-
6.ForRESPAwechoosey1i=7andnZ=6whichgivesanoverallsavingof7=20.ThedirectmeasurementofqZgives3.3sothatthepredictionisviqZ=23. IV.ASTIFFOSCILLATORDISSOLVEDINASOFTFLUID Considerafluidmixtureconsistingofoneverystiffdiatomicmoleculedissolvedin862LJatomswherethemolecularatomsinteractwiththesolventatomsandthesolventatomsinteractwitheachotherthroughthesameLJ(12-6)potential.WehavealreadyshownhowNAPAandRESPAcanbeusedtotreatthemultipletimescaleproblem.’Thesestudiesyielduptoeightfoldelerationsoverstandardmethodsforthesameenergyconservation.Hereweshowhowtoimproveuponthisbyincludingthebreakupoftheforceintoshortandlongponents. Theequationofmotionfortherelativecoordinateroftheoscillatortakestheform pi:=./w+F(r), (4.1) wherepisthereducedmass,f(r)istheoscillatoryforce,andP(r)istheforceduetothesurroundingsolventatoms.Wehaveshownthatchoosingareferencesystembasedsolelyonf(r)canleadtoafactorof8elerationincputimewhenafrequencyof300isusedfortheoscillator. However,ifwebinethisoscillatoryreferencesystemwithshortrangeforcereferencesystem,asubstantialimprovementcanbeachieved.Allthesolvent-solventandsolvent-soluteforcesaresubdividedordingtoEq.(2.1).Theequationofmotionfortheoscillatornowtakestheform pi”=.m+E(r)+E;lcr>. (4.2) Therelativecoordinateriswrittenasthesumofareferencesystemtrajectoryr,andacorrectionr,.r,andr,satisfytheequationsofmotion pys=f(r,)+I;:(<)+Fl
(0), (4.3) pi;l=f(r,+rr>-f(r,>+
F,(r,+rl’I) -
F,(r,)+4(r,-kc)---
F,(O), (4.4) whereF,(0)denotesthevalueofthelongrangepartoftheforceatthebeginningofatimestep.Theinitialconditionsarechosentobe r,
(0)=r(O),fs
(0)=L(O), (4.5) 21(O)=f(O)=
0. (4.6) Thereferencesystemtrajectoryr,isfurthersubdividedordingtor,(t)=ri”)(t)+
S,(t),whereri”’andS,satisfy pF2°’=f[ri”], (4.7) pZs=f[rZ”+Ss]-f[pj’)
J -t-
F,[r~~‘+4] --
F,(O) (4.8) withinitialconditions rJO’(O)=r,
(0), i$O’(O)=i;(O),
S,(O)==&(O)=o. (4.9)(4.10) TheprocedureistointegrateEq.(4.7)forn,timestepsSt,subjecttotheinitialconditions,Eq.(4.9)andthensimultaneouslytocorrectforri”’usingEq.(4.8)andevolvethereferencesystem,Eq.(2.2)forthesolventatomsusingatimestepSt,.Theinitialconditionsonry)arereset,andtheprocedureisrepeatedn2timestogeneratethefullreferencesystemtrajectoryr,(t)andx,(t)fortheoscillatorandsolventatoms,respectively.ThenthecorrectionsputedordingtoEqs.(4.4)and(2.3)foronebigtimestepAt=n,St,=n,n,St,.Ifv1isthecpuelerationfactorfortheoscillatorreferencesystemaloneand7,isthecpuelerationfactorfortheshortrangeforcereferencesystemalone,thenweexpecttheoverallspeeduptobetheproduct 771%.Wehavetestedthispredictiononanoscillatorforwhich f(r)=-,ud(r-a)withp=1,w=300,anda=1.25in abathof864LJatomsattemperatEi-e1.0anddensity0.9.TheenergyconservationissetatAE=2x10-5whichrequiresaVerlettimestepof2.5X10v4.IntheRESPAsimulation,weuse‘
2,=8andn2=6andalargetimestepAt=1.39X10-*.Theseparametersgiveacpusavingfactorof7=22.Frompreviouswork,wehavedeterminedthat-fr=7.9whileQ=3.2whichgivesapredictionofthesavingfactorofv1Q=25incloseagreementwithourfinding.
V.CONCLUSION Thereferencesystemmethods(RESPA)leadtoadramaticelerationofmoleculardynamicsforsystemswithmultipletimescalesandshortandlongrangeforces.Thesemethodsaresimpletousearidarecapableofgeneralizationtoplicatedsystems.TheunderlyingequationsofmotionusedinRESPAareexactandcanbesolvedusinganyofthestandardnumericalintegrators.TheworkpresentedhereusesthevelocityVerletintegrator,”butwearepresentlytryingtoapplyRFSPAtodynamicalsystemswithbondlengthandbondangleconstraintsusingSHAKE.’Wearealsopresentlytryingtoextendthesemethodstotreatlargemoleculeswithmanystiffcoupledinternaldegreesoffreedom. ACKNOWLEDGMENT ThisworkwassupportedbyagrantfromtheNationalScienceFoundationandfromthePetroleumResearchFund,administeredbytheAmericanChemicalSociety, ’
M.Tuckerman,
G.Martyna,andB.J.Berne.J.Chem.Phys.(1990). *
M.Tuckerman,
B.J.Berne,andA.Rossi,
J.Chem.Phys.94,14653M.Tuckerman,
G.Martyna,andB.J.Berne,
J.Chem.Phys. (1991).’
W.CSwope,
H.C.Andersen,
P.H.Berens,andK.R.Wilson. Phys.72,435O(1980).‘
J.P.Ryckaert,G.otti,andH.J.C.Berendson,
J.Comput.327(1977). 93,1287 (1991).94,6811
J.Chem.Phys.23,
J.Chem.Phys.,Vol.95,No.11,iDecember1991