Download Presentation
## Javier Junquera

- - - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - - -

**Advanced Computing**Starting-up: Howtosetuptheinitialconditionsfor a Molecular DynamicSimulation Javier Junquera**The initial configuration**In molecular dynamicsimulationsitisnecessarytodesign a startingconfigurationforthefirstsimulation • Initial molecular positions and orientations • Initialvelocities and angular velocities Forthefirstrun, itisimportanttochoose a configurationthat can relax quicklytothestructure and velocitydistributionappropriatetothe fluid Thisperiod of equilibrationmust be monitoredcarefully, sincethedisapperance of theinitialstructuremight be quite slow**The initial configuration:**more usual approach, start from a lattice Almostanylatticeissuitable. Historically, theface-centeredcubicstructure has beenthestartingconfiguration Thelatticespacingischosen so theappropriateliquidstatedensityisobtained Duringthecourse of thesimulation, thelatticestructurewilldisappear, to be replacedby a typicalliquidstructure Thisprocess of “melting” can be enhancedbygivingeachmolecule a smallrandomdisplacementfromitsinitiallatticepoint**The initial configuration:**more usual approach, start from a lattice Almostanylatticeissuitable. Historically, theface-centeredcubicstructure has beenthestartingconfiguration A supercellisconstructedrepeatingtheconventionalcubicunitcell of the FCC lattice times alongeachdirection Thenumber of atoms in thesimulation box, , isaninteger of theform , whereisthenumber of FCC unitcells in eachdirection**Units for the density**Forsystemsconsisting of justonetype of atomormolecule, itis sensible to use themass of themolecule as a fundamental unit • Withthisconvention: • Particlemomenta and velocitiesbecomenumericallyidentical • Forces and accelerationsbecomenumericallyidentical In systemsinteractingvia a Lennard-Jone potential, thedensityisoftenlyquoted in reducedunits Assumingan FCC lattice (4 atoms in theconventionalcubic), and , thenwe can compute thelatticeconstantfromthereduceddensity And then, thelatticeconstantwill come in thesameunits as theonesusedto determine**Units for the length of the supercell**Ifthelatticeconstant of theconventionalcubicunitcell of an FCC latticeis , thenthelength of theside of thesupercellis However, theimplementation of theperiodicboundaryconditions and thecalculation of minimumimagedistancesissimplifiedbythe use of reducedunits: thelength of the box istakento define the fundamental unit of length of thesimulation, In particular, theatomiccoordinates can be definedusingthisunit of length, so nominallytheywill be in therange**Implementation of a fcc lattice**Thesimulation box is a unit cube centred at theorigin Thenumber of atoms in thesimulation box, isaninteger of theform , whereisthenumber of FCC unitcells in eachdirection**Initial velocities**For a molecular dynamicsimulation, theinitialvelocities of allthemoleculesmust be specified Itis usual tochooserandomvelocities, with magnitudes conformingtotherequiredtemperature, corrected so thatthethereis no overallmomentum Thedistribution of molecular speedsisgivenby For a derivation of thisexpression, readFeynmanlecturesonPhysics, Volume 1, Chapter 40-4 Probabilitydensityforvelocitycomponent Similar equationsapplyforthey and zvelocities**Normal distributions**The normal distributionwith mean and varianceisdefined as A randomnumbergeneratedfromthisdistributionisrelatedto a numbergeneratedfromthe normal distributionwithzero mean and unitvarianceby The Maxwell-Boltzmanndistributionis a normal distributionwith Ifwetakethemass of theatomsormolecules as**Units of temperature and velocity**Thetemperatureisusuallygiven in reducedunits (rememberthatisusuallytabulated as , in K) In thissystem of units, thevelocities are given in units of squareroot of thetemperature**Implementation of random numbers following a Gaussian**distribution**Implementation of random numbers uniformly distributed**between 0 and 1**Typical system sizes**Thesize of thesystemislimitedbythe: - availablestorageonthe host computer - speed of execution of theprogram Thedoubleloopusedtoevaluatetheforcesorthepotentialenergyisproportionalto Specialtechniquesmay reduce thisdependencyto , buttheforce /energyloopalmostinevitablydictatestheoverallspeed.**How can we simulate an infinite periodic solid? Periodic**(Born-von Karman) boundary conditions We should expect that the bulk properties to be unaffected by the presence of its surface. A natural choice to emphasize the inconsequence of the surface by disposing of it altogether Supercell + Born-von Karman boundary conditions**Periodic (Born-von Karman) boundary conditions**The box isreplicatedthroughoutspacetoformaninfinitelattice In thecourse of thesimulation, as a moleculemoves in the original box, itsperiodicimage in each of theneighboring boxes moves in exactlythesameway Thus, as a moleculeleavesthe central box, one of itsimageswillenterthroughtheoppositeface There are no walls at theboundary of the central box, and no surfacemolecules. This box simplyforms a convenient axis systemformeasuringthecoordinates of themolecules Thedensity in the central box isconserved Itisnotnecessarytostorethecoordinates of alltheimages in a simulation, justthemolecules in the central box**Periodic (Born-von Karman) boundary conditions: influence on**the properties Itisimportanttoaskiftheproperties of a small, infinitelyperiodic, systemand themacroscopicsystemwhichitrepresents are thesame For a fluid of Lennard-Jones atoms, itshould be possibletoperform a simulation in a cubic box of side without a particlebeingableto “sense” thesymmetry of theperiodiclattice**Periodic (Born-von Karman) boundary conditions: drawbacks**and benefits Benefits Drawbacks • Inhibitstheoccurrence of long-wave lengthfluctuations. For a cube of side , theperiodicitywillsuppressanydensity wave with a wave lengthgreaterthan . • Be careful in thesimulation of phasetransitionswheretherange of criticalfluctuationsismacroscopic, and of phonons in solids. • Commonexperience in simulationworkisthatperiodicboundaryconditionshavelittleeffectontheequilibrium of thermodynamicproperties and structures of fluids: • - awayfromphasetransitions • - wheretheinteractions are short-range Iftheresources are available, itisalways sensible toincreasethenumber of molecules (and the box size, so as tomaintainconstantdensity) and rerunthesimulations**Periodic (Born-von Karman) boundary conditionsand external**potentials Up tonow, wehaveassumedthatthereis no externalpotential (i.e. no term in theexpansion of thepotentialenergy) Ifsuch a potentialispresent: - itmusthavethesameperiodicity as thesimulation box or - theperiodicboundaryconditionsmust be abandoned**Periodic (Born-von Karman) boundary conditionsin 2D and 1D**In some cases, itisnotappropriatetoemployperiodicboundaryconditions in each of thethreecoordinatedirections In thestudy of surfaces(2D) In thestudy of wiresortubes (1D) Thesystemisperiodiconly in the planes paralleltothesurfacelayer Thesystemisperiodiconlyalongthedirection of thewireortube**Calculating properties of systems subject to periodic**boundary conditions Heartof a Molecular Dynamicor Monte Carlo program: - calculation of thepotentialenergy of a particular configuration - in Molecular Dynamics, compute theforcesactingonallmolecules Howwouldwe compute theseformolecule 1 Wemustincludeinteractionsbetweenmolecule 1 and everyothermolecule (orperiodicimage) Thisisaninfinitenumber of terms!! Impossibletocalculate in practice!! For a short-rangepotentialenergyfunction, wemustrestrictthissummationbymakinganapproximation**Cutting the interactions beyond a given radius to compute**the potential energy and forces Thelargestcontributiontothepotential and forces comes fromneighboursclosetothemolecule of interest, and for short-rangeforceswenormallyapply a sphericalcutoff Thatmeans: Molecules 2 and 4Econtributetotheforceon 1, sincetheir centers lie insidethecutoff Molecules 3E and 5F do notcontribute In a cubicsimulation box of side , thenumber of neighboursexplicitlyconsideredisreducedby a factor of approximately Theintroduction of a sphericalcutoffshould be a smallperturbation, and thecutoffdistanceshould be sufficientlylargetoensurethis. Typicaldistance in a Lennard-Jones system:**Difficulties in defining a consistent POTENTIAL in MD**method with the truncation of the interatomic pot Thefunctionused in a simulationcontains a discontinuity at : Whenever a pair of moleculescrossesthisboundary, the total energywillnot be conserved We can avoidthisbyshiftingthepotentialfunctionbyanamount Thesmalladditiontermisconstantforanypairinteraction, and doesnotaffecttheforces However, itscontributiontothe total energyvariesfrom time stepto time step, sincethe total number of pairswithincutoffrangevaries**Difficulties in defining a consistent FORCE in the MD method**with the truncation of the interatomic potent Theforcebetween a pair of moleculesisstilldiscontinuous at For a Lennard-Jones case, theforceisgivenby And themagnitude of thediscontinuityisfor It can cause instability in thenumericalsolution of thedifferentialequations. Toavoidthisdifficulty, a “shiftedforcepotential” has beenintroduced Thediscontinuitynowappears in thegradient of theforce, not in theforceitself**Difficulties in defining a consistent FORCE in the MD method**with the truncation of the interatomic potent Thedifferencebetweentheshifted-forcepotential and the original onemeansthatthesimulation no longercorrespondstothedesiredmodelliquid However, thethermodynamicpropertieswiththeunshiftedpotential can be recoveredusing a simple perturbationscheme**Computer code for periodic boundaries**Letusassumethat, initially, themolecules in thesimulation lie within a cubic box of side , withtheorigin at its centre. As thesimulationproceeds, thesemoleculesmoveabouttheinfiniteperiodicsystem. When a moleculeleavesthe box bycrossingone of theboundaries, itis usual toswitchtheattentiontotheimagemoleculeenteringthe box, simplyaddingorsubtractingfromtheappropriatecoordinate**Loops in Molecular Dynamic (MD) and Monte Carlo (MC)**programs Looponallthemolecules(from ) For a givenmolecule , loopoverallmoleculestocalculatetheminimumimageseparation Ifmolecules are separatedbydistancessmallerthanthepotentialcutoff Ifmolecules are separatedbydistancesgreaterthanthepotentialcutoff Compute potentialenergy and forces Skiptotheend of theinnerloop, avoidingexpensivecalculations Time requiredto examine allpairseparationsisproportionalto**Neighbour lists: improving the speed of a program**Innerloops of the MD and MC programsscaleproportionalto Verlet: maintain a list of theneighbours of a particular molecule, whichisupdated at intervals Betweenupdates, theprogramdoesnotcheckthroughallthemolecules, butjustthoseappearingonthelist • Thenumber of paisseparationsexplicitlyconsideredisreduced. • Thissaves time in: • Looping through • Minimumimaging, • Calculating • Checkingagainstthecutoff**The Verletneighbour list**Thepotentialcutoffsphere, of radius , around a particular moleculeissurroundedby a “skin”, togive a largersphere of radius At firststep, a listisconstructed of alltheneighbours of eachmolecule, forwhichthepairseparationiswithin Theseneighbours are stored in a largeone-dimensional array, LIST Thedimension of LIST isroughly • At thesame time, a secondindexingarray of size , POINT, isconstructed: • POINT (I) pointstothe position in thearray LIST wherethefirstneighbour of molecule I can be found. • Since POINT(I+1) pointstothefirstneighbour of molecule I+1, then POINT(I+1)-1 pointstothelastneighbour of molecule I. • Thus, using POINT, we can readilyidentifythepart of thelast LIST arraywhichcontainsneighbours of I.**The Verletneighbour list**Overthenextfewsteps, thelistisused in force/energyevaluationroutine Foreachmolecule I, theprogramidentifiestheneighbours J, byrunningover LIST frompointto POINT(I) to POINT(I+1)-1 Itisessentialtocheckthat POINT(I+1) isactuallygreaterthen POINT(I). Ifitisnotthe case, thenmolecule I has no neighbours From time to time, theneighbourlistisreconstructed and thecycleisrepeated.**The Verletneighbour list**Overthenextfewsteps, thelistisused in force/energyevaluationroutine Foreachmolecule I, theprogramidentifiestheneighbours J, byrunningover LIST frompointto POINT(I) to POINT(I+1)-1 Itisessentialtocheckthat POINT(I+1) isactuallygreaterthen POINT(I). Ifitisnotthe case, thenmolecule I has no neighbours From time to time, theneighbourlistisreconstructed and thecycleisrepeated.**The Verletneighbour list**Thealgorithmissuccessfuliftheskinaroundischosento be thickenough so thatbetweenreconstructions: A molecule, such as 7, whichisnotonthelist of molecule 1, cannotpenetratethroughtheskinintotheimportantsphere Molecules, such as 3 and 4 can move in and out of thissphere, butsincethey are onthelist of molecule 1, they are alwaysconsidereduntilthelistisnextupdated.**Parameters of the Verletneighbour list: the interval between**updates Oftenfixed at thebeginning of theprogram Intervals of 10-20 steps are quite common Animportantrefinement: allowtheprogramtoupdatethelistautomatically: - Whenthelistisconstructed, a vector foreachmoleculeis set tozero - At subsequentsteps, the vector isincrementedwiththedisplacement of eachmolecule. - Thus, itstores, thedisplacement of eachmoleculesincethelastupdate - Whenthe sum of the magnitudes of thetwolargestdisplacementsexceeds , theneighbourlistshould be updatedagain**Parameters of the Verletneighbour list: the list sphere**radius Is a parameterthatwe are free tochoose As isincreased, thefrequency of updates of theneighbourlistwilldecrease However, with a largelist, theefficiency of the non-updatedstepswilldecrease Thelargerthesystem, the more dramatictheimprovement**Algorithm to check whether the update of the list is**required**Algorithm of the inner loop to compute forces and potential**energy