NEWTONIAN MOLECULAR DYNAMICS

1. SMA5233 Particle Methods and Molecular DynamicsLecture 4: Integration MethodsA/P Chen Yu ZongTel: 6516-6877Email: phacyz@nus.edu.sghttp://bidd.nus.edu.sgRoom 08-14, level 8, S16 National University of Singapore

2. NEWTONIAN MOLECULAR DYNAMICS Fi=force on particle i, m=mass, a=acceleration and U=potential energy • Some important properties of Newton's equation of motion are: • Conservation of energy • Conservation of linear momentum • Conservation of angular momentum • Time reversibility

3. Computational algorithms: Numerical procedure to integrate the differential equation Finite-difference approach: molecular coordinates and velocities at time t=t+Dt, obtained from coordinates and velocities at time t (Dt is the time interval). Taylor expansion: Rewrite in discrete form (where n is at time t and n+1 at time t+Dt): Integration algorithm: implementation of the above equation

4. Runge-Kutta method • This is a very straightforward but less accurate method if restricted to a particular order • Method • Our task is to solve the differential equation: dx/dt = f(t, y), x(t0)= x0 • That is, given the function f(t,y), find the function x(t) such that it passes through the point (t0, x0), whose derivative is the given function. Any set of ordinary differential equations can be cast in this form if we allow x to be a vector of certain dimensions. • For example, a second-order differential equation can be rewritten as a set of two first-order differential equations. Thus it is suffice to just discuss how to solve the above equation, keeping in mind that x could mean (x1, x2, ... , xd).

5. Runge-Kutta method • Clearly, the most obvious scheme to solve the above equation is to replace the differential by finite differences: dt = h dx = x(t+h) - x(t) • We then get the Euler method or first-order Runge-Kutta formula: x(t+h) = x(t) + h f(t, x(t)) + O(h2) • The term first order refers to the fact that the equation is accurate to first order in the small step size h, thus the (local) truncation error is of order h2. The Euler method is not recommended for practical use, because it is less accurate in comparison to other methods and it is not very stable.

6. Runge-Kutta method Improvement of accuracy • The accuracy of the approximation can be improved by evaluating the function f at two points, once at the starting point, and once at the midpoint. This lead to the second-order Runge-Kutta or midpoint method: k1 = h f(t, x(t)) k2 = h f(t+h/2, x(t)+k1/2) x(t + h) = x(t) + k2 + O(h3) • However, the most popular Runge-Kutta formula is the following fourth-order one: k1 = h f(t, x(t)) k2 = h f(t+h/2, x(t)+k1/2) k3 = h f(t+h/2, x(t)+k2/2) k4 = h f(t+h, x(t)+k3) x(t + h) = x(t) + k1/6 + k2/3 + k3/3 + k4/6 + O(h5)

7. Runge-Kutta method • The fourth-order Runge-Kutta method is suitable for many applications for systems involving few degrees of freedom. But it was almost never used for molecular dynamics involving many degrees of freedom. Reasons for this: • Computationally demanding • Four evaluations of the right-hand side of the equations • Lacks the time-reversal symmetry of the Newton equations.

8. VERLET ALGORITHM This algorithm is particularly suited for molecular dynamics. It has been widely used in many areas from simulation of liquids and solids to biological molecules. Method Our problem is to solve the set of Newton's equations: m d2x / dt2 = F where x is the position vector of a particle, m is the mass of this particle, and F is the total force acting on this particle. We'll discuss the form of the forces later on. To get approximate difference formula for derivatives, let's look at the Taylor expansion of the function x(t) with h=dt: x(t + h) = x(t) + h x'(t) + h2/2 x''(t) + h3/6 x'''(t) + O(h4) Replacing h by -h, we get: x(t - h) = x(t) - h x'(t) + h2/2 x''(t) - h3/6 x'''(t)  + O(h4)

9. VERLET ALGORITHM Sum the forward and backward expansions: • Use rn to calculate Fn • Use rn, rn-1 and Fn (step1) to calculate rn+1 Subtract the forward and backward expansions: Propagate velocities

10. VERLET ALGORITHM • Advantages: • 1. Integration does not require the velocities, only position information is taken into account. • 2. Only a single force evaluation per integration cycle. (Force evaluation is the most computationally expensive part in the simulation). • 3. This formulation, which is based on forward and backward expansions, is naturally reversible in time (a property of the equation of motion). • Disadvantages: • The velocities, which are required for energy evaluation are calculated in an approximate manner only through the equation: vn = (rn+1 - rn-1)/2 Dt . (large errors) • Need to know rn+1 to calculate vn

11. Verlet Algorithm Flow Diagram t-Dt t t+Dt r v F Given current position and position at end of previous time step

12. Verlet Algorithm Flow Diagram t-Dt t t+Dt r v F Compute the force at the current position

13. Verlet Algorithm Flow Diagram t-Dt t t+Dt r v F Compute new position from present and previous positions, and present force

14. Verlet Algorithm Flow Diagram t-2Dt t-Dt t t+Dt r v F Advance to next time step, repeat

15. Instantaneous velocity at time t LEAP FROG INTEGRATOR Define: Evaluate velocities at the midpoint of the position evaluations and Vice versa. Where v n+1/2 is the velocity at t+(1/2)Dt • Use rn to calculate Fn • Use Fn and vn-1/2 to calculate v n+1/2 • Use rn and vn+1/2 to calculate r n+1

16. Leap Frog • Advantages: • 1. Improved evaluation of velocities. • 2. Direct evaluation of velocities gives a useful handle for controling the temperature in the simulation. • 3. Reduces the numerical error problem of the Verlet algorithm. Here O(dt1) terms are added to O(dt0) terms. • Disadvantages: • 1. The velocities at time t are still approximate. • 2. Computationally a little more expensive than Verlet.

17. Leapfrog Algorithm Flow Diagram t-Dt t t+Dt r v F Given current position, and velocity at last half-step

18. Leapfrog Algorithm Flow Diagram t-Dt t t+Dt r v F Compute current force

19. Leapfrog Algorithm Flow Diagram t-Dt t t+Dt r v F Compute velocity at next half-step

20. Leapfrog Algorithm Flow Diagram t-Dt t t+Dt r v F Compute next position

21. Leapfrog Algorithm Flow Diagram t-2Dt t-Dt t t+Dt r v F Advance to next time step,repeat

22. VELOCITY VERLET INTEGRATOR Stores positions, velocities and accelerations all at the same time t

23. Velocity Verlet Algorithm Flow Diagram t-Dt t t+Dt r v F Given current position, velocity, and force

24. Velocity Verlet Algorithm Flow Diagram t-Dt t t+Dt r v F Compute new position

25. Velocity Verlet Algorithm Flow Diagram t-Dt t t+Dt r v F Compute velocity at half step

26. Velocity Verlet Algorithm Flow Diagram t-Dt t t+Dt r v F Compute force at new position

27. Velocity Verlet Algorithm Flow Diagram t-Dt t t+Dt r v F Compute velocity at full step

28. Velocity Verlet Algorithm Flow Diagram t-2Dt t-Dt t t+Dt r v F Advance to next time step, repeat

29. Symplectic algorithms • What is symplectic algorithm? • In the Verlet algorithm, our basic dynamic variables are the coordinates of the particles. And we have a somewhat inconvenient situation that the initial conditions are specified as the initial positions and initial velocities, while the difference scheme involves position only. It is possible to reformulate the Newtonian mechanics so that both positions and momenta are the basic dynamic variables. • This formulation is called Hamiltonian dynamics (which is equivalent to Newtonian dynamics). A Hamiltonian system is completely described by the Hamiltonian function (kinetic energy plus potential energy): • H(p,q) = K + V

30. Symplectic algorithms • A Hamiltonian system is completely described by the Hamiltonian function (kinetic energy plus potential energy): • H(p,q) = K + V • together with Hamilton's equations of motion: • dpi/dt = - dH/dqi,    dqi/dt = dH/dpi • where H is considered to be a function of the generalized coordinates qi and generalized momentum pi. The derivative on H is a partial derivative. The index i labels the degree of freedom, not the particle. • Algorithms for solving this set of  Hamilton's equations of motion are symplectic algorithms.

31. Symplectic algorithms • Systems for which a symplectic algorithm can be derived: • The symplectic algorithm can be derived if the Hamiltonian takes a simple form: • H = K + V = sumi pi2/2mi + V(q1, q2, ...) • This is certainly the case if we use Cartesian coordinates and the interactions between particles depend only on the positions of the particles. • With this particular form of Hamiltonian, Hamilton's equations can be written as: • dpi/dt = - dV/dqi = Fi • dqi/dt = dK/dpi = pi/mi • This set of equations reduce to Newton's equations if one eliminates the variable pi.

32. Symplectic algorithms • Algorithm A • We use an Euler type approximation for the momentum derivatives, and then a similar Euler algorithm for the position derivatives. However, for the positions, we use immediately the new result for the momenta. Thus, it is not exactly an Euler algorithm. For brevity, we'll use a vector notation: • p = ( p1, p2, ... ) • q = ( q1, q2, ... ) • F = ( F1, F2, ... ) • and assume that all particles have the same mass m. With these notations and assumption, we have: • p(t + h) = p(t) + hF(q(t)) • q(t + h) = q(t) + h/m p(t + h) • Clearly, the algorithm is accurate to first order in h.

33. Symplectic algorithms • Algorithm B • We do exactly the same thing as in algorithm A, except that we reverse the order of calculation. The new position is found first. • q(t + h) = q(t) + h/m p(t) • p(t + h) = p(t) + hF(q(t + h)) • If you compare algorithms A and B, you will find that it is not exactly the same. It is very interesting to note that both algorithm A and B are in fact mathematically equivalent to the Verlet algorithm if you eliminate the momenta from the equations. However, they have better roundoff error control over the plain Verlet algorithm.

34. Symplectic algorithms • Algorithm C (second-order symplectic algorithm or velocity form of the Verlet algorithm): • This is the most widely used algorithm that we are really interested. • This algorithm is derived by dividing the step size h into two, for the first half step size h/2, we use the algorithm A, for the second half of the step, we use the algorithm B. In total, we have full step size h. That is, use the formula for algorithm A but set h/2. We get equations relating the quantities at time t to time t + h/2. Then we use the equations for algorithms B with step size again h/2, but the time is for t+h/2 to t + h. We get four equations: • p(t + h/2) = p(t) + h/2 F(q(t)) • q(t + h/2) = q(t) + h/2m p(t + h/2) • q(t + h) = q(t + h/2) + h/2m p(t + h/2) • p(t + h) = p(t + h/2) + h/2 F(q(t + h))

35. Symplectic algorithms • Now we eliminate the intermediate quantities evaluated at time t + h/2, we get the final result for the symplectic algorithm: • q(t + h) = q(t) + h/m p(t) + h2/2m F(q(t)) • p(t + h) = p(t) + h/2 [ F(q(t)) + F(q(t + h)) ] • Note that the coordinates q have to be evaluated first since the new values at time t+h are used to evaluate the forces. • The algorithm C is also second-order in step size h. From the derivation it is not obvious that it is necessarily superior than Verlet algorithm. However, the symplectic algorithm has several advantages: • It is time-reversible • The system can be started naturally, with initial position q0 and • initial momentum p0 = mv0 • The symplectic algorithms A, B, and C have one important property that share with the original Hamiltonian system---they preserve Poincare invariants.

36. Symplectic algorithms This last property of Hamiltonian system is usually discussed in advanced mechanics or dynamical system courses. It is about the phase space (p,q) flow property. A solution with initial condition (p0,q0) of the Hamiltonian system p=p(t,p0,q0), q = q(t,p0, q0) can be viewed as a map in the phase space from the point (p0,q0) to the point (p,q), where the time parameter t is taken as some constant. The map is also called a canonical transformation, or symplectic transformation. One of the series of Poincare invariants is the phase space volume. Certain region of points of volume V0 can be mapped to a set of new points forming volume V. Invariance means V=V0. This is called Liouville theorem. Not only is the volume in phase space invariant under the mapping of the dynamics, but also there are other lower dimensional objects (hypersurfaces) which are invariant. The symplectic algorithms preserve this invariant property exactly.

37. Example MD Integration Code Verlet.C: A C++ program that implements the Verlet algorithm to solve the Harmonic Oscillator problem. To run this program, it must first be compiled. To compile this program, name it as a file called "Verlet.C". On most machines, the following command will work to compile this program: "g++ Verlet.C -o Verlet -lm". This command should produce an executable file called "Verlet". This file can be run from the command line: "./Verlet 0.01". This will run the Verlet program with a time-step of 0.01 and output the results to standard output. To plot the results, the output should first be placed into a file. This can be done on Unix-based machines with the following command: "./Verlet 0.01 > Verlet.out". This will place the output into a file "Verlet.out", which can be plotted.

38. Example MD Integration Code #include <iostream> #include <fstream> #include <iomanip> #include <cmath> using namespace std; int main(int argc, char* argv[]) { double alpha = 1.0; // force constant double f; // force double m = 1.0; // mass double t = 0.0; // current time double x = 4.0; // current pos double xold; double xnew; double v = 0.0; // current velocity double dt = 0.1; // timestep

39. Example MD Integration Code double tot_time = 10000.0; // total time unsigned long nstep = (int)(tot_time/dt); //number of integration steps dt = strtod(argv[1],0); nstep = (int)(tot_time/dt); cout << t << " " << x << endl; xold = x - v * dt; for (unsigned long step=0ul; step<nstep; step++) { t = t + dt; f = -alpha * x; // force at n xnew = 2.0*x - xold + f/m * dt * dt; // position at n+1 xold = x; x = xnew; // we could save mem here by reusing xold to // the xnew value cout << t << " " << x << endl; } return 0; }

40. Example MD Integration Code #ifndef __FORCE_h_ #define __FORCE_h_ #include "VEC3D.h" #include "MDUTIL.h" using namespace std; void compute_force(const VEC3D& L, const vector<VEC3D>& x, vector<VEC3D>& f, float& pe); #endif

41. Example MD Integration Code #include "FORCE.h" void compute_force(const VEC3D& L, const vector<VEC3D>& x, vector<VEC3D>& f, float& pe) { static VEC3D dx, dy; static int i, j; static float invDistSq, invDistSix, invDistTwelve; static VEC3D fpart; static int Np; // set number of particles Np Np = x.size(); // zero all the energies/forces pe = 0.0;

42. Example MD Integration Code for (i=0;i<Np;i++) f[i]=VEC3D(0.0); // Loop over all pairs of particles for (i=0;i<Np;i++) for (j=i+1;j<Np;j++) { dx = x[i] - x[j]; // vector separation dy = closest_image(dx,L); // modulo box length: nearest image convention invDistSq = 1.0/mag2(dy); // 1/dy^2 invDistSix = invDistSq*invDistSq*invDistSq; // 1/dy^6 invDistTwelve = invDistSix*invDistSix; // 1/dy^12 pe += invDistTwelve - 2.0*invDistSix; // add contribution to potential energy fpart = 12.0 * (invDistTwelve-invDistSix) * invDistSq * dy; f[i] += fpart; f[j] -= fpart; } };