Transcript Lecture 4: Integration Methods - BIDD

SMA5233
Particle Methods and Molecular Dynamics
Lecture 4: Integration Methods
A/P Chen Yu Zong
Tel: 6516-6877
Email: [email protected]
http://bidd.nus.edu.sg
Room 08-14, level 8, S16
National University of Singapore
NEWTONIAN MOLECULAR DYNAMICS
..
Fi  iU  mi ai mi r i
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
2
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:

1 
1
2
r t  Dt   r t   r t Dt  r t Dt  ...  r t   vt Dt  at Dt 2  ...
2
2
Rewrite in discrete form (where n is at time t and n+1 at time t+Dt):
 
1  Fn  2
rn1  rn  vn Dt   Dt  O Dt 3 ...
2 m 
Integration algorithm: implementation of the above equation
3
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).
4
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.
5
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 fourthorder 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)
6
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.
7
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)
8
VERLET ALGORITHM
 
1  Fn  2
rn1  rn  vn Dt   Dt  O Dt 3 ...
2 m 
 
1  Fn  2
rn1  rn  vn Dt   Dt  O Dt 3 ...
2 m 
Sum the forward and backward expansions:
 
 Fn  2
1.
rn1  2rn  rn1   Dt  O Dt 4 ...
2.
m
Use rn to calculate Fn
Use rn, rn-1 and Fn (step1) to calculate rn+1
Subtract the forward and backward expansions:
 
rn1  rn1
vn 
 O Dt 2
2Dt
Propagate velocities
9
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:
1. 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)
2. Need to know rn+1 to calculate vn
10
Verlet Algorithm Flow Diagram
t-Dt
r
v
F
t
t+Dt
Given current position
and position at end of
previous time step
 
 Fn  2
rn1  2rn  rn1   Dt  O Dt 4 ...
m
11
Verlet Algorithm Flow Diagram
t-Dt
r
v
t
t+Dt
Compute the force at
the current position
F
 
 Fn  2
rn1  2rn  rn1   Dt  O Dt 4 ...
m
12
Verlet Algorithm Flow Diagram
t-Dt
r
v
F
t
t+Dt
Compute new position
from present and previous
positions, and present
force
 
 Fn  2
rn1  2rn  rn1   Dt  O Dt 4 ...
m
13
Verlet Algorithm Flow Diagram
t-2Dt
r
v
t-Dt
t
t+Dt
Advance to next time
step, repeat
F
 
 Fn  2
4


rn1  2rn  rn1   Dt  O Dt ...
m
14
LEAP FROG INTEGRATOR
Define:
vn 1
2
 Dt  r (t )  r (t  Dt ) rn  rn1
 v t   

2
Dt
Dt

 Dt  r (t  Dt )  r (t ) rn1  rn
vn 1  v t   

2
Dt
Dt
2

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
rn1  rn  vn 1 Dt
vn  1
2
2
 Fn 
 vn  1   Dt
2  m
1. Use rn to calculate Fn
2. Use Fn and vn-1/2 to calculate v n+1/2
3. Use rn and vn+1/2 to calculate r n+1
vn


 vn  1  vn  1 
2
2

2
Instantaneous velocity at time t
15
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.
16
Leapfrog Algorithm Flow Diagram
t-Dt
t
r
t+Dt
Given current position,
and velocity at last
half-step
v
F
rn1  rn  vn 1 Dt
2
vn  1
2
 Fn 
 vn  1   Dt
2  m
17
Leapfrog Algorithm Flow Diagram
t-Dt
t
r
t+Dt
Compute current force
v
F
rn1  rn  vn 1 Dt
2
vn  1
2
 Fn 
 vn  1   Dt
2  m
18
Leapfrog Algorithm Flow Diagram
t-Dt
t
r
t+Dt
Compute velocity
at next half-step
v
F
rn1  rn  vn 1 Dt
2
vn  1
2
 Fn 
 vn  1   Dt
2  m
19
Leapfrog Algorithm Flow Diagram
t-Dt
t
r
t+Dt
Compute next
position
v
F
rn1  rn  vn 1 Dt
2
vn  1
2
 Fn 
 vn  1   Dt
2  m
20
Leapfrog Algorithm Flow Diagram
t-2Dt
t-Dt
t
t+Dt
r
Advance to next time
step,repeat
v
F
rn1  rn  vn 1 Dt
2
vn  1
2
 Fn 
 vn  1   Dt
2  m
21
VELOCITY VERLET INTEGRATOR
Stores positions, velocities and accelerations all at the same time t
1  Fn  2
rn1  rn  vn Dt   Dt
2 m 
1  Fn Fn 1  2
Dt
vn 1  vn   
2 m
m 
22
Velocity Verlet Algorithm
Flow Diagram
t-Dt
r
v
t
t+Dt
Given current position,
velocity, and force
F
1  Fn  2
rn1  rn  vn Dt   Dt
2 m 
23
Velocity Verlet Algorithm
Flow Diagram
t-Dt
t
t+Dt
r
Compute new position
v
F
1  Fn  2
rn1  rn  vn Dt   Dt
2 m 
24
Velocity Verlet Algorithm
Flow Diagram
t-Dt
r
v
t
t+Dt
Compute velocity at
half step
F
25
Velocity Verlet Algorithm
Flow Diagram
t-Dt
r
v
t
t+Dt
Compute force at new
position
F
26
Velocity Verlet Algorithm
Flow Diagram
t-Dt
r
v
t
t+Dt
Compute velocity at
full step
F
27
Velocity Verlet Algorithm
Flow Diagram
t-2Dt
r
v
t-Dt
t
t+Dt
Advance to next time
step, repeat
F
28
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
29
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.
30
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.
31
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.
32
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.
33
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))
34
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
35
Poincare invariants.
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.
36
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.
37
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
38
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;
}
39
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
40
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;
41
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;
}
};
42