Javier Junquera

1. Advanced Computing Starting-up: Howtosetuptheinitialconditionsfor a Molecular DynamicSimulation Javier Junquera

2. 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

3. 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

4. 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

5. 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

6. 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

7. 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

8. 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

9. 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

10. 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

11. Implementation of the initial velocities

12. Implementation of random numbers following a Gaussian distribution

13. Implementation of random numbers uniformly distributed between 0 and 1

14. Typical system sizes Thesize of thesystemislimitedbythe: - availablestorageonthe host computer - speed of execution of theprogram Thedoubleloopusedtoevaluatetheforcesorthepotentialenergyisproportionalto Specialtechniquesmay reduce thisdependencyto , buttheforce /energyloopalmostinevitablydictatestheoverallspeed.

15. Boundary conditions

16. 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

17. 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

18. 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

19. 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

20. 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

21. 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

22. 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

23. 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:

24. 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

25. 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

26. 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

27. 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

28. Computer code for periodic boundaries

29. 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

30. 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

31. 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.

32. 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.

33. 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.

34. 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.

35. 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

36. 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

37. Algorithm to compute the Verletneighbour list

38. Algorithm to check whether the update of the list is required

39. Algorithm to save the list for future checkings

40. Algorithm of the inner loop to compute forces and potential energy