Dynamics Molecular I: Principles
Basics of the molecular-dynamics (MD) method1-3 are described, along with corresponding data structures in program, md.c.
Newtons Second Law of Motion
TRAJECTORY, COORDINATE, AND ACCELERATION Physical system = a set of atomic coordinates: r {ri = (xi ,yi ,zi ) | xi ,yi ,zi " #,i = 0,...,N $ 1} ,
where ! is the set of r numbers (in the program, represented by a double precision variable) and we real use a vector notation, ri . (Data strcutures in md.c) int nAtom: N, the number of atoms. NMAX: Maximum number of atoms that can be handled by the program. double r[NMAX][3]: r[i][0], r[i][1], and r[i][2] are the x, y, and z coordinates of the i-th atom, where i = 0, .., N-1.
z 1 y 2
N
x
r Trajectory: A mapping from time to a point in the 3-dimensional space, t " # a ri (t) " # 3 . In fact, a trajectory of an N-atom system is regarded as a curve in 3N-dimensional space. A point on the r curve is then specified by a 3N-element vector, r N = (x 0 , y0 ,z0 , x1 , y1 ,z1 ,..., x N!1, yN !1,zN !1 ).
vi(t=0) ri(t=0)
vi(t=t1) ri(t=t1)
Velocity: Short-time limit of an average speed (how fast and in which direction the particle is moving), r r r r r (t) = dr " lim ri (t + #) % ri (t) . v i (t) = ri dt # $0 #
rv[i][0], rv[i][1], and rv[i][2] are the x, y, and z components of the r velocity vector, vi , of the i-th atom. Acceleration: Rate ! which a velocity changes (whether the particle is accelerating or decelerating), at double rv[NMAX][3]:
r r r 2r r r (t) = d r = dvi " lim vi (t + #) % vi (t) . ai (t) = ri dt #$0 # dt 2
1
!
ra[i][0], ra[i][1], and ra[i][2] are the x, y, and z components of the r acceleration vector, ai , of the i-th atom. double ra[NMAX][3]:
vi(t+!) vi(t) vi(t)
ai(t)! vi(t+!)
The acceleration of a particle can be estimated from three consecutive positions on its trajectory separated by small time increments: r r r v i (t + " /2) $ v i (t $ " /2) ai (t) = lim " #0 " r r r r ri (t + ") $ ri (t) ri (t) $ ri (t $ ") $ " " = lim " #0 " r r r ri (t + ") $ 2 ri (t) + ri (t $ ") = lim " #0 "2 NEWTONS EQUATION Newtons equation states that the acceleration of a particle is proportional to the force acting on the ! particle, r r m(t) = Fi (t), ri where m is the mass of the particle. For a heavier particle, the same force causes less acceleration (or deceleration). r r Initial value problem: Given initial particle positions and velocities, {( ri (0),vi (0)) i = 0,...,N " 1} , r r obtain those at later times, {( ri (t),vi (t)) i = 0,...,N " 1;t # 0} . Note that both positions and velocities must be specified in order to predict future trajectories. ! Potential Energy
r ! We consider forces which are derived from a potential energy, V( r N ), which is a function of all atomic positions.
r $ #V #V #V ' r # Fk = " r V ( r N ) = "& , , ). #rk % #x k #y k #zk (
Here partial derivative is defined as
"V V (x 0 , y 0 ,z0 ,..., x k + h, y k ,zk ,...x N $1, y N $1,zN $1 ) $ V (x 0 , y 0 ,z0 ,..., x k , y k ,zk ,...x N $1, y N $1,zN $1 ) = lim , ! "x k h#0 h
i.e., what is the rate of change in V when we slightly change one coordinate of an atom while keeping all the other coordinates fixed.
!
LENNARD-JONES POTENTIAL To model certain materials such as neon and argon liquids, the Lennard-Jones potential is often used by scientists: 2
N #2 r V ( r N ) = " u(rij ) = " i< j i= 0
N #1
" u(| r
j= i+1
r
ij
|)
rrr where rij = ri ! rj is the relative position vector between atoms i and j, and
!
+$ # '12 $ # ' 6 . u(r) = 4"-& ) * & ) 0 . %r( / 0 ,% r (
Physical meaning of the potential: An atom consists of a nucleus and electrons surrounding it. Two atoms interact with each other in the following ways. > Short-range repulsion! (the first term): the Pauli exclusion principle states that electrons cannot occupy the same position, resulting in repulsion between the atoms. Note that for r 0, this term is dominant.
> Long-range attraction (the second term): Electrons around nuclei polarize (change distribution), creating electrostatic attraction between atoms. Note that for r , this term becomes dominant.
For Argon atoms, > m = 6.6 10-23 gram > = 1.66 10-14 erg (erg = gramcm2/(second2) is a unit of energy) > = 3.4 10-8 cm
* Dont program with these numbers (in the CGS unit); it will cause floating-point under- or overflows. NORMALIZATION
r We introduce normalized coordinates, ri ! , potential energy, V', and time, t', that are dimensionless and are defined as follows. rr r ' ri = ri# = 3.4 $10%8 [cm] $ ri" " ) (V = V "& = 1.66 $10%14 [erg] $ V " ) %12 * t = # m /& t " = 2.2 $10 [sec] $ t "
With these substitutions, Newtons equation becomes r " d 2 ri$ " &V $ m # 2 =% r. ! 2 m# dt $ # &ri$ * In the derivation above, note that r r r r d 2 ri ri (t + ") $ 2 ri (t) + ri (t $ ") m 2 = m lim ! " #0 dt "2 r r r % [ ri&(t + ") $ 2 ri&(t) + ri&(t $ ")] = m lim 2 2 " #0 % m /' (" &)
[
]
3
!
In summary, the normalized Newtons equation for atoms interacting with the Lennard-Jones potential is given below. (From now on, we always use normalized variables, and omit primes.) r d 2 ri "V r = ! r = ai 2 dt "ri rN V( r ) = # u(rij )
i<j
$1 1' u(r) = 4& 1 2 ! 6 ) %r r(
The figure in the right shows the normalized Lennard-Jones potential, u(r), as a function of interatomic distance, r. The distance, r0, at which the potential takes its minimum value, u(r0), is obtained as follows.
du = 0 ! r0 = 21/ 6 " 1.12 dr r =r0 $1 1' u(r0 ) = 4& 12/ 6 # 6 / 6 ) = #1 %2 2(
0
u(r)/!
-0.5
-1
21/6 ~ 1.12 0 0.5 1 1.5 r/" 2 2.5 3
Figure: Lennard-Jones potential.
ANALYTIC FORMULA FOR FORCES
r # ak = " r $ u(rij ) #rk i< j
#r du = "$ rij i< j #rk drij
Lets evaluate the two factors in the above expression:
"rij # " " " & ,( r =% !, "rk $ "x k "y k "zk '
=
i j i 2
(x
j j
i
) x j ) + ( y i ) y j ) + (zi ) z j )
i j
2
2
2
1.
r rij = (*ik ) * jk ) rij
(2(x ) x ),2(y ) y ),2(z ) z )) * ( 2 ( x ) x ) + ( y ) y ) + (z ) z )
2 2 i j i i j
ik
) * jk )
#d 1 1/ 2 )1/ 2 df & %Q [ f (x)] = [ f (x)] ( $ dx 2 dx ' # 12 6 & du 48 # 1 1& = 4% " 13 + 7 ( = " % 12 " 6 ( $r dr r' r $r 2r ' # d "n "n"1 & %Q r = "nr ( $ dr '
2.
!
In the first result,
!
4
#1 (i = k) !ik = $ %0 (i " k)
is called Kroneckers delta expression.
r Substituting these two results back into the expression for ak , we obtain
r r " 1 du % ak = ( rij $ ! ' () ! ) jk ) # r dr & r= rij ik i< j
where
!
1 du 48 " 1 1% = 2 $ 12 ! 6 ' r dr r # r 2r &
Summary of Molecular-Dynamics Equation System
r r Given initial atomic positions and velocities, {( ri (0),vi (0)) i = 0,...,N " 1} , obtain those at later times, r r {( ri (t),vi (t)) i = 0,...,N " 1;t # 0} , by integrating the following differential equation:
!
r r r (t) = a (t) = " # $ u( r ) =$ r (t)%" 1 du ( rk ' * (+ " + jk ) !r k ij ij & r dr ) r= rij (t ) ik #rk i< j i< j
where
" 1 du 48 # 1 1& = 2 % 12 " 6 ( r dr r $ r 2r '
! The sum over i and j to evaluate the acceleration is implemented in a program as follows: r for i = 0 to N-1, ai = 0 ! for i = 0 to N-2 for j = i+1 to N-1 r r $ 1 du ' compute " = rij &# ) r % r dr ( r= rij r r ai + = ! r r aj ! = "
! Discretization
We need to discretize the trajectories in order to solve the problem on a computer: Instead of considering r r (ri (t), vi (t)) for t ! 0 for continuous time, we consider a sequence of states r r r r r r (ri (0), vi (0)) a (ri (!), vi (!)) a (ri (2!), vi (2!)) aL
ri(!) ri(2!)
ri(0)
r r r r The question is: How to predict the next state, (ri (t + !), vi (t + !)) , from the current state, (ri (t), vi (t)) ?
5
double DeltaT:
in md.c.
VERLET DISCRETIZATION Lets apply the Taylor expansion
f (x 0 + h) = #
hn dn f n n=0 n! dx
"
= f (x 0 ) + h
x= x 0
df dx
+
x= x 0
h2 d2 f 2 dx 2
+
x= x 0
h3 d3 f 3! dx 3
+ L,
x= x 0
which predicts a function value at a distant point, x0+h, from derivatives at x0, to an atomic coordinate.
!
! To obtain velocities, !
r r r 1r 1r ri (t + ") = ri (t) + v i (t)" + ai (t)"2 + (t)"3 + O("4 ) ri 2 6 r r r 1r 1r + ri (t # ") = ri (t) # v i (t)" + ai (t)"2 # (t)"3 + O("4 ) ri 2 6 r r r r ri (t + ") + ri (t # ") = 2 ri (t) + ai (t)"2 + O("4 ) r r r r " ri (t + #) = 2 ri (t) $ ri (t $ #) + ai (t)#2 + O(#4 )
r r r 1r 1r ri (t + ") = ri (t) + v i (t)" + ai (t)"2 + (t)"3 + O("4 ) ri 2 6 r r r 1r 1r # ri (t # ") = ri (t) # v i (t)" + ai (t)"2 # (t)"3 + O("4 ) ri 2 6 r r r ri (t + ") # ri (t # ") = 2v i (t)" + O("3 ) r r r ri (t + #) $ ri (t $ #) ! " vi (t) = + O(#2 ) 2# ! The above two results constitute the Verlet discretization scheme.
!
(Verlet discretization scheme) ! r r r r $ r (t + ") = 2 r (t) # r (t # ") + a (t)"2 + O("4 ) (1) i i i &i r r % r ri (t + ") # ri (t # ") vi (t) = + O("2 ) (2) & ' 2" The Verlet algorithm for MD integration is a straight-forward application of these two equations. (Verlet algorithm) ! r r Given ri (t ! ") and ri (t), r r 1. Compute ai (t) as a function of {ri (t)}, r r r r 2. ri (t + !) " 2 ri (t) # ri (t # !) + ai (t)!2 , r r r 3. vi (t) " [ ri (t + #) $ ri (t $ #)] (2# ) , A problem of this discretization scheme is that the velocity at time t cannot be calculated until the coordinate at time t+ is calculated. In the slight variation, called the velocity-Verlet scheme, of the r r r r ! above discretization scheme, we are given (ri (t), vi (t)) and predict (ri (t + !), vi (t + !)) .
6
VELOCITY-VERLET DISCRETIZATION
r r (Theorem) The following algebraic equation gives the same sequence of states, (ri (n!), vi (n!)) , as that obtained by the Verlet discretization.
#r r r 1r 2 (3) % ri (t + ") = ri (t) + v i (t)" + 2 ai (t)" r r $ r r %v i (t + ") = v i (t) + ai (t) + ai (t + ") " (4) & 2
(Proof) 1. Using Eq. (3) of the velocity-Verlet discretization scheme for two consecutive time steps, t and t+, ! r r r 1r ri (t + 2") = ri (t + ") + v i (t + ")" + ai (t + ")"2 2 r r r 1r # ri (t + ") = ri (t) + v i (t)" + ai (t)"2 2 r r r r r r r r ai (t + ") # ai (t) 2 ri (t + 2") # ri (t + ") = ri (t + ") # ri (t) + [v i (t + ") # v i (t)]" + " 2
! Lets eliminate the velocities using Eq. (4),
r r r r $ ai (t + ") + ai (t) ' r r r r ai (t + ") # ai (t) 2 ri (t + 2") # ri (t + ") = ri (t + ") # ri (t) + & ")" + " % ( 2 2
!
i.e.,
!
r r r r r ri (t + 2!) " ri (t + !) = ri (t + !) " ri (t) + ai (t + !)!2 r r r r ! ri (t + 2") = 2ri (t + ") # ri (t) + ai (t + ")"2
This is nothing but the Verlet rule, Eq. (1), applied to t+2 instead of t+. We have thus derived the first equation of the Verlet scheme from the two equations in the velocity-Verlet scheme. 2. Using Eq. (3)
r r r 1r ri (t + !) = ri (t) + vi (t)! + ai (t)!2 2 r r r 1r + ri (t) = ri (t " !) + vi (t " !)! + ai (t " !)!2 2 r r r r r r ai (t) + ai (t # ") 2 ri (t + ") = ri (t # ") + [vi (t) + vi (t # ")]" + " 2
By multiplying to Eq. (4) with the time shifted by -, r r r r ai (t) + ai (t " #) 2 0 = [vi (t) " vi (t " #)]# " # ! 2 By adding Eqs. (5) and (6),
(5)
(6)
! and this is equivalent to the second equation, Eq. (2), of the Verlet scheme. In summary, by using only the two equations in the velocity-Verlet scheme, Eqs. (3) and (4), we have derived both equations, Eqs. (1) and (2), of the Verlet scheme. Thus the two schemes are equivalent, given atomic positions at the first two time steps.//
7
r r r ri (t + !) = ri (t " !) + 2 vi (t)!
VELOCITY-VERLET ALGORITHM For notational convenience, lets define r ! r !r vi (t + ) " vi (t) + ai (t). 2 2 Eqs. (3) and (4) are then rewritten as
"r r r ! (8) $ ri (t + !) = ri (t) + v i (t + )! $ 2 # r r r $v (t + !) = v (t + ! ) + ! a (t + !) (9) i $i % 2 2i
(7)
The velocity-Verlet algorithm is just a re-statement of Eqs. (7)(9). (Velocity-Verlet algorithm) r r Given (ri (t), vi (t)) , r r 1. Compute ai (t) as a function of {ri (t)}, r r " "r 2. v i (t + ) # v i (t) + ai (t), 2 2 r r r ! 3. ri (t + !) " ri (t) + v i (t + )! , 2 r r 4. Compute ai (t + !) as a function of {ri (t + !)}, ! r r ! !r 5. vi (t + !) " vi (t + ) + ai (t + !) 2 2 r r Computing ai (t) as a function of {ri (t)} is done in function computeAccel(). In the program, we perform an MD simulation for StepLimit steps, and computeAccel() is called StepLimit+1 times. (Velocity-Verlet algorithm for StepLimit steps) rr ! ! Initialize (ri , v i ) for all i r r Compute ai as a function of {ri} for all ifunction computeAccel() for stepCount = 1 to StepLimit do the followingfunction singleStep() r rr vi ! vi + ai" / 2 for all i r rr ri ! ri + vi " for all i r r Compute ai as a function of {ri} for all ifunction computeAccel() r rr vi ! vi + ai" / 2 for all i endfor Now we have learned the skeleton of md.c program. The next lecture will deal with some details of the program. These include: Periodic boundary conditions Shifted potential Energy conservation and temperature Initialization: fcc lattice, random velocity * The velocity-Verlet algorithm is more elegantly derived using a split-operator technique in the Ph.D. thesis by M. Tuckerman at Columbia University. See Reversible multiple time scale molecular dynamics, M. Tuckerman, B. J. Berne, and G. J. Martyna, J. Chem. Phys. 97, 1990 (1992). 8
Molecular Dynamics II: Implementation Details
Periodic Boundary Conditions
(There are only two places you see this in the program) This technique allows us to simulate bulk solid and liquid properties with a small number of atoms (N < 1,000) by eliminating surface effects. Note that for 1,000 atoms arranged in a 10 10 10 cube, 103 83 = 488 atoms (nearly half the atoms) are at the surface. Since these surface atoms are in very different environment from the bulk atoms, the average properties computed for this system are very different from bulk properties. Periodic boundary conditions are implemented by replicating a simulation box of size Lx Ly Lz to form an infinite lattice. Note that we assume atoms to be contained in a rectangular simulation box.
double Region[3]:
Three-element array containing the box lengths, Lx, Ly, and Lz
Ly
Lx
Figure: Periodically repeated images of an original simulation box (shaded) in 2 dimensions. As an atom moves in the original simulation box, all the images move in a concerted manner by the same amount. Since all the images are just shifted copies of an original atom, we need to keep only the coordinates of the original (central) image as a representative of all images. When an atom leaves the central box by crossing a boundary, attention is switched to the image just entering the central box.
Ly
Lx
Figure: Change of the representative atom due to boundary crossing of the atom. 9
At every MD step, we enforce that an atomic coordinate satisfies 0 xi < Lx, 0 yi < Ly, 0 zi < Lz. r r If ri (t) is in the central box and time increment is small, ri (t + !) is at most in the neighbor image (if this assumption becomes invalid, then exceptions may occur during program execution): -Lx xi(t+) < 2Lx, etc. In this case, an atom can be pulled back to the central box by
* $ Lx ' $ Lx ' , x i " x i # SignR& , x i ) # SignR& , x i # Lx ) %2 ( %2 ( , $ Ly ' $ Ly ' , + y i " y i # SignR& , y i ) # SignR& , y i # Ly ) %2 ( %2 ( , $ Lz ' $ Lz ' , , zi " zi # SignR& 2 ,zi ) # SignR& 2 ,zi # Lz ) % ( % ( where
!
and this means
* Lx " Lx % , 2 SignR$ , x i ' = + # 2 & ,( Lx -2
$ & Lx < x i & %0 < x i # Lx & & x #0 i ' xi "
xi > 0 xi ) 0
!
Lx Lx " = x i " Lx 2 2 L L xi " x + x = xi 2 2 Lx Lx xi + + = x i + Lx 2 2
double RegionH[3]:
An array (Lx/2, Ly/2, Lz/2)
Minimum Image Convention
! In principle, an atom interacts with all the images of another atom (and even with its own images except for itself).
In practice, we make the simulation box larger enough so that u(r > Lx/2) etc. can be neglected. Then an atom interacts with only the nearest image of another atom.
i
Ly
j
Lx
10
Or we can imagine a box of size Lx Ly Lz centered at atom i, and this atom interacts only with other atoms in this imaginary box. Therefore,
$ Lx Lx &" 2 # xij < 2 &L Ly &y # yij < %" 2 &2 Lz Lz &" # z < ij &2 2 '
during the computation of forces. In the program, this is achieved by
* $ Lx $ Lx Lx ' Lx ' , xij " x!# SignR& , xij + ) # SignR& , xi # ) ij %2 %2 2( 2( , $ Ly Ly ' , $L L' + yij " yij # SignR& , yij + ) # SignR& x , yij # x ) %2 2( 2( %2 , $ Lz $ Lz , Lz ' Lz ' , zij " zij # SignR& 2 ,zij + 2 ) # SignR& 2 ,zij # 2 ) % ( % ( -
Shifted Potential
With the minimal image convention, interatomic interaction is cut-off at half the box length. If the ! box is not large enough, this truncation causes a significant in discontinuity the potential function and our derivation of forces breaks down (we cannot differentiate the potential).
u(r)
Potential energy jump
Lx/2
r
To remedy this, we introduce a cut-off length rc < Lx/2 etc. and modify the potential such that both u(r) and du/dr are continuous at rc.
% du ' u(r) # u(rc ) # (r # rc ) u"(r) = & dr r=rc ' 0 (
r < rc r $ rc
For the Lennard-Jones potential, rc = 2.5 is conventionally used. This causes some shift in the potential minimum but little shift in the minimum position. ! 11
0
u(R)/!
u(R) u'(R) -0.5
Rc
-1
21/6 ~ 1.12 0 0.5 1 1.5 R/" 2 2.5 3
Figure: Truncated Lennard-Jones potential.
Energy
In addition to the potential energy, V, introduced before, we introduce the kinetic energy, K, and the total energy, E, as a sum of both.
E =K +V
N "1
=#
i=0
mr2 ri + # u(rij ) 2 i< j
where
!
r2 ri = xi2 + yi2 + zi2 .
ENERGY CONSERVATION We can prove that the total energy, E, does not change in time, if Newtons equation is exactly ! integrated. In order to prove the energy conservation, let us consider
E =K +V
N "1% dx $V d2 dy $V dzi $V ( xi + yi2 + zi2 + # ' i +i + * dt $yi dt $zi ) i=0 2 dt i=0 & dt $xi N "1 m N "1% $V $V $V ( x y z = # 2( xi i + yi i + zi i ) + # ' xi + yi + zi * $yi $zi ) i=0 2 i=0 & $xi N "1 m
=#
(
)
N "1 r r N "1 r $V r = # mri + # ri r i $ri i=0 i=0 N "1 r rr = # ri m " Fi ri i=0
(
)
From Newtons equation, the last expression is zero. In the derivation, vector inner product is defined as ! 12
rr a b = (ax ,ay ,az ) (bx ,by ,bz ) = ax bx + ay by + azbz
and note that
d ( fg ) = fg + f g . dt
With discretization approximation, the total energy is not rigorously conserved anymore. How well the total energy is conserved is a good measure for how small we should use. = 0.005 in the normalized unit is typically used for the Lennard-Jones potential (see the figure below).
Temperature
Temperature, T, is related to the kinetic energy as follows:
3Nk B T N "1 m r 2 = # ri , 2 2 i= 0
where kB = 1.38062 10-16 erg/KelvinBoltzmann constant ! By using the normalized coordinates,
3Nk B T # N %1 m r$ 2 = m" 2 2 & ri 2 " m i= 0 2
!
13
"
where /kB = 120 Kelvin
T 1 N %1 r$ 2 = & ri # /kB 3N i= 0
! is the Lennard-Jones temperature unit. So the average temperature 0.5 means 60 K << 273 K = 0C.
Initialization
The program starts to initialize the coordinates and velocities of all atoms. FACE CENTERED CUBIC (FCC) LATTICE Atoms are initially arranged to form a regular lattice. Namely they occupy all the corners and face centers of a cube called a unit cell.
cz cy cx
double gap[3]:
(cx, cy, cz), the edge lengths of the unit cell
Since our unit cell is cubic, cx = cy = cz = c. Each unit cell contains
1 1 ! 8(corners) + ! 6(faces) = 4atoms , 8 2
c = (4 / ! )
1 /3
and the number density of atoms is given by ! = 4 / c 3 , or .
The four coordinates are given (in unit of c) as: (0, 0, 0), (0, 1/2, 1/2), (1/2, 0, 1/2), (1/2, 1/2, 0). These four coordinates are stored in
double origAtom[4][3]
The simulated system is constructed by repeating the unit cells.
int InitUcell[3]:
Stores the number of unit cells in the x, y, and z directions N = 4 InitUcell[0] InitUcell[1] InitUcell[2]
The total number of atoms is therefore given by RANDOM VELOCITY We generate random velocities of magnitude
2 v0 = Tinit " v 0 = 3Tinit 3
double InitTemp:
Initial temperature
!
14
For each atom, the velocity vector is then given by r r vi = v0 (! 0 ,!1 ,! 2 ) = v 0! r where ! is a rondomly oriented vector of unit length. RANDOM POINT IN A UNIT CIRCLE The question is how to randomly generate a point on a unit sphere. We first start from a simpler problem, random point in a unit circle.
y 1
-1
1
x
-1
(Algorithm) Let rand() be a uniform random-number generator in the range [0, 1] 1. i = 2*rand() - 1 (i = 0, 1) so that -1 i < 1 2 2. s 2 = ! 0 + !12 v 3. if s2 < 1, accept ! = (! 0 ,!1 ) 4. else reject and go to 1 (Probability distribution)
y 1 d!0 !1 -1 !0 d!1 1 x
-1
Probability to find a point in the shaded area = (# of points in the shaded area)/(total # of generated points) = P(0, 1)d 0d 1 = (d 0d 1 ~ shaded area)/(12 ~ area of unit circle), if there is no bias
! P(" 0 ,"1 ) = 1/ #
(Polar coordinate system)
(! 0 ,!1 ) = (s cos", ssin" ) 0 # s < 1;0 # " < 2$
15
y 1 d! ! -1 s ds 1 x
-1
Again the probability density is given by the ratio of the shaded area divided by the total circular area. ds sd! (if uniform) P(s,! )dsd! = " 12 s (1) ! P(s," ) = # RANDOM POINT ON A UNIT SPHERE The above distribution is closely related to the uniform distribution on a unit sphere. (Spherical coordinate)
z d! sin(!)d" ! " y
x
(!0 ,!1 ,!2 ) = (sin" cos#, sin" sin#,cos" ) 0 $ " < % ,0 $ # < 2%
(Probability distribution)
P(!,")d!d" =
d! sin!d" 4# 12 sin" ! P(",#) = 4$
(2)
" If we identify s ! sin , Eqs. (1) and (2) are equivalent. 2 "# " Q0 ! " < # $ 0 ! < $ 0 ! s % sin < 1 22 2 #$1 #' # ds # # sin d#d! sin & cos ) d#d! 2sin cos d#d! sin#d#d! sdsd! 2%2 2( 2 d# 2 2 // = = = = " " " 4" 4"
16
Note that
' # # 2 &0 2 )" 0 = sin# cos$ = 2sin 2 cos 2 cos$ = 2s 1% s s = 2 1% s & 0 ) # # & ( "1 = sin# sin$ = 2sin cos sin$ = 2s 1% s2 1 = 2 1% s2 &1 2 2 s ) 2# 2 ) " 2 = cos# = 1% 2sin = 1% 2s * 2
(AlgorithmGenerating random points on a unit sphere) Let rand() be a uniform random-number generator in the range [0, 1] 1. i =! 2*rand() - 1 (i = 0, 1) so that -1 i < 1 2 2. s 2 = ! 0 + !12 3. if s2 < 1, accept (" 0 ,"1," 2 ) = (2 1# s2 $ 0 ,2 1# s2 $1,1# 2s2 ) 4. else reject and go to 1 LINEAR CONGRUENTIAL RANDOM NUMBER GENERATOR Sequence ! seemingly random integers between 1 and m-1 is obtained by starting from some of nonzero I1 and
I j +1 = aI j (mod m)
We choose m = 231-1 = 2147483647 a = 16807 Uniform random number in the range [0, 1] is obtained as
r j= Ij / m .
17
Molecular Dynamics III: O(N) Linked-List Cell Algorithm
This chapter explains the linked-list cell MD algorithm, the computational time of which scales as O(N) for N atoms. We will replace function computeAccel() in the O(N2) program md.c by this algorithm. Recall that the naive double-loop implementation to compute pair interactions scales as O(N2). In fact, with a finite cut-off length, rc (see #define RCUT 2.5 in md.c), an atom interacts with only a limited number of atoms ~ (4/3)rc3(N/V), where V is the volume of the simulation box. The linked-list cell algorithm explained below computes the entire interaction with O(N) operations. CELLS First divide the simulation box into small cells of equal size. The edge lengths of each cell, (rcx, rcy, rcz) (double rc[3]), must be at least rc; we use rc = L/Lc, where Lc = L/rc ( = x, y, z) and L is the simulation box length in the direction (double Region[3]). Here x is the floor function, i.e., the largest integer that is less than or equal to x. An atom in a cell interacts with only the other atoms in the same cell and its 26 neighbor cells. The number of cells to accommodate all these atoms is LcxLcyLcz, r where Lc = L/rc ( = x, y, z) (int lc[3]). We identify a cell with a vector cell index, c = (cx, cy, cz) (0 cx Lcx1; 0 cy Lcy1; 0 cz Lcz1), and a serial cell index (see the figure below), c = cxLcyLcz + cyLcz + cz or cx = c/(LcyLcz) cy = (c/Lcz) mod Lcy cz = c mod Lcz.
r An atom with coordinate r belongs to a cell with the vector cell index,
c = r/rc ( = x, y, z).
!
!
18
LISTS The atoms belonging to a cell is organized using linked lists. The following describes the data structures and algorithms to construct the linked lists and compute the interatomic interaction using them.
DATA STRUCTURES An array implementation of the linked lists. lscl[i] holds the atom index to which the ith atom points. head[NCLMAX]: head[c] holds the index of the first atom in the c-th cell, or head[c] = EMPTY (= 1) if there is no atom in the cell.
lscl[NMAX]:
ALGORITHM 1: LIST CONSTRUCTOR
/* Reset the headers, head */ for (c=0; c<lcxyz; c++) head[c] = EMPTY; /* Scan atoms to construct headers, head, & linked lists, lscl */ for (i=0; i<nAtom; i++) { /* Vector cell index to which this atom belongs */ for (a=0; a<3; a++) mc[a] = r[i][a]/rc[a]; /* Translate the vector cell index, mc, to a scalar cell index */ c = mc[0]*lcyz+mc[1]*lc[2]+mc[2]; /* Link to the previous occupant (or EMPTY if you're the 1st) */ lscl[i] = head[c]; /* The last one goes to the header */ head[c] = i; }
where lcyz = lc[1]*lc[2]; lcxyz = lcyz*lc[0].
19
ALGORITHM 2: INTERACTION COMPUTATION
/* Scan inner cells */ for (mc[0]=0; mc[0]<lc[0]; (mc[0])++) for (mc[1]=0; mc[1]<lc[1]; (mc[1])++) for (mc[2]=0; mc[2]<lc[2]; (mc[2])++) { /* Calculate a scalar cell index */ c = mc[0]*lcyz+mc[1]*lc[2]+mc[2]; /* Scan the neighbor cells (including itself) of cell c */ for (mc1[0]=mc[0]-1; mc1[0]<=mc[0]+1; (mc1[0])++) for (mc1[1]=mc[1]-1; mc1[1]<=mc[1]+1; (mc1[1])++) for (mc1[2]=mc[2]-1; mc1[2]<=mc[2]+1; (mc1[2])++) { /* Periodic boundary condition by shifting coordinates */ for (a=0; a<3; a++) { if (mc1[a] < 0) rshift[a] = -Region[a]; else if (mc1[a]>=lc[a]) rshift[a] = Region[a]; else rshift[a] = 0.0; } /* Calculate the scalar cell index of the neighbor cell */ c1 = ((mc1[0]+lc[0])%lc[0])*lcyz +((mc1[1]+lc[1])%lc[1])*lc[2] +((mc1[2]+lc[2])%lc[2]); /* Scan atom i in cell c */ i = head[c]; while (i != EMPTY) { /* Scan atom j in cell c1 */ j = head[c1]; while (j != EMPTY) { if (i < j) { /* Avoid double counting of pair (i, j) */ rij = ri-(rj+rshift); /* Image-corrected relative pair position */ if (rij < rc2) Compute forces on pair (I, j) ... } j = lscl[j]; } i = lscl[i]; } } }
A neighbor cell can be outside of the central simulation box, i.e., the cell-index vector component c can be 1 or L ( = x, y, z). To find the atoms that belong to such an outside cell, the modulo operator (%) pulls back the cell index into the range [0, L1] ( = x, y, z):
c1 = ((mc1[0]+lc[0])%lc[0])*lcyz + ((mc1[1]+lc[1])%lc[1])*lc[2] + ((mc1[2]+lc[2])%lc[2]);
The atomic positions in an outside cell, on the other hand, must be treated as they are and must not r be pulled back into the central simulation box. The shift vector rshift (double rshift[3]) to shift the central image of such an atom to the appropriate image position by a simulati...