Introduction into Molecular Dynamics Ralf Schneider, Abha Rai, Amit Rai Sharma

1. Introduction into Molecular DynamicsRalf Schneider, Abha Rai, Amit Rai Sharma Outline: 1. Basics 2. Potentials 3. History 4. Numerics 5. Analysis of MD runs (6. Physics extensions= 7. Numerical extensions 8. Summary

2. 1. Molecular Dynamics • Solve Newton’s equation for a molecular system: • Or equivalently solve the classical Hamiltonian equation:

3. 1. Molecular Dynamics method • deterministicmethod: state of the system at any future time can be predicted from its current state • MD cycle for one step: 1) force acting on each atom is assumed to be constant during the time interval 2) forces on the atoms are computed and combined with the current positions and velocities to generate new positions and velocities a short time ahead

4. 1. Molecular Dynamics method K. Nordlund U. Helsinki

5. 1. Molecular Dynamics method K. Nordlund, U. Helsinki

6. 1. Molecular Dynamics method K. Nordlund U. Helsinki

7. 1. Molecular Dynamics method Important for converting CO2 to organic forms of carbon and in the photosynthetic process. Even though the pocket is closed, a CO2 molecule escapes, which was a surprise. RuBisCO protein simulations Paul Crozier, Sandia

8. 1. Molecular Dynamics method ATP Synthase: within the mitochondria of a cell a rotary engine uses the potential difference across the bilipid layer to power a chemical transformation of ADP into ATP H. Wang U. California, Santa Cruz

9. 1. Motivation: Why atomistic MD simulations? • MD simulations provide a molecular level picture of structure and dynamics (biological systems!)  property/structure relationships • Experiments often do not provide the molecular level information available from simulations • Simulators and experimentalists can have a synergistic relationship, leading to new insights into materials properties

10. 1. Motivation: Why atomistic MD simulations? MD simulations allow prediction of properties for • Novel materials which have not been synthesized • Existing materials whose properties are difficult to measure or poorly understood • Model validation

11. 1. Timescales Bond vibrations - 1 fs Collective vibrations - 1 ps Conformational transitions - ps or longer Enzyme catalysis - microsecond/millisecond Ligand Binding - micro/millisecond Protein Folding - millisecond/second Molecular dynamics: Integration timestep - 1 femtosecond Set by fastest varying force. Accessible timescale about 10 nanoseconds.

12. 1. MD dynamics

13. 2. Molecular Dynamics: Potential We need to know The motion of the atoms in a molecule, x(t) and therefore, the potential energy, V(x)

14. 2. Molecular Dynamics: Potentials How do we describe the potential energy V(x) for a molecule? Potential Energy includes terms for Bond stretching Angle Bending Torsional rotation Improper dihedrals

15. 2. Molecular Dynamics: Potentials Potential energy includes terms for (contd.) Electrostatic Interactions van der Waals Interactions

16. 2. Molecular Interaction Types – Non-bonded Energy Terms • Lennard-Jones Energy. • Coloumb Energy.

17. 2. Molecular Interaction Types – Bonded Energy Terms • Bond energy: • Bond Angle Energy:

18. 2. Molecular Interaction Types – Bonded Energy Terms • Improper Dihedral Angle Energy: • Dihedral Angle Energy:

19. 2. Scaling • Scaling by model parameters • size s • energy e • mass m taken from Dr. D. A. Kofke’s lectures on Molecular Simulation, SUNY Buffalo http://www.eng.buffalo.edu/~kofke/ce530/index.html

20. 2. L-J: dimensionless form • Dimensions and units - scaling • Lennard-Jones potential in dimensionless form • Dimensionless properties must also be parameter independent • convenient to report properties in this form, e.g. P*(r*) • select model values to get actual values of properties • Equivalent to selecting unit value for parameters taken from Dr. D. A. Kofke’s lectures on Molecular Simulation, SUNY Buffalo http://www.eng.buffalo.edu/~kofke/ce530/index.html

21. 3. Historical Perspective on MD

22. 3. First molecular dynamics simulation (1957/59) Hard disks and spheres (calculation of collision times) IBM-704: solid phase liquid phase liquid-vapour-phase N=32  6.5x105 coll.  4 days N=500  107 coll.  2.3 years N=32: 7000 collisions / h N=500: 500 collisions / h Production run ~20000 steps

23. 3. First MD simulation using continuous potentials (1964) CDC-3600 VACF RDF MSD 864 particles Time / Step ~ 45s Production run ~20000 steps  10 days ! (standard PC [Pentium 1.2 GHz]: ½ hour)

24. 3. MD – development (aus T. Schlick, „Molecular Modelling and Simulation“, 2002)

25. 4. Verlet algorithm then Let Starting from and all subsequent positions are determined For the kinetic energy we need the velocities • Note: the velcoities are one step behind. Therefore: • Specify positions and • Compute the forces at timestep n: • 3. Compute the positions at timestep (n+1) as in (1.1): • 4. Compute velocities at timestep n as in (1.2); then increment n and goto 2.

26. 4. A widely-used algorithm: Leap-frog Verlet • Using accelerations of the current time step, compute the velocities at half-time step: v(t+Dt/2) = v(t – Dt/2) + a(t)Dt v t-Dt/2 t t+Dt/2 t+Dt t+3Dt/2 t+2Dt

27. 4. A widely-used algorithm: Leap-frog Verlet • Using accelerations of the current time step, compute the velocities at half-time step: v(t+Dt/2) = v(t – Dt/2) + a(t)Dt • Then determine positions at the next time step: X(t+Dt) = X(t) + v(t + Dt/2)Dt v X t-Dt/2 t t+Dt/2 t+Dt t+3Dt/2 t+2Dt

28. 4. A widely-used algorithm: Leap-frog Verlet • Using accelerations of the current time step, compute the velocities at half-time step: v(t+Dt/2) = v(t – Dt/2) + a(t)Dt • Then determine positions at the next time step: X(t+Dt) = X(t) + v(t + Dt/2)Dt v v X t-Dt/2 t t+Dt/2 t+Dt t+3Dt/2 t+2Dt

29. 4. Verlet algorithm- velocity form

30. 4. Advantages • Positions and velocities are now in step => kinetic and potential energies are in step • Numerical stability is enhanced • Eq (1.2) gives velocity as difference of 2 r’s of the same order of magnitude => round-off errors • important for long runs • With reasonable h, Verlet’s algorithm conserves energy

31. 4. Beeman algorithm • Better velocities, better energy conservation • More expensive to calculate

32. 4. Predictor-corrector algorithms • 1. Predictor. From the positions and their time derivatives up to a certain order q, all known at time t, one ``predicts'' the same quantities at time by means of a Taylor expansion. Among these quantities are, of course, accelerations . • 2. Force evaluation. The force is computed taking the gradient of the potential at the predicted positions. The resulting acceleration will be in general different from the ``predicted acceleration''. The difference between the two constitutes an ``error signal''. • 3. Corrector. This error signal is used to ``correct'' positions and their derivatives. All the corrections are proportional to the error signal, the coefficient of proportionality being a ``magic number'' determined to maximize the stability of the algorithm. Fifth-order Gear (requires more calculational effort and memory than Verlet, but needs only one calculation of the force per time step, wheras Verlet needs two!

33. 4. Evaluate integration methods • Fast, minimal memory, easy to program • Calculation of force is time consuming • Conservation of energy and momentum • Time-reversible • Long time step can be used

34. 4. Choosing the time step • Too small: covering small conformation space • Too large: instability • Suggested time steps • Translation, 10 fs • Flexible molecules and rigid bonds, 2fs • Flexible molecules and bonds, 1fs

35. 4. How do you run a MD simulation? • Get the initial configuration • Assign initial velocities At thermal equilibrium, the expected value of the kinetic energy of the system at temperature T is: This can be obtained by assigning the velocity components vi from a random Gaussian distribution with mean 0 and standard deviation (kBT/mi):

36. 4. Periodic Boundary Conditions • infinite system with small number of particles • remove surface effects • shaded box represents the system we are simulating, while the surrounding boxes are exact copies in every detail • whenever an atom leaves the simulation cell, it is replaced by another with exactly the same velocity, entering from the opposite cell face (number of atoms in the cell is conserved) • rcut is the cutoff radius when calculating the force between two atoms

37. 4. Minimum Image These two are same distance from central atom, yet: Black atom interactsblue atom does not • Bulk system modeled via periodic boundary condition • not feasible to include interactions with all images • must truncate potential at half the box length (at most) to have all separations treated consistently • Contributions from distant separations may be important Only interactions considered

38. 4. Potential cut-offs Bonded interactions: local, therefore O(N), where N is the number of atoms in the molecule considered) Non-bonded interactions: involve all pairs of Atoms, therefore O(N2) Reducing the computing time: use of cut-off in UNB The cutoff distance may be no greater than ½ L (L= box length)

39. 4. Potential truncation common approach: cut-off the at a fixed value Rcut problem: discontinuity in energy and force possibility of large errors

40. 4. Speed-up Tamar Schlick, “Molecular Modeling and Simulation”, Springer

41. 4. Cutoff schemes for faster energy computation wij: weights(0< wij <1). Can be used to exclude bonded terms, or to scale some interactions (usually 1-4) S(r) :cut-off function. Three types: 1)Truncation: b

42. 4. Cutoff schemes for faster energy computation 2. Switching a b with 3. Shifting or b

43. 4. Neighbor lists • Verlet requires O(N) operations • Force needs O(N2) operations at each step • Most of these are outside range and hence zero • Time reduced by counting only those within range listed in a table (needs to be updated) • - Verlet (1967) suggested keeping a list of the near neighbors for a particular molecule, which is updated periodically • - Between updates of the list, the program does not check through all the molecules, just those on the list, to calculate distances and minimum images, check distances against cutoff, etc.

44. 4. Simulating at constant T: the Berendsen scheme system • Bath supplies or removes heat from the system as appropriate • Exponentially scale the velocities at each time step by the factor : where  determines how strong the bath influences the system Heat bath T: “kinetic” temperature Berendsen et al. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81:3684 (1984)

45. system 4. Simulating at constant P: Berendsen scheme • Couple the system to a pressure bath • Exponentially scale the volume of the simulation box at each time step by a factor : where  : isothermal compressibility P : coupling constant pressure bath where u : volume xi: position of particle i Fi : force on particle i Berendsen et al. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81:3684 (1984)

46. 4. MD scheme

47. 5. Analysis of MD Configurations Averages Fluctuations Time Correlations

48. 5. Time averages and ensemble averages macroscopic numbers of atoms or molecules (of the order of 1023,Avogadro's number is 6.02214199 × 1023 ): impossible to handle for MD statistical mechanics (Boltzmann, Gibbs): a single system evolving in time is replaced by a large number of replications of the same system that are considered simultaneously time average is replaced by an ensemble average:

49. 5. Ergodic hypothesis • Classical statistical mechanics integrates over all of phase space {r,p}. • The ergodic hypothesis assumes that for sufficiently long time the phase trajectory of a closed system passes arbitrarily close to every point in phase space. • Thus the two averages are equal:

50. 5. Statistical Mechanics Extracting properties from simulations • static propertiessuch as structure, energy, and pressure are obtained from pair (radial) distribution functions DME-water and DMP-water solutions