1 Molecular Dynamics Basics Basics of the molecular-dynamics (MD) method1-3 are described, along with corresponding data structures in program, md.c. Newton’s Second Law of Motion TRAJECTORY, COORDINATE, AND ACCELERATION • Physical system = a set of atomic coordinates: {!ri= (xi, yi, zi) | xi, yi, zi! ",i = 0,..., N #1}, where ! is the set of real numbers (in the program, represented by a double precision variable) and we use a vector notation, . (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. xyz12N • Trajectory: A mapping from time to a point in the 3-dimensional space, t ! " !"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 curve is then specified by a 3N-element vector, !rN= (x0, y0, z0, x1, y1, z1,..., xN!1, yN!1, zN!1). ri(t=0)vi(t=0)ri(t=t1)vi(t=t1) • Velocity: Short-time limit of an average speed (how fast and in which direction the particle is moving), !vi(t) =!"ri(t) =d!rdt! lim"#0!ri(t + ")$!ri(t)". double rv[NMAX][3]: rv[i][0], rv[i][1], and rv[i][2] are the x, y, and z components of the velocity vector, !vi, of the i-th atom. • Acceleration: Rate at which a velocity changes (whether the particle is accelerating or decelerating),2 !ai(t) =!""ri(t) =d2!rdt2=d!vidt! lim"#0!vi(t + ")$!vi(t)". double ra[NMAX][3]: ra[i][0], ra[i][1], and ra[i][2] are the x, y, and z components of the acceleration vector, !ai, of the i-th atom. vi(t)vi(t+!)vi(t)vi(t+!)ai(t)! The acceleration of a particle can be estimated from three consecutive positions on its trajectory separated by small time increments: !ai(t) = lim!"0!vi(t + ! / 2 ) #!vi(t # ! / 2)!= lim!"0!ri(t + !)#!ri(t)!#!ri(t) #!ri(t # !)!!= lim!"0!ri(t + !)# 2!ri(t) +!ri(t # !)!2 NEWTON’S EQUATION Newton’s equation states that the acceleration of a particle is proportional to the force acting on the particle, m!""ri(t) =!Fi(t), where m is the mass of the particle. For a heavier particle, the same force causes less acceleration (or deceleration). • Initial value problem: Given initial particle positions and velocities, !ri(0),!vi(0)( )i =1,..., N{ }, obtain those at later times, !ri(t),!vi(t)( )i =1,..., N;t ! 0{ }. Note that both positions and velocities must be specified in order to predict future trajectories. Potential Energy We consider forces which are derived from a potential energy, V (!rN), which is a function of all atomic positions. !Fk= !!!!rkV!rN( )= !!V!xk,!V!yk,!V!zk"#$%&'. Here partial derivative is defined as !V!xk= limh!0V (x0, y0, z0,..., xk+ h, yk, zk,...xN"1, yN"1, zN"1) " V (x0, y0, z0,..., xk, yk, zk,...xN"1, yN"1, zN"1)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.3 LENNARD-JONES POTENTIAL To model certain materials such as neon and argon liquids, the Lennard-Jones potential is often used by scientists: V!rN( )= u(rij)i< j!=i=0N"2!u(|!rij|)j=i+1N"1! where !rij=!ri!!rj is the relative position vector between atoms i and j, and u(r) = 4!"r!"#$%&12'"r!"#$%&6()**+,--. • 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 = gram•cm2/(second2) is a unit of energy) > σ = 3.4 × 10−8 cm * Don’t program with these numbers (in the CGS unit); it will cause floating-point under- or over-flows. NORMALIZATION • We introduce normalized coordinates, !!ri, potential energy, V′, and time, t′, that are dimensionless and are defined as follows. !ri=!!ri!= 3.4 "10#8[cm]"!!riV =!V!=1.66 "10#14[erg]"!Vt ="m /!!t = 2.2 "10#12[sec]"!t$%&'& With these substitutions, Newton’s equation becomes m!m"2"d2!!rid!t2= "!"#!V#!!ri. * In the derivation above, note that4 md2!ridt2= m lim!"0!ri(t + !)# 2!ri(t) +!ri(t # !)!2= m lim!"0!!$ri(t + !)# 2!$ri(t) +!$ri(t # !)[ ]!m /"%&'(2$!( )2 In summary, the normalized Newton’s equation for atoms interacting with the Lennard-Jones potential is given below. (From now on, we always use normalized variables, and omit primes.) d2!ridt2= !!V!!ri=!aiV (!rN) = u(rij)i< j"u(r) = 41r12!1r6#$%&'( 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. dudrr=r0= 0 ! r0= 21/6" 1.12u(r0) = 41212/6#126/6$%&'()= #1 Figure: Lennard-Jones potential. ANALYTIC FORMULA FOR FORCES !ak= !!!!rku(rij)i< j"= !!rij!!rki< j"dudrij Let’s evaluate the two factors in the above expression: 1. !rij!!rk=!!xk,!!yk,!!zk!"#$%&xi' xj( )2+ yi' yj( )2+ zi' zj( )2=2(xi' xj), 2(yi' yj), 2(zi' zj)( )2 xi' xj( )2+ yi' yj( )2+ zi' zj( )2"ik'"jk( )=!rijrij"ik'"jk( )!ddxf (x)[ ]1/2=12f (x)[ ]'1/2dfdx!"#$%& 2. dudr= 4 !12r13+6r7"#$%&'= !48r1r12!12r6"#$%&'!ddrr!n= !nr!n!1"#$%&'5 In the first result, !ik=1 (i = k)0 (i ! k){ is called Kronecker’s delta expression. Substituting these two results back into the expression for !ak, we obtain !ak=!rij!1rdudr"#$%&'r=riji< j(!ik!!jk( ) where !1rdudr=48r21r12!12r6"#$%&' Summary of Molecular-Dynamics Equation System Given initial atomic positions and velocities, !ri(0),!vi(0)( )i =1,..., N{ }, obtain those at later times, !ri(t),!vi(t)( )i =1,..., N;t ! 0{ }, by integrating the following differential equation: !""rk(t) =!ak(t) =!rij(t) !1rdudr"#$%&'r=rij(t )i< j(!ik!!jk( ) where !1rdudr=48r21r12!12r6"#$%&' The sum over i and j to evaluate the acceleration is implemented in a program as