Molecular Dynamics

1. Molecular Dynamics A brief overview

2. Notes - Websites • "A Molecular Dynamics Primer", F. Ercolessi http://www.fisica.uniud.it/~ercolessi/md/ • http://cacs.usc.edu/education/cs596.html "Scientific computing and visualization", A. Nakano University of Southern California • "The Art of Molecular Dynamics Simulation", D. C. Rapaport, CUP, 1997 • “Computer simulation of liquids”, M. P. Allen, D. J. Tildesley, OUP, 1990

3. What is molecular dynamics? • Solving the classical equations of motion • For a system of N (N>>3) particles • Which interact through a “given” potential • And then apply some “tricks” … • Deterministic technique  Monte Carlo

4. What is molecular dynamics? • However, errors in trajectories always accumulate: MD is astatistical mechanics method  thermodynamic properties • To obtain set of configurations according to statistical ensemble • Microcanonical ensemble (NVE) • Canonical ensemble (NVT) • Isobaric-Isothermal Ensemble (NPT) • Grand canonical ensemble (also number of particles can change, mVT) • Ergodic hypothesis • Also used for the optimization of structures (simulated annealing)

5. Applications of molecular dynamics • Properties of liquids • Plasma physics • Defects in solids • Fracture • Surface properties • Friction • Molecular clusters • Biomolecules • Dynamics of galaxies • Formation of stellar clusters

6. Overview I • Model • System Hamiltonian • Interaction potentials: bonded and non-bonded interactions • Finite system – infinite system • Integrator Symplecticity? • Statistical ensemble • Collecting results

7. Overview II Collecting data

8. First MD simulation (1957/1959) weeks

9. First MD simulation usingcontinuous potentials (1960) IBM 704

10. First MD simulation usingLennard-Jones potential (1964)

11. Limitations • Use of classical forces  quantum effects • MD only valid if (at 300 K, dt > 6 ps) • Realism of forces • Time and size limitations • Thousands to millions of atoms • Time step dt should be as large as possible while conserving total energy • In general, dt ≈0.01 x fastest behavior of your system Atoms oscillate about once every 10-12 s in a solid  dt≈10-14 s • Total time: picoseconds to hundreds of nanoseconds • Simulation only reliable if simulation time is much longer than relaxation time of quantities of interest • Idem for correlation length

12. Initialisation • Positions • Random positions • Regular pattern, e.g. fcc lattice • Previous run • Velocities • Random velocity or from Maxwell distribution • Previous run • Linked with temperature • No drift condition • Rescale velocities to realize desired temperature

13. Interaction potentials • Origins are quantum mechanical • Easiest model: Lennard-Jones potential • Truncated (but with energy conservation)

14. Interaction potentials • Distinguish between long range  short range interactions • Distinguish between intermolecular  intramolecular forces

15. Interaction potentials • Stretch energy • Bending energy

16. Interaction potentials • Interactions between charge inhomogeneities Approaches • Point charges • Point multipoles • Screened Coulomb interaction

17. Interaction potentials • Multi-body interactions e.g. Tersoff and Brenner potentials

18. Infinite systems Number of interacting pairs increases enormously • Periodic boundary conditions • Minimum image criterion for short range potentials • At most one among all pairs formed by a particle i in the box and the set of all periodic images of another particle j will interact rc Central Simulation box

19. Infinite systems • Ewald method for Coulomb interaction

20. Integrators • How should a good integration scheme look like? • High accuracy (reproduces true trajectories well) • Good stability (conservation of energy) • Time reversible • Robust (allow for large time steps) • Conservation of phase space density (Liouville’s theorem)(symplecticity) • Simple Euler method is not time reversible and not symplectic.

21. Integrators • Verlet algorithm • Positions • Velocities • Properties • Time reversible • Symplectic • Does not suffer from energy drift • But no info on velocity untill the next step is made +

22. Integrators • Comparison between Euler and Verlet Test system consists of 7 Lennard-Jones atoms (Ar) Time step is 10 fs Time step is 1 fs

23. Integrators • Velocity Verlet • Properties • Velocity calculated explicitly • Possible to control the temperature • Stable • Most commonly used algorithm

24. Neighbour lists • Complexity of force calculations ~O(N2) But there are only interactions, often n << N

25. Neighbour lists • Verlet lists • Idea: introduce a list, where particles are included which are located within interaction sphere • Also introduce a “reservoir”, where particles outside Rc are stored, so that unknown particles cannot become neighbors in next steps • Do an update of the list every n steps • Either statically with fixed n • Or dynamically with an update criterion • There exists an optimal Rskin

26. Neighbour lists • Linked lists method ~ O(N) Interacts with atoms in 26 neighbour cells

27. Measuring • Kinetic energy + Potential energy = Total energy • Temperature per degree of freedom • The caloric curve E(T) • Mean square displacement !Periodic boundary conditions • Diffusion coefficient • Correlation functions

28. MD (and MC) as optimization tool Cooling schemes • Simulated annealing • Start at high T, decrease T in small steps (cooling schedule) • Easy to understand & implement • Drawback: might be easily trapped in local minima

29. Parallel strategies • Atom decomposition • Atoms are distributed among processors • All coordinates are exchanged before computing forces • OK for long range interactions • Easy to implement • Force decomposition • Each processor calculates the interactions for certain atom pairs • Spatial decomposition • Subdivides space and assigns each processor a particular subregion • Atoms are allowed to move from one processor to its neighbours • More complex to implement (similar to linked lists)

30. Analyzing the sequential code • md.c and lmd.c "Scientific computing and visualization", A. Nakano University of Southern California • Code description, details at http://cacs.usc.edu/education/cs596.html • “Basic molecular dynamics algorithms” • “Linked-list cell MD algorithm”

31. Guidelines for final report • Test energy conservation as function of dt • Compare speed of md.c and lmd.c • Add possibility to save configuration, and to start from a previous run. Allow temperature rescaling and equilibration. • Study the behavior of the caloric curve E(T) by means of constant energy runs at a fixed density, starting from a crystalline arrangement (r=0.6-0.8, Tmax=1.5). • Insert calculations of the total linear and angular momenta. Check their conservation. • Insert calculation of the mean square displacement. • Remove periodic boundary conditions and study a free cluster.