wwwelseviernlTabulatedpotentialsinmoleculardynamicssimulationsDWolffWGRuddDepartmentofPhysicsOregonStateUniversityCorvallisOR97331USA AbstractTheuseoftablestodetermineinterparticlepotential ID: 189120
Download Pdf The PPT/PDF document "ComputerPhysicsCommunications120" is the property of its rightful owner. Permission is granted to download and print the materials on this web site for personal, non-commercial use only, and to display it on your personal computer provided you do not modify the materials and that you retain all copyright notices contained in the materials. By downloading content from our website, you accept the terms of this agreement.
ComputerPhysicsCommunications120 www.elsevier.nlTabulatedpotentialsinmoleculardynamicssimulationsD.Wolff,W.G.RuddDepartmentofPhysics,OregonStateUniversity,Corvallis,OR97331,USA AbstractTheuseoftablestodetermineinter-particlepotentialsandforcescanspeedupdeterminationsofthesefunctionsinclassicalmoleculardynamicssimulationsbyasmuchasafactoroffourfortypicalpotentialfunctions.Doingsoresultsinalossofaccuracyintheresultingtrajectoriesandrequiresincreaseduseofmainmemorytostorethetables.Ananalysisof 1.IntroductionMoleculardynamicscs113]isaprovenmeansforsimulatingsystemsofatoms,molecules,andotherparticles.Themethodinvolvesintegratingtheequationsofmotionderivedfrominter-particlepotentialfunctions.Statisticalaveragesthenyield E-mail:wolff@physics.orst.eduE-mail:rudd@cs.orst.eduforthesequantitiesWefocusonpotentialsthatincludetermsthatcanbecomputedtoanydesiredaccuracyaspartoftheini-tiationbeforeactualsimulationsbegin,withempha- Henceforth,wewilluseªpotentialºtorepresentpotentialenergy, D.Wolff,W.G.Rudd/ComputerPhysicsCommunications120(1999)20latedpotentialshavealsobeenusedtoacceleratethecalculationoflongrangeelectrostaticforcesintheEwaldsumum9].Memorycycletimesoncommonworkstationsareabout60ns,whichisapproximatelythetimerequiredforasingledouble-precision¯oating-pointcal-culation.Wethereforewouldexpectthatlookingupvaluesforpotentialswouldbefasterthancomputingthemforallbutthemosttrivialformsofthepotential.However,thisisnotalwaysthecase.ProcessorsusedinseriousMDcalculationsincludeinternalpipelinedprocessorsfor¯oating-pointarithmetic.Manyoffervectorprocessorsaswell.Boththeseprocessorenhancementsrequireasteady¯owofoperandsfrommemory.Ef®cientcachinghardwareandalgorithmsandoptimizingcompilersaretunedtopreventdelaysduetointerruptionsinthemovementofdatafrommemory.Thesetechniquesworkbestwhenthedataaccesspatternsarelocalizedinmemoryandoccurinapredictableorder.Ontheotherhand,usingtabulatedpotentialsisequivalenttoperformingrandomaccessesintolargearrays,sothataccessesareneitherlocalizednorpredictable.Wethereforeexpectandobserveatrade-offbetweenusingtablestosavearithmeticcalculationsanddegradationinrawprocessorperformance.Nevertheless,weexpecttobene®tfromtabulatingpotentialsforthosethatrequirerelativelyintensearithmeticcomputationandrelativelylittleassociatedlogic.Thisstudysupportsthatexpectation.Evenwhentabulatingpotentialsoffersbetterper-formancethancomputingthem,therearestilltwoothertrade-offstoconsiderindecidingwhichmethodtouse.Theuseoftablesguaranteesalossofaccu-racyincomputingthepotentials.WeuseDParith-meticinMDsimulationstoeliminatepossibleinac-curaciesthatcouldcloudtheinterpretationofresults.Itisnotclearthatweneedthatlevelofaccuracyinmanysimulations.Oftenwedonotknowvaluesofparametersinpotentialfunctionswithanydegreeofcertainty;inmuchoftheworkinmaterialsscienceandotherapplicationareasweusesemi-empiricalpoten-tialsthatinvolveconstantsthatareknowntoonlyoneortwodecimalplacesinaccuracy.Thesevaluesfre-quentlyappearinpower-laworexponentialtermsinpotentialfunctions.Therefore,inaccuraciesindeter-miningvaluesofpotentialsunderthesecircumstancesareswampedbylackofprecisioninourknowledgeoftheexactpotential.Notethatthiscommentdoesnotapplytonumericaltechniquesusedtointegratetheequationsofmotion;thesetechniquesmustbechosencarefullytoensurethattheyaresuf®cientlyaccurate.Otherwise,numericalinstabilitiescanarisethatleadtoerroneousordivergentsolutions.Thethirdtrade-offbetweenusingtablesandcom-putingpotentialsandforcesisthespace-timetrade-off,whichisalsolinkedwiththeaccuracyissue.Tomakethecomputationsfaster,wemustusememorytostoretablesthatwouldotherwisebeavailableforotherpurposes.Thisfactorbecomesparticularlysig-ni®cantwhenweconsidersimulatingsystemswithmorethanonespeciesofatom.Fortwobodypoten-tials,thenumberoftablesincreaseswiththenumberofpossiblepairinteractions,tableswhereisthenumberofatomicspecies.Whenthecalculationsaretobedoneonparallelordistributedprocessors,wemustkeepcopiesofthetablesoneachprocessortoavoidextrainter-processorcommunica-tions.Memoryisanimportantfactorinlimitingthesizesofsimulationwecando.Asmainmemorysizescontinuetogrow,therelativeoverheadfromthetablesbecomessmaller.Ourresultsdescribedbelowindicatethatundernormalcircumstancestheuseofmemoryfortablesisnotanimportantconsideration.Inthisreport,weexplorealloftheseissues.Webeginwithadiscussionofhowonemightconstructandaccesstablesef®ciently.Wethenpresentourre-sultsontheanalysisoferrorsresultingfromtheuseoftablesandtheireffectsontheresultsofsimulations.Followingacomparisonofperformancebetweenthetwomethodsforseveraldifferentkindsofpotentials,wepresentourconclusionsandrecommendations.2.ImplementingtablesWell-designedmoleculardynamicscodescodes2,5,6,10112]includesinglemodulesthatcomputepotentials.Tochangeacodetoobtainpotentialsandforcesfromtables,wesimplyreplacethecodethatcomputesthesevalueswithcodethatretrievesthemfromdatastructuresinwhichtheyhavebeenstoredbeforethesimulationsarestarted.Inorderforthisto D.Wolff,W.G.Rudd/ComputerPhysicsCommunications120(1999)20beworthwhile,anef®cientmethodofindexingthetablemust®rstbefound.2.1.HashfunctionsTheobjectiveistodeterminethecontributionstotheforceandpotentialactingonparticlebyalltheotherparticles.Weassumetheseforcesdependonlyonseparationsbetweentheparticles.ThepotentialforceroutinesinMDsystemsusuallyin-cludecheckstoseeifparticlesarewithinacertaindistancefromparticle,asforexample,whenthepo-tentialfunctionincludesashort-rangecutoff.Thesecomparisonsaredoneusingtoavoidcomputingasquareroot.Therefore,itmakessensetousethemeasureofthedistanceusedtodeterminewhichelementsofthetabulatedforcesandpotentialstouse.Sincepointersandindexesforarraysareintegers,weneedtoconstructmappingsbetweendoublepreci-sionmeasuresofsquareddistancesandpointersorindexes.Wehaveinvestigatedtwokindsofmappings,castsofscaledvaluesofthedistancemeasureontoin-tegers,andlogicalextractionofindexesfromtheDPdistancemeasures.Considerthepseudo-codebelowforthecomputa-tionofthepotentialandforceinaMDsimulation:for:=computeEnergy:=computeForcetotalE:=totalENotethatthisisforatwo-bodypotential.Forathree-bodypotentialtheloopwouldbeovertripletsijk.ThecomputeEnergyandcomputeForcerou-tinesinpracticearenotseparateroutines,butsincesomeofthecomputationisredundant,theyareoftenjustonefunctionorareinlineddirectlyintotheloop.Thetableversionoftheabovecodemightlooklikethefollowing:for:=energyTableablehash:=forceTableablehashtotalE:=totalEHere,hashisafunctionwhichmapspairwisedis-tancestoindexesinthetables.Anexampleisafunc-tionthatsimplycastsDPvaluestointegers:functionhashdoubletreturnint*10000Thisyieldstheintegerpartoftheargumentaftershift-ingthedecimalpointfourplacestotheright.Wecouldsimulatearoundingoperationwiththefollowinghashfunction:functionhashdoubletreturnintt*10000Theabovetwomappingsaresimple,butrequirecon-vertingaDPnumbertoaninteger,aprocessthatisnotinexpensive.2.2.DoubleprecisionextractionTheothermethodweconsiderfordeterminingin-dexesintopotentialarraysistheextractionofanindexdirectlyfromaDPrepresentationofsomemeasureofdistancebetweenatoms.Theadvantageofthistech-niqueisthatiteliminatesthearithmetic,logic,andimplicitfunctioncallsinvolvedinconvertingDPval-uestointegers.Foragivendoubleprecisionexponent,wecanusethemantissa,truncatedatsomeprecision,asanintegertoaccessthetable.Thetablesizemustbeofsizeatleast2whereisthenumberofbitstakenfromtheleft-handsideofthemantissa.However,inmostcasesweareinterestedinarangeofnumberswhichincludeseveraldifferentpossibledoubleprecisionexponents.Inthiscase,wecanuseafewoftheloworderbitsintheexponentasthe®rstfewhighorderbitsintheintegerindexseeFig.1.Ofcourse,caremustbetakentomakesurethatenoughbitsaretakenfromtheexponenttoensurethattheyareuniqueforallpossibleexponentsthatwouldbeencounteredinasimulation.Wewillassumethatadoubleprecisionnumberis64bitslong.AfunctionforextractionoftheIN- BITSofFig.1mightbethefollowingthecodeiswrittenusingCsyntax,andlinenumbersarein-sertedforfuturereferenceintextractdoubletint*int intindex; D.Wolff,W.G.Rudd/ComputerPhysicsCommunications120(1999)20int int*&t;index=*int p&MASKindex=index20-MANTISSA BITSreturnindexWehaveshownextractasafunctionand,forclarity,wehavewrittenthecomputationsinseveralsteps.Inpractice,itshouldbein-linedwithinthecodeinordertoavoidtheoverheadofthefunctioncall.Intheextractfunction,MANTISSA BITSisthenumberofbitsofthemantissawhichwewouldliketouse.Inotherwords,wearetruncatingthenumberataprecisionof2MANTISSA BITSwhereisthevalueoftheexponent.ThetotallengthofanindextothetableisINDEX BITS,whichisequaltoMAN-TISSA BITSplusthenumberofbitstakenfromtheexponent BITSinFig.1.ThevalueofIN- BITSdeterminesthesizesofthetables,whichshouldbeoflengthatleast2INDEX BITSItisimportanttonotethatthetablesizeisdeter-minedsolelybythelengthoftheintegerusedtoac-cessit.However,allofthetablewillnotnecessarilybeused.Theentiretablewillonlybeusedwhengivenasetofpossibledoubleprecisioninputs,INDEX BITSincludesallpossiblebinarystrings.Clearly,thisde-pendsonanumberoffactorsincludingwhetherornotallpossiblebinarystringsareexpectedinthe BITS,andwhetherornotthefullrangeofman-tissasforagivenexponentisexpected.HerearewhatthenumberedstepsinextractWeestablishapointertoanobjectoftypeintandsetittopointtotheDPargument.Thisal-lowsustoignorethetyperulesthatordinarilywouldrestrictustousingDPoperationsontheargument.Sinceanintegeristypically32bitsandadoubleisusually64,wearerestrictedtotheleftmost32bitsofthedouble.Withanexpo-nentof11bits,thereare20bitsofthemantissatoworkwith.Onecouldusetwointegerpointerstopointtothe®rstandsecondhalfofthedoubleinordertousemorebitsofthemantissa.How-ever,thiswouldrequiresomeextrainstructionstocombinethetwo.MASKisusedinanoperationtoremovetheleft-handbitsfromtheargument.MASKhasthevalueINDEX BITSMANTISSA BITSShifttheresultingintegertotherighttoproduceanindexoftherequiredlength.Asanexamplecase,considerindexingatablewhenourpairwisedistancesrangefrom0.5to127.9Thiswouldbereasonableforashortrangepairpo-tentialwithacutoffofabout11.3.Theexponentsofthe¯oatingpointrepresentationsrangefrom1to6.Inbias127orbias1023notationtheseexponentsareuniqueintheirthreelowestorderbits.Notonlyaretheyunique,butforthisrangeofvaluestheycom-pleteallpossiblethreebitstrings.Hence,ifwecreateatableusingthreeasthevalueforEXP BITSandsay17forMANTISSA BITS,wecancreateatableofsize2withnearlyeveryentry®lled.Wewouldalsobeaccurateinthevalueoftoaprecisionofaboutintheworstcase.Thekeyhereistomakethemappingfunctionaseasytocomputeaspossible.Whileitappearsthattheextractfunctionismorecomplexthanthehashfunction,thefactisthattheintegercastisanexpensiveoperation.Wehavefoundthehashfunctionrequiresabout5timesasmuchtimeastheextractfunction.2.3.MultipletablesMDcodesoftendealwithmorethanonespeciesofatom,makingitnecessarytousemorethanonetable.Thiscanbedonequitereadilyforexamplebymakinganarrayoftables.The®rsttwodimensionsofthearraymightbeindexedbyacodeforthespeciesoftheatom,andthethirddimensionwouldbethetableitself.Wewouldarrangethepointersinthetoptwolevelsofthearraysothataninteractionbetweenatomswouldaccessthesametable,independentoftheorderinwhichthespeciesarespeci®ed.Anaccesstothetablemightbethefollowing:e=tablesstype1ype1type2ype2hashwheretype1andtype2refertotheatomicspecies.2.4.InterpolationFinally,onecouldinterpolatebetweensuccessivepointsinthetableformoreaccurateresultsfortheforcesandpotentials.Alinearinterpolationismostoftenallthatisrequired,althoughmorecomplicatedinterpolationshavebeenuseded13].Thesameeffect D.Wolff,W.G.Rudd/ComputerPhysicsCommunications120(1999)20 Fig.1.Representationofdoubleprecisionindexextraction.INDEX BITSarethebitsthatareactuallyextractedtoindexthetable.canoftenbeachievedatlittleextracostbyincreasingthesizeofthetables.2.5.ConstructingtablesConstructionofthetableisdonebeforethesimula-tionstarts.Wehavefoundthatforrealisticproblems,computingthetablestakesanegligibleamountoftimeandsowecomputethemwitheachrun.Itwouldmakesensetostoretablesthatareexpensivetocomputeondisk.Hereissamplecodetocomputeatable.maxrsqr:=index:=hashrsqrenergyTableableindex:=computeEnergyrsqrforceTableableindex:=computeForcersqrTheroutineabovebuildsatableformaxwheremaxisthecutoffdistanceforthepotential,andissomeminimumdistancewhereitisknownthattwoatomswillneverbelessthanapart.This,ofcourse,dependsonthedetailsofthesimulation.Whenusingthehashfunctionssuppliedabove,someofthistablewouldbeleftunusedbecausethedistancesthatmaptothelowerindicestotheta-blewillneverbeneeded.ThisisalsooftenalthoughnotnecessarilythecasewiththeDPex-tractionroutine.Dependingonthesizeofthetablethatisneededthewastedspacemaynotbealargeconcern.Onecouldmodifythehashfunctionsothatnospaceiswasted,butthatwouldbeatthecostofadditionalarithmeticwhenaccessingthetable.Itispossibletousehashfunctionsalongwithchain-ingmethodsforcollisionresolutiontohelptoelimi-natethewastedspace.Weexperimentedwiththisap-proachanddeterminedthattheextracodingandcom-putationaleffortrequiredprobablydoesnotjustifythe 2.03.04.05.06.07.08.0 -120.0-80.0-40.00.040.0120.0240.0280.0LJ Energy Fig.2.TheLennard-Jonestable-basedpotential.Thespacingofthetablevaluesismuchlargerthanwouldnormallybeusedinpractice.smallsavingsinmemory.3.ErrorsanderroraccumulationTheuseoftableswithoutinterpolationchangesthepotentialsandforcefunctionstosequencesofstepfunctionsseeFig.2.Again,atthecostofsomearithmetic,wecanimprovetheaccuracyoftabulatedpotentialsviainterpolation.Inourcodes,westoreboththepotentialandtheforcefunctionssothatweneednotdifferentiatethepotentialfunctionnumericallytoobtainforces.3.1.EstimationoferrorinthepotentialInordertodetermineamorequantitativepictureoftheerrorsintroducedbytheuseofatabulatedforce®eld,weconsiderageneralexample.Mostofthecom-monpotentialsinusetodayhavetermsthatdependei-theronaninversepowerofthedistanceorex- D.Wolff,W.G.Rudd/ComputerPhysicsCommunications120(1999)20ponentiallyondistance.Wewillconsiderbothcases.Firstweexpresstherelativeerrorforaninversepowerterminasingletablelook-upastable (r)=1�1p 1=1�rp whereistheactualpositionoftheparticleandthetruncatedorroundedvalueusedforatablelookup.Wecanaccessthetablebyeitherroundingthevaluetothenearesttableentry,orbysimplytruncatingroundingdown.Inthecaseoftheroundingmethodofaccessingthetable,becauseoftheincreasing,neg-ativeslopeof1,themaximumpossibleerroroc-curswhen2whereisthetablespacingseeFig.3.Hence Sinceispositiveandisalwaysgreaterthan2,wecanwritetheupperboundas Thisgivesanupperboundfortheerrorin1typepotential.Themaximumvaluethatcantakeonde-pendsonthesmallestpossiblevalueforforagiven.Todeterminethiswemustlooktotheparticularsubstancebeingsimulated.Thepairdistributionfunc-tionforasubstancegivesusinformationontheenvironmentofagivenatom.Sincemostatomsactlikehardspheresforsmallvaluesof,thepairdistri-butionfunctionmustgotozeroatsomesmallvalue,letuscallit.Hence,wecanexpectthatnoatomwilleverbewithinofanother.ThisgivesusalowerboundoninEq..Hencewecanwritethemaximumpossibleerrorinthepotentialasmax Forourargoncasestudy,thepairdistributionfunctiongoestozeroatabout3.0.Soforaof0.001=6,max001.Thisisrespectablewhencomparedwiththeaccuracytowhichmanypotentialparametersareknown.Also,itshouldbenotedherethatthismayalsobeontheorderoftheaccuracygivenbytheintegrationmethodused.Mostintegrationmethodsconservetotalenergytowithintheorderofthetimestepsize.Andperhapsmoreimportantly,mosttablelookupswillhavesubstantiallysmallererror.Inaddition,fortheroundingmethod,wecanerroneithersideoftheactualvalueforthepotentialFig.3.Thatis,thetabulatedvaluemightbelargerorsmallerthantheactualvalue.Hence,wecanexpectthatsomeoftheoverallerrorwillbecanceledoutbythiseffect.Asimilaranalysiscanbedoneforanexponentialpotential,whichyields3.2.AccumulationoferrorsForapairpotential,thetotalpotentialenergyforasingletimestepissimplyasumofpairwisepotentialvalues,whereisthedistancebetweenatomandatom.Forthecaseoftheexponentialpotential,therelativeparticleparticleerrorisindependentofthepairwisedistance.givesanexpressionforthe1potentialthatisalsoindependentoftheparticleparticledistance.Now,sincethemaximumrelativeerrorofeachterminthesumisindependentofthepairwisedistance,therelativeerrorinthesumisthesameastherelativeerrorofeachterm.Hence,Eqs.alsoholdforthetotalpotential,Next,weconsiderthetotalforceonagivenatom, D.Wolff,W.G.Rudd/ComputerPhysicsCommunications120(1999)20 Fig.3.Representationofasingletableentry.Notethattheerrorislargestwhenfortruncation,andwhen2fortheroundingmethod.where ItcanbeshownthattherelativeerrorinforthepotentialtypeisthesameasEq.exceptwithanexponentof1.Similarly,itiseasytoshowthattherelativeerrorinfortheexponentialtypepotentialisthesameasEq.Now,theforceisavectorquantity,sothesumma-tionisalittledifferenthereasopposedtothetotalpotentialenergy.However,onecanstillputanupperboundonthetotalerrorbyassumingthatalloftheer-rorspointinthesamedirection.Againthesameargu-mentappliesasabove.Therelativeerrorintheforceforasingleparticleis forthe1potential.Fortheexponentialpotentialasimilaranalysisrevealsthesameupperboundas3.3.ErrorsinmeasuredvaluesNowthatwehaveanupperboundontherelativeerrorfortheforceandpotentialforagiventimestep,letusconsiderhowthoseerrorspropagateintosome 246810 0e+002e-044e-046e-048e-04Energy (kJ/mol * 100) Actual error Product 246810 12g(r) Fig.4.Thescaled,actualerrorandpairdistributionfunctiontakenfromanMDsimulationofLJliquidargon,whereisthedepthofthepotentialwell.Thesolidlineistheproductoftheerrorandthepairdistributionfunction.ofthequantitiesthataremeasured.Thepressure,tem-perature,stress,andtotalenergyalldependonvarioussummationsinvolvingtheforce,positions,andveloc-itiesoftheparticles.Wehavealreadydiscoveredhowtheerrorsinindividualtablelookupspropagatetothetotalforceandpotentialenergy.Thetotalenergyisgivenby 2Xjmjv2j:( D.Wolff,W.G.Rudd/ComputerPhysicsCommunications120(1999)20 100200300400500 -1.2-1.1-1-0.9-0.8 Fig.5.Totalenergyperatomversustimeforseveralsizesofpotentialtables.Wehavealreadyputanupperboundontheerrorinthe®rstterm,itremainstoputanupperboundontheerrorinthesecond.Thevelocityisobtainedfromthetotalaccelerationbyintegratingoversomeshorttimeinterval.Therearemanymethodsofdoingthiss1],butforsimplicityandinordertoconsiderupperbounds,wewillconsiderjustaTaylorexpansionofthevelocity,Foragiventimestepthereissomeerrorintheforceduetotheuseofatabulatedforce®eld.Thisisequivalenttoaddingatermtotheacceleration,whereistherelativeerrorboundintheaccelerationwhichisequaltotherelativeerrorboundintheforce.Now,thistermisofordert"whichissimilarorderastheerrorinthepotentialenergyterm.Intheexpres-sionforthetotalenergy,thistermissquared,whichmakestheerrorofordert".Regardless,wecanconsidertheerrorinthetotalenergytobeproportional.Weuseinsteadofbecausethebinomialisraisedtoahigherpowerin.Inanycase,aswillbeseeninthenextsection,wearemostinterestedinhowtheerrorchangeswithrespecttothetablesize,andforthispurpose,theaboveargumentissuf®cient.3.4.ErrorwithrespecttothetablesizeWeexpandthebinomialinEq. 2r�p(p�1) 2Dr wherewehavedroppedtheprimeon.Sinceforallreasonable,weknowthat Now,wewouldliketotruncatetheseriesto®rstorder,butwemust®rstconvinceourselvesthattheupper-boundstillholdsinthatcase.Inotherwords,wemustshowpx;for01.Since1,thisistruefor=1and=0.Considerthe®rstderivativeofbothsideswithrespecttoSinceweareonlyconcernedwiththerange01andsince1,thisisclearlytrue.Havingestab-lishedthattheslopeoftheleft-handsideisalwayssmallerorequaltotheslopeoftheright-handside,andthattheendpointsofourrangeofinterestsatisfy,wecansaythatEq.holdsforall0TruncatingEq.to®rstorderin,wehave 2rpk k=Nwheremaxistherangeofvaluesofinthetable,andisthenumberoftableentriesorthesizeofthetable.Thisindicatesalinearrelationshipbetweenthein-versetablesizeandtherelativeerror. Byreasonable,wemeananyvalueofwhichwillnormallybeencounteredinasimulation.Thisisboundedatthelowendbythepointatwhichthepairdistributionfunctionbecomesnegligible. D.Wolff,W.G.Rudd/ComputerPhysicsCommunications120(1999)203.5.ErrorsinpracticeInpractice,theupper-boundsthatwehavederivedintheprevioussectionsareveryhigh.Infact,theactualerrorsthatareincurredoverallwillbeofamuchlowermagnitude.Thisisbecauseofseveralreasons.First,theerrorinatablelookupcanbeoneithersideoftheactualvalue.Hence,wecanexpectsomeªwashingoutºoftheerrorswhenmultiplevaluesfromthetablearesummed.Second,atleastfor1potentials,theerrorisre-ducedforatomsthatarefartherapartthan.Aswasstatedpreviously,istheminimumdistancethattwoatomsmightapproacheachother.Thisisindicatedbythepointatwhichthepairdistributionfunctionforthesubstancebeingsimulatedgoestozero.Fig.4showsagraphofthepairdistributionfunc-tionforliquidargonoverlaidwiththemaximumer-rorasafunctionof.Theregionoflargesterrorintheforceiswherethepotentialissteepest.Alsoin-cludedinthegraphistheproductofthetwofunc-tions.Thepairdistributionfunctionispropor-tionaltotheprobabilityof®ndinganatomatadis-tancearoundanygivenatom.Theproductofthetwofunctionsshowsamaximumatapproximatelythelo-cationwherethetwocurvescross,andgoestozero.Thisdemonstratesthatdistancesbetweenneighboringatomsareseldomintheregionoflargesterror.Finally,fortheforce,itisveryunlikelythatalloftheforcesonanatomwillpointinthesamedirec-tion.Hence,weshouldexpecttheabsoluteerrorstosuminaªrandomwalkºfashion.Forpurposesofanalysiswehavemadethesimpli®-cationthatthepotentialenergyismadeupofasingletermthatisalwayspositive.Inreality,thepotentialenergywillbeasumofpositiveandnegativecon-tributionsfromtheseterms.Insuchcaseswewouldneedtobemuchmorecarefulinusingtherelativeer-rorsincethepotentialenergyandtheforcemaybezero.Nevertheless,wecanstilldoasimilaranalysisoneachtermofthepotentialandthencombinetheerrorsafterthesummation. 0.0000.0050.0100.015 0.000.100.200.30Relative Error Fig.6.RelativeerrorintotalenergyversusinversetablesizefortheexperimentinFig.5.4.Experimentsonsamplesystems4.1.VariationoftablesizeFig.5showstheresultsoftestsinwhichwevar-iedthesizesofthepotentialandforcetables.Eachcurveisa50000-steprunforLJliquidargon.Thesys-temsconsistedof10648atoms.WeusedRapaport'ss2]andaleapfrogintegratorwithreducedtimestepsof0.001.Thetablesizesaredeterminedbytheparameter usedintheextractmappingdescribedabove.Weseethatresultsusingtablesofsizes128,256,and512elementsforeachtablemovetowardtheaggregateofresultsusinglargertables,labeledªOth-ersºinthegraph,whichareessentiallyindistinguish-able.Thelargertablescontained2elements,for10,12,14,16,18,and20.Wecomputedtheaverageenergyforthelast25000timestepsforeachtablesizeinFig.5.Takingtheav-eragevalueofthethreecurves=16,18,and20whichareindistinguishableastheªexactºenergy,weplottedrelativeerrorversus1whereisthetablesizeinFig.6.Thisdemonstratesthelinearrelation-shippredictedinEq..Hence,wecanseethattheexpressionforsingleparticleerrorcontributesinasimilarwaytooverallerrorsintheenergy.Theseresultsindicatethattablesofsize2orlargeraremorethanadequateatleastforthetotalenergyincreasingthetablesizebeyondthathasnonoticeable D.Wolff,W.G.Rudd/ComputerPhysicsCommunications120(1999)20effectontheresultsforthetotalenergy.Thisalsoindicatesthatinterpolationmaynotbenecessary.Weobtainedsimilarresultsinsimulationsof216siliconatomsintheliquidstateusingatableforthepairpartoftheStillingerWeberpotentiall14].Onewouldexpectthatsolongasthetablesarelargeenoughtoprovidethesameresultsascomputedpotentials,itmakesnodifferencehowlargethetablesare.However,todosowouldbetoignorethecacheeffect.IntheLJargoncase,forexample,wefoundthatthetimepercalltotheroutinethatcomputesalltheforcesinatimestepincreasedby50%whenwasincreasedfrom14to16,whichrepresenttablesizesof16Kand64Kentries,respectively,onanultrasparcprocessorwitha1Mbytecache.Storagerequirementsforthesetableswere256Kand1Mbytes,respectively.Therewasnofurtherincreaseinexecutiontimeforlargertablesupto2entries.4.2.ComparisonofoutputWenowdescribeourcomparisonsoftheresultsfromsimulationsthatusetabulatedpotentialsandsim-ulationsinwhichthepotentialsandforcesarecom-putedwitheverystep.Table1showstheresultsforLJliquidargonusingtheMoldyy10]code.Wenotethat,despitetherelativelyshorttimesusedforequilibrationandforaccumulatingaverages,thevariationswithinenergyarewellwithinacceptablelimitsandthestan-darddeviationsinenergyandtemperatureshowtheexpectedbehaviorsasthesizeofthesystemincreases.Fig.7showsthesimulationtimedependenceofthetotalenergyperparticleforLJliquidargonusingthetwomethods.Thesimulationwasfor50000stepsfor10648atomsusingtheRapaportcode.Whilethe¯uctuationsinenergybetweenthetwosimulationsarenotsynchronizedafterabout10000steps,thereisnosigni®cantdifferencebetweenthem.Fig.8showsthatthereisnosigni®cantdifferenceinstructurebetweentheresultsproducedbythetwomethods.Theinsetshowsasmalldifferenceintheheightsofthesecondpeak.Thisisclosetobeingwithinthenoiseofthefunction,butcouldbesigni®cantfordetailedstructuralstudies.Weconcludethat,atleastfortheequationofstateandstructurestudiesinwhichgeneralstructuralinformationisneeded,thereisnoimportantdifferenceinresultsbetweenthetwomethods.5.PerformanceWeturnnowtoananalysisoftherelativeperfor-manceofthetwomethods.Wehavefocusedonthreeformsforpotentials:polynomialsin1,suchastheLJpotential,eralizedºpotentials,whicharecombinationsofpoly-nomialsandexponentials,andthree-bodypoten-tials,suchastheStillingerWeberpotentialforsili-con.Thelatterclassinvolvesconsiderablelogicandarithmetictocomputethethree-bodycontributiontothepotentialinadditiontothepolynomialsandexpo-nentialsinthetwo-bodycontributioninthepotential.Table2showstheperformanceresultsfortheLJpotential.Thesimulationswereonasystemofar-gonatoms.Kerneltimesrefertotheamountoftimespentinthesubroutinethatcalculatesthepotentialandforce.Therelativeperformanceisreproducibleacrosscodesandmachines.WeusedtheRapaportcodeandand10]onHP,Sun,andDECworkstations,withlittlevariationindifferencesinrelativeperformance.Theseresultsmightseemdisappointingat®rst,sinceitwouldseemthattheLJpotentialwouldre-quireamoderateamountofarithmetic.However,awell-codedLJcalculatorrequiresonlythreemultipli-cationsandasubtraction,ifisalreadyavailable,andagoodcompilercanrearrangetheloopthroughtheatomsandmanagethecachetomakesurethatthepipelineisalwaysfull.Furthermore,onmodernpro-cessors¯oating-pointarithmeticisoftenfasterthantheintegerarithmeticandrandomlogicourextractroutineexecutes.Hence,onewouldexpecttoseeverylittleornoimprovementwithtabulatedpotentialsforaLJpotential.Table3showstheresultsforageneralizedpoten-tialwithparametersforliquidargon.Thegeneralizedpotentialhasthefollowingform: rD r4�E r6�F ThisexperimentwasdonebyinsertingatableintoapreviouslyestablishedMDcode.ThecodethatweusedwascalledªMoldyºdyº10].Thetimingsforthekernelroutineshowthatthetablemethodisaboutafactoroffourfasterinkerneltimethandirectcomputation.Ourtestsononeofthemachinesonwhichthesecalculationsweredone,a170MHzul-trasparc,showedthat¯oating-pointmultiplications, D.Wolff,W.G.Rudd/ComputerPhysicsCommunications120(1999)20 15253545 -1.15-1.1-1.05-1 Computed Fig.7.TotalenergyperatomusingcomputedandtabulatedpotentialsandforcesforaLJargonsystem.Table1Resultsfortableandnon-tablebasedrunsofliquidargondensity0andinitialtemperature120KwiththeLennard-Jonespotential.Averagesareover45000timestepsafter5000timestepsofequilibrationwithatimestepsizeof0.01ps.Noscalingofthetemperatureisperformed NumberofatomsEnergymolStd.Dev.TemperatureStd.Dev.Type 152.640.27140.535.53Computed152.750.28140.675.52Tabulated275.850.42150.384.06Computed275.970.41150.524.06Tabulated785.900.74140.282.57Computed785.730.72140.092.48Tabulated1715.211.11132.871.73Computed1715.051.08132.851.78Tabulated8447.272.43134.130.80Computed8448.922.55134.220.82Tabulated15930.543.38139.340.54Computed15931.173.17139.330.57Tabulated Table2KerneltimeandtimepercalltotheforceroutinefortheLJliquidargonexperiments ExperimentTimeinTimeperKernelforcecall Table,2elements87714.9Table,2elements94616.3Table,2elements161929.5Computed99817.3 squareroots,andtranscendentalfunctionevaluationsrequire0.006,0.21,and1.8s,respectively.Theex-ponentialfunctionsmakethedifference.Thespeedupisconsiderablylesswhenoneconsid-erstotalexecutiontime.Thisisbecausethecompu-tationofthepotentialisonlyonepartofthecom-putationsthatneedtobedoneinagiventimestep.Integrationoftheforce,measurementofdifferentpa-rameters,movementofatomsbetweencells,andoc-casionalneighborlistcalculationmakeuptherestofasingletimestep.Forexample,considerthespeedupofagiventimestepwhenthespeedupinthekernelis4,Speedup=kernel 1 kernelIftheothercalculationsinatimestepareequallyex-pensiveasthekernelcomputations,thenwehaveSpeedup=kernelkernel kernel kernel So,overallwewouldseeaspeedupofonly1.6asopposedto4.Table4showstimesforsimulationsoftwospeciesofatoms,whichrequires3tables.Hereweusedapotentialofthe6-expform, D.Wolff,W.G.Rudd/ComputerPhysicsCommunications120(1999)20Table3Timingsfortableandnon-tablebasedrunsofliquidargonusingageneralizedpotential.Thekernelreferstothesubroutinethatcalculatesthepotentialandforceforalistofpairwisedistances.Alltimesareinseconds NumberofAtomsTimeinkernelTotaltimeTimepertimestepType 100002969.2824759.334.95Tabulated1000010944.7633467.906.69Computed50002204.5317783.343.56Tabulated50007775.0923188.804.64Computed1000370.913310.070.66Tabulated10001617.314436.660.89Computed500203.491463.600.29Tabulated500689.001869.680.37Computed (r)=�A whichissimilartothegeneralizedpotential.Whenonecomparesthetimingsforthispotentialtothetim-ingsfortheLennard-Jonespotentialthesigni®canceoftheexponentialtermwithrespecttospeedupbe-comesclear.Weobservethattheresultsat®rstglancedonotindicateaspeedupaslargeasexpected,sincewithoneortwomoremultiplicationswewouldhavethegeneralizedpotentialforwhichwemeasuredafac-torof4speedup.Thisisbecausethereisnowmorethanonetable,andextralogicisincludedinthetablelookupstagetodeterminewhichtabletouse.Therea-sonthisisnecessaryisduetospeci®csoftheMoldycode.Hence,someoptimizationsbythecompilerfortheinnermostlooparenotimplemented,nottomen-tiontheadditionalinstructionsneededtoexecutecon-ditionalstatements.Nevertheless,thespeedupdemon-stratesthesigni®cantcontributionoftheexponential.Table5showstheperformanceofthemethodsonsilicon,4096atomsfor50000timestepsusingtheRa-paportcode.Ingeneratingthecodefortheforcecal-culations,wewereabletotabulateonlythetwo-bodyterms.Thethree-bodyterms,asusuallywritten,re-quirethecalculationofthecosineoftheanglebetweeneachtripletofatomsthatarewithinacutoffdistanceofeachother.Thecosineisadotproductbetweentwoofthevectorslinkingtheatoms.Thereislittlepointinconstructingtablesforthesebecausetheywouldre-quireatableindexedby,andcos.Hence,thetablewouldhavetobethree-dimensionalwhichwouldincreasethestoragerequirementsbyaboutapowerofthree.Also,verylittlearithmeticissavedsincemuchoftheworkisindoingthedotproduct. 051015 0g(r) Fig.8.Pair-correlationfunctionfromcomputedandtabulatedpotentialsandforcesforaliquidsiliconsystem.Nakanoetal..15]havederivedameansforsepa-ratingthethree-bodyforceintosetsoftwo-bodyin-teractions.Itwouldbeinterestingtoseethegainsthatcanbemadethroughusingtheirtechniquewithatab-ulatedpotential.6.ConclusionsInTable6weprovideapproximateratiosoftimesinforcecalculationsforcomputedandtabulatedpo-tentials.Onemessageisclear.Oneshouldcertainlyusetablesforpotentialswhenthepotentialinvolvestranscendentalfunctioncallsandotherintensivecom-putationswithrelativelylittlelogic.Thecaseformorecomplexpotentials,suchasthree-bodypotentialsandproblemsthatrequiremultipoleexpansionsisstill D.Wolff,W.G.Rudd/ComputerPhysicsCommunications120(1999)20Table4Timingsfortableandnon-tablebasedrunsofaSiOsystemwith192Siliconatomsand384Oxygenatoms.Weseeanimprovementinspeedofaboutafactorof1.18inexecutiontimeand1.72inkerneltime.EvaluationoftheCoulombinteractionwasomittedforthepurposeofthisstudy. TimeinTotalSecondsperTypekerneltimetimestep 370.801529.640.77Tabulated639.151810.450.91Computed Table5Kerneltimeandtimepercalltotheforceroutineforthesiliconexperiments ExperimentMillisecondsperforcecallperatom Table,216atoms,5000steps0.79Table,4096atoms,50000steps0.75Computed,216atoms,5000steps0.98 Table6Approximateratiosforcomputingpotentialsversustablelookupsforthreeclassesofpotential PotentialComputeTabletime Lennard-Jones16-exp7Generalized4StillingerWeber5 cloudy.WehaveshownthattheuseoftabulatedpotentialscanbeeffectiveinspeedingupMDsimulations.Wehaveoutlineddifferenttradeoffsthatneedtobetakenintoaccountbeforeimplementingatableintoapar-ticularcode.Thenumberofdifferentspeciestobesimulated,complexityofthepotential,andcachesizeareallimportantconsiderations.Wehavealsoinves-tigatedtheerrorsintroducedbytabulatedpotentials.Surprisinglysmalltablesofafewthousandelementsproduceresultsthatareindistinguishablefromthoseobtainedusingdoubleprecisionarithmetic.Referencesnces1]M.P.Allen,J.P.Tildesley,ComputerSimulationofLiquidsOxfordSciencePublications,Oxford,198719872]D.C.Rapaport,TheArtofMolecularDynamicsSimulationCambridgeUniv.Press,Cambridge,199519953]R.W.Hockney,J.W.Eastwood,ComputerSimulationUsingParticlesAdamHilger,Bristol,198819884]Privatecommunicationswithseveralinvestigators..5]P.Clapp,J.Rifkin.Personalcommunication;seehttp://www.ims.uconn.edu/centers/simul/l/6]D.Turner,J.Morris,seehttp://cmp.ameslab.gov/ Theory/cmd/alcmd 7]D.E.SandersM.S.Stave,L.S.Perkins,A.E.DePristo,Comput.Phys.Commun.70708]M.S.StaveD.E.Sanders,T.J.Raeker,A.E.DePristo,J.Chem.Phys.9339]A.Y.Toukmaji,J.JohnA.Board,Comput.Phys.Commun.un.10]K.Refson,MoldyUsersManual,seehttp://www.earth.ox.ac.uk/~keith/moldy.htmltml11]M.Nelson,W.Humphrey,A.Gursoy,A.Dalke,L.Kale,R.D.Skeel,K.Schulten,J.Supercomput.Appl.&HighPerformanceComputingg12]D.M.Beazley,P.S.Lomdahl,Parall.Comput.20013]T.A.Andrea,W.C.Swope,H.C.Andersen,J.Chem.Phys..14]F.H.Stillinger,T.A.Weber,Phys.Rev.B313115]A.Nakano,P.Vashishta,R.K.Kalia,Comput.Phys.Commun.77