Introduction To Molecular Dynamics Simulations. Prabal Maiti
Introduction To Molecular Dynamics Simulations. Prabal Maiti
Introduction To Molecular Dynamics Simulations. Prabal Maiti
6.1 Introduction
Laplace vision . . . An intelligence which could, at any moment, comprehend all the forces by which nature is animated and respective positions of the beings of which it is composed, and moreover, if this intelligence were far-reaching enough to subject these data to analysis, it would encompass in that formula both the movements of the largest bodies in the universe and those of the lightest atom: to it nothing would be uncertain, and the future, as well as the past would be present to its eyes. The human mind oers us, faint sketch of this intelligence [1]. Molecular dynamics (MD) can be termed as Laplaces vision of Newtonian mechanics on supercomputers. Molecular dynamics is the term used to describe the solution of the classical equation of motion (Newtons equations) for a set of molecules. In MD particles move at constant velocity between perfectly elastic collisions, and it is possible to solve the dynamic problem without making any approximations, within the limits imposed by machine accuracy. Computer simulation has been very powerful tool to attack many body problems in Statistical physics, Physical chemistry and Biophysics. Although the theoretical description of complex system in the framework of statistical physics is rather well developed and experimental techniques for detailed microscopic are sophisticated, it is often possible to study specic aspect of those system in great details via simulation only. There are only handfuls of non-trivial, exactly
174
soluble problem in statistical mechanics, 2-D Ising model and ideal gas being example. Some problems in statistical mechanics, while not being exactly soluble, are analyzed based on straightforward approximation scheme. Problem is that for many problems of interest these straightforward approximations do not work. For many systems it may not be even clear that how to begin constructing an approximate theory. In all such cases computer simulation can help a lot by providing sucient input about microscopic properties of complex systems of interest. Computer simulations provide exact results for the problems which otherwise only be soluble by approximate methods, or might be quite intractable. Thus computer simulations can be used to test approximate theories. Results of computer simulations can also be compared with real experiment than it can be used to test underlying model. Simulation oers insights to experimentalist to interpret new result on the basis of microscopic details of system. Computer simulations act as a bridge between microscopic world and macroscopic world of laboratory. We provide microscopic input about system (masses of constituent atoms and interaction between them), we get macroscopic properties measurable in experiments, like equation of state, diusion constant, correlation functions to name a few as simulation output. Simulation are also used to study properties under extremes of temperature and pressure, which may be dicult or impossible to carry out experimentally, however can provide details useful for lot of technological applications as well as for academic interest. Simulation has also been used to improve our understanding of phase transition and behavior at interfaces. Before we go into the details of computer simulation technique we give a brief history of the development of this eld.
175
well as many other condensed matter systems. Computer simulation has also advanced to non-equilibrium system, stochastic dynamics and incorporation of quantum eects. This article is by no means a exhaustive treatment of the subject. There exist several excellent text books[11, 13] which give a very comprehensive understanding of the various advanced concepts of the simulation methodology. The reader is referred to these for further details. The rest of the article is organized as follows: In section 2 we describe the model building, in section3 we describe the basic of Hamiltonian dynamics as a basis for the equation of motion and will discuss various schemes for integration used in simulations. Section4 is about various complexities arising during simulation and way of dealing with those complexities. In Section5 simulation techniques for dierent ensembles has been discussed. Finally we conclude this article by giving two case studied with MD simulations.
176
(a)
(b)
(c)
(d)
Figure 6.1: Examples of molecular models at various levels. (a) Atomistic model of a class of banana liquid crystal (b) coarse-grained representation of the banana molecules using spherocylinder of aspect ratio L/D (c) atomistic model of a two-tail surfactant (d) coarse-grained representation of various class of surfactants using simple bead-spring model. Moreover for many problems of interest these approximate simple models seem to be adequate enough. With the power of modern computer and recent algorithm development it is possible to take care of all quantum eects. To describe molecular charge, a set of ctitious charges are distributed in such a way that this distribution reproduce known multipole moments. Microscopic state of a system may be specied in terms of the positions and momenta (assuming classical description is adequate) of constituent set of particles. Within adiabatic approximation, we neglect electronic motion and consider only nucleic motion. This approximation is excellent for a wide range of system but is unsuitable for reaction involving electronic rearrangement such as bond formation and cleavage, polarization and chemical bonding of metal ions. Quantum approaches are used for such studies. For a classical system, specifying the instantaneous positions and momenta of all the particles constituting the system can specify the microstate at any time t. For N particles there are 3N coordinates q1 , q2 , , q3N and 3N conjugate momenta p1 , p2 , , p3N . The Hamiltonian of the system can be written in term of these qi s and pi s as follows H(q, p) = K(p) + V (q) q = (q1 , q2 , ...........qN ) p = (p1 , p2 , ..........pN ) (6.1)
K and V being kinetic and potential energy respectively. q is the generalized coordinate, it may be Cartesian coordinate of each nucleus or Cartesian coordinate of each center of mass (COM) with orientation parameter in case of molecules and rigid bodies p is conjugate momenta. Kinetic energy is dened as
177
K=
i=1
p2 /2mi i
(6.2)
where i stands for particle number and represents various components of momenta. Potential energy V depends on intermolecular interaction. We will describe various potential for dierent kind of systems.
V =
i
v1 (ri ) +
1 2 i i>j
v2 (ri, rj ) +
i j>i k>j>i
where v1 represents eect of external eld on individual atoms, v2 is pair potential, and is most important in computer simulation. In general we consider only pair potential. Once we have information of Hamiltonian, we have equations of motion given by
d2 r = V (r) dt2
(6.4)
This governs the evolution of the system and is at the heart of Molecular dynamics simulation. Numerically there are various schemes to solve the above dierential equation. In the next section we will discuss few of them. A major challenge in the eld of molecular dynamics is to describe the inter-atomic potential V (r) (inter and intra-molecular potential in case of molecular systems). The accuracy and validity of the simulation results critically depends on the quality of the potential V (r). There exist several well-known empirical force-elds which give forms as well as parameters for a variety of intra and inter-molecular potential. Below we discuss those.
Where Etotal , EvdW , Eelec , Ebond , Eangle , and Etorsion are the total energy, the van der Waals, electrostatic, bond stretching, angle bending, and torsion energy components, respectively.
178
The van der Waals interaction is given by the 12 6 LJ interaction of the following form EvdW (R) = D0 R0 R
12
R0 R
(6.6)
where D0 is the strength of the interaction and R0 is the range of interaction. These constants can be found either by rst principle quantum calculation or by experimental input. In general these parameters are available for various atom types from various well known force elds such as Dreiding, AMBER and CHARMM. Interaction parameters between dierent kind of atom types (like between C & O) are calculated by taking either arithmetical mean or geometrical mean. Other commonly used form of non-bond interaction is the so-called Buckingham potential and has the following form EBuckingham = nonbonded pairs A exp(crik ) B 6 rik (6.7)
Again constants appearing in these potential are given in available force elds. In case one needs to maintain special coordination using non-bond interaction Morse potential is the choice and has the following form EMorse = nonbonded pairs where Do , Ro and are constants. Electrostatic interaction between charged particles can be computed using Coulombs law and given by ECoulomb = nonbonded pairs qi qk rik (6.9) Do 1 e(RRo )
2
(6.8)
qi s being charge on each atom and being dielectric constant. In fact computationally this is the most expensive calculation because of its long range character. In Complexity of force calculation section, we will discuss various algorithms to compute Coulombs interaction eciently. Apart from these various non-bond interactions discussed above for molecular systems we have also the following bonded interaction (as shown in Fig.6.2) present in the systems to maintain molecular topology.
179
rij
ijk
Figure 6.2: Dierent intra-molecular interactions. The intra-molecular interaction corresponding to the bond stretching potential can be described by the harmonic potential of the following form 1 2 Kb (rij r0 ) (6.10) 2 The harmonic potential is adequate only for small deviations from reference value. For large deviations Morse potential is the choice. This is particularly useful for maintaining specic coordination in a given geometry. This has functional form: Ebond (rij ) = EMorse =
1,2pairs
Db {1 exp[Km (r r0 )]}
(6.11)
Where the parameters Db and Km characterize the well depth and well width respectively. Similarly angle bending potential can be described by the following harmonic form 1 K ( ijk 0 )2 (6.12) 2 Many available force elds however use a cosine harmonic form rather than simple theta harmonic as given above. In case of cosine harmonic theta the angle bending potential has the following form Eangle (ijk ) = Eangle (ijk ) = 1 K (cos ijk cos 0 )2 2 (6.13)
Apart from the bond stretching and angle bending potential to maintain certain topology of the molecular system torsional potentials are included. Mainly two types of torsional potentials are used: dihedral angle potential and improper torsions. The dihedral angle potential is mostly used to constrain the rotation around a bond and involves four consecutive bonded atoms (i-j-k-l) as shown in Fig.6.2. Improper torsion is used to maintain planarity of certain atoms and involves four atom which are not bonded in the same sequence as i-jk-l. Out-ofplane bending is incorporated through the improper torsion potential. In many force elds instead of improper torsion an inversion potential is used. One of the most commonly used dihedral potential is the cosine form given by
180
E(ijkl ) =
n
Vn [1 + cos(nijkl 0 )] 2
(6.14)
Where n is an integer. For each such rotational sequence described by torsion angle , n denotes the periodicity of the rotational barrier, and Vn is the associated barrier height. 0 is reference torsion angle. Commonly used improper torsion is harmonic in nature and is described by following functional form
I E(ijkl ) = kijkl (ijkl 0 ) I Where kijkl is the force constant. 2
(6.15)
181
6.3.1 Verlet AlgorithmTaylor Expansion for position around time t is r(t + t) = r(t) + v(t)t + 1/2a(t)t2 + b(t)t3 + Ot4 If we reverse the time we have r(t t) = r(t) v(t)t + 1/2a(t)t2 b(t)t3 + Ot4
2
(6.16)
(6.17)
Adding the above two equations and keeping terms up to t we getr(t + t) = 2r(t) r(t t) + a(t)t2 + Ot4 v(t) = r(t + t) r(t t)/2t (6.18)
This is the update equation for position. Similarly for velocity update we have (6.19)
Equation (6.18) and (6.19) are the update equations in Verlet algorithm[8]. Thus Verlet algorithm takes position and acceleration at time t and position at previous time step as input. Note that for Verlet integrator we need position at previous time step which can be obtained by simple approximation as follows r(t0 t) = r(t0 ) v(t0 )t (6.20)
It is easy to see that Verlet integrator indeed satises the above listed criteria. Time Reversible Forward time step equation r(t + t) = 2r(t) r(t t) +
2 1 m F (t)t
(6.21)
2 1 m F (t)(t)
(6.22) (6.23)
2 1 m F (t)t
Thus we have same algorithm to move system backward with same force and same position. Accuracy Verlet algorithm calculates position accurate up to order of t4 and velocities are accurate up to t2 . However it is not very accurate in terms of numerical precision as it adds large and small numbers in same equation. r(t + t) r(t) = r(t) r(t t) +
2 1 m F (t)t
(6.24)
Here term on left hand side is of order twhile on right hand side rst term is of the order of t0 second term of t1 and third term is of t2 .
182
Memory usage and eciency It requires 9N (3N position of previous step and 6N for current step) words of memory so it is compact. Verlet algorithm is simple and easy to implement in a MD code. Area preserving property This algorithm preserves area in phase which is essential criteria for an integrator[14]. It has excellent energy conservation even at large time step.
6.3.2 Leap-Frog AlgorithmLeap- Frog algorithm can be obtained by simple algebraic manipulation of the Verlet integrator. It eliminates one of the major disadvantage of the Verlet algorithm namely the addition of small numbers O(t2 )to dierences in large ones O(t0 ). In this scheme rst velocity at half time-step is calculated which in turn is used to update the position at full time step. The equations for Leap-Frog algorithm is as follows r(t + t) = r(t) + v(t + 1 t)t 2
1 v(t + 2 t) = v(t 1 t) + 2 1 m F (t)t
(6.25) (6.26)
The Leap-Frog algorithm is equivalent to the Verlet algorithm as can be seen easily as follows. Substituting equation (6.26) in equation (6.25) we have
1 r(t + t) = r(t) + v(t 2 t) + 1 m F (t)t
(6.27)
From equation (6.25), we can also get r(t) evaluated at previous time step as follows
1 r(t) = r(t t) + v(t 2 t)t
Substituting the value of v(t 1 t) from above equation in into equation (6.26), 2 we get r(t + t) = r(t) + (r(t) r(t t)) + . = 2r(t) r(t t) +
2 1 m F (t)t
2 1 m F (t)t
(6.28)
Equation (6.28) is equivalent to the original Verlet equation as given by equation (6.18). Note that for Leap-Frog algorithm we need to have velocity at the previous time step which can be obtained by simple approximation as follows v(t0 t) = v(t0 )
1 1 m F (t0 ) 2 t
(6.29)
Velocity at current time step can be obtained using simple interpolation like
183
(6.30)
One of main advantage of the Leap-Frog integration scheme is that we need not add numbers which are dierent order in t.
(6.31) (6.32)
[F (t) + F (t + t)] t
Basic Verlet Scheme can be recovered from these equations by eliminating velocity. This algorithm also requires storage of 9N words i.e. position, velocity and acceleration at current time step. Although it is not implemented same as above, it involves one intermediate step. First position is updated with equation (6.18) and velocities at mid step are computed using
1 v(t + 2 t) = v(t) + 1 2m
[F (t)] t
(6.33)
The force and acceleration at time t+t are then computed and velocity move is computed with
1 v(t + t) = v(t + 2 t) + 1 2m
[F (t + t)] t
(6.34)
Due to stability, accuracy and simplicity this is the most preferred choice of integrator.
184
2. Force evaluation: Use the predicted position to compute the force and acceleration at the predicted positions. The resulting acceleration will be in general dierent from the predicted acceleration in previous step. 3. Corrector: use the new acceleration to correct the predicted position, velocities and acceleration. Use the Taylor expansion of the position at time t+dt t2 t3 t4 a(t) + b(t) + c(t) + ... 2 6 24
(6.35)
Where v is the velocity, a is the acceleration, b is the third derivative of position, c is the fourth derivative etc. Using the Taylor expansion for velocity, acceleration etc. we have v(t + t) = v(t) + a(t)t + t3 t2 b(t) + c(t) + ... 2 6 t2 c(t) + ... 2 (6.36)
(6.37) (6.38)
The dierence between the predicted (step 1) and calculated (step 2) acceleration is given by a(t + t) = ac (t + t) a(t + t) (6.39)
and is used to correct the positions and velocities in the correction step as follows rc (t + t) = r(t + t) + c0 a(t + t) v c (t + t) = v(t + t) + c1 a(t + t) ac (t + t) = a(t + t) + c2 a(t + t) bc (t + t) = b(t + t) + c3 a(t + t) (6.40) (6.41) (6.42) (6.43)
The values of the coecients depend on the order of the Taylor series expansion. Gear[15] has suggested the best values of the set of coecients c0 , c1 , c2 , c3 .
185
Let us consider a simple one-particle system in one dimension with a Hamiltonian. H= The equations of motion are q= p m p= dU = F (x) dx (6.46) p2 + U (x) 2m (6.45)
Now Liouvilles theorem states that any phase space function A(x, p) evolves according to the following equation dA = {A, H} dt Where {A, H} is the Poisson bracket and is given by {A, H} = H A H A p x x p (6.48) (6.47)
The evolution equation (6.47) gives back Hamiltons equation of motion: To see this consider A(x, p) = x. Then dx = x = {x, H} dt {x, H} = p x dU x p x = since =0 m x dx p m p (6.49)
So we have x = p/m x = p/m. Similarly if we consider A(x, p) = p, we have the equation of motion for p as
186
(6.50)
Thus we get same equations of motion as in equation (6.46). Now dene a two-dimensional phase space vector = (x, p). Hamiltons equation of motion for is given by d = {, H} (6.51) dt We dene Liouville operator L such that iL ={,H} The equation of motion given by equation (6.51) in Liouville operator form is given by, d = iL dt Solution of this equation of motion is given by (t) = eiLt (0) (6.52)
(6.53)
The operator exp (iLt ) is called the classical propagator and acts as phase space evolution operator. Notice that the presence of i allows one to make an analogy with the QM propagator exp (-iHt/h). In general it is dicult to evaluate exp (iLt) because of the following reason. iL can be written as iL = p + F (x) = iL1 + iL2 m x p (6.54)
p Where, iL1 = m x and iL2 = F (x) p The diculty in any computation arises from the fact that iL1 and iL2 do not commute: [iL1 , iL2 ] = 0. Since they dont commute,
(6.55)
Trotter Theorem helps us to evaluate evolution operator in the following way: e(iL1 +iL2 )t = lim eiL2 t/2M eiL1 t/M eiL2 t/2M
M
(6.56)
187
e(iL1 +iL2 )t = eiL2 t/2M eiL1 t/M eiL2 t/2M This implies:
(6.57)
(6.58)
The expression on the left looks like approximate propagation of the system up to time t by M application of the operator in the bracket. If we interpret t/M as single time step, t, then we have eiLt = e(iL1 +iL2 )t eiL2 t/2 eiL1 t eiL2 t/2 (6.59)
This is the propagator U(t) for time step t. U(t) is unitary and preserve the time reversibility of the dynamics. U + (t)U (t) = I (6.60)
The unitarity of the propagator implies time reversal symmetry in the equations of motion. If the system is propagated forward in time up to a time t and then the clock is allowed to run backwards for a time -t, the system will evolve according to the same equations of motion but the direction of the velocities will be reversed, so that the system will simply return to its initial condition.To see this we note that U (t) = exp (iLt) Now apply U (t) on (0) to get (t) followed by U (t) : (t) = U (t) (0) Operating U (t) will give back the U (t) U (t) (0) = eiLt eiLt (0) = (0) So we have U (t) U (t) = I (6.61) Equation (6.61) implies time reversibility. Another important property of the unitary operator U(t) is that its determinant is one. This is consistent with the fact that volume in phase space remains conserved under Hamiltons equation. So we have demonstrated that evolution propagator U (t) satises two fundamental criteria of the Newtons equation of motion namely the time reversibility and the area preserving property. So equation of motion derived from such evolution operator will have these properties in-built. Now we will demonstrate that evolution operator gives back the original Verlet equation as given in equation(6.19)-(6.20). To see this we rst apply equation (6.59) on x, the position of particle at time t to get position x(t+dt)
188
U (t)x = e 2 F (x) p et m x e 2 F (x) p x t p = e 2 F (x) p et m x x t t = e 2 F (x) p (x + m p) t t = x + m (p + 2 F (x) So we have the update equation for the position given by t2 t p(t) + F (x, t) m 2m
x(t + t) = x(t) +
(6.62)
Similarly when we apply equation (6.59) on p, the velocity of particle at time t, we get velocity of the particle at time t+dt to as follows U (t)p = e 2 F (x) p et m x e 2 F (x) p p t p = e 2 F (x) p et m x (p + t F (x)) 2 t = e 2 F (x) p (p + t F (x + t p)) 2 m t = p + t F (x) + t F (x) + m (p + 2 2 So we have the velocity update equation as follows
t p t
t 2 F (x))
t 2 [F (x, t) t 2 [F (x(t))
+ F (x +
t mp
t2 2m F (x))]
(6.63)
+ F (x(t + t))]
Equations (6.62) and (6.63) are same as the Verlet equation as given by equation (6.18) and (6.19). So Verlet equations have both the time reversible and area preserving properties.
6.3. Various Schemes of Integration Vibrational mode Wave number (1/) cm1 3200-3600 3000 2100 1700 1600 700 PeriodTp (/c) fs 9.8 11.1 15.9 19.6 20.8 47.6 Tp /2 (fs) 3.1 3.5 5.1 6.2 6.4 15
189
O-H, N-H stretch C-H stretch CC, CN stretch C=C stretch H-O-H bend O-C-O bend
Table I: Typical time scale associated with various vibrational modes in molecular system. The table is based on the values given in reference 13. So for simulation of systems having hydrogen bond we can use 1 fs time step without having any problem with energy conservation. Another solution to this problem is to restrain these fast degrees of freedoms while solving the un-constrained degrees of freedom. Bonds involving H have highest frequency hence they are constrained during dynamics so that larger time step can be used. Several algorithms exist for this purpose. The SHAKE algorithm for bond constraints was introduced by Ryckaert et al [16]and is widely used in molecular simulation. A full description of SHAKE algorithm is outside the scope of this article. The reader is referred to the original article for further details. Basic idea of SHAKE is to use Lagrange multiplier formalism to enforce bonds distances constant. Suppose we have Nc such constrained given by
2 2 k = rk1 k2 Rk1 k2 = 0,
where
k = 1, 2, 3 . . . . . . . . . ..Nc
Rk1k2 being constrained distant between atoms k1 and k2 atoms. This leads to modied constrained equation of motion mi d2 ri (t) = [V (r1 ....rN ) + dt2 ri
Nc
(6.64)
Where mi is mass of ith particle and k is the Lagrange multiplier(unknown) for k th constraint. This equation can be solved for unknown multiplier by solving Nc quadratic coupled equations. And we get the following equation of motion
uc rk1 (t + t) = rk1 (t + t) 2(t)2 m1 k (t)rk1 k2 (t) k1 uc rk2 (t + t) = rk2 (t + t) 2(t)2 m1 k (t)rk1 k2 (t) k2
(6.65)
Where ruc is position updates with unconstrained force only. This procedure is repeated till dened tolerance given as |rk1 k2 (t + t) Rk1 k2 | Rk1 k2 (6.66)
190
Figure 6.3: In periodic boundary condition central simulation cell is replicated in all direction to form an innite lattice.
191
Figure 6.4: Illustration of minimum image convention. The red particle(A) interacts with particles within box (dark boundary). Though yellow(B) and green(C) particles are at same distance, red(A) particles interacts with green(C) only, as yellow(B) particle is out of box
box. There are N 1 such term. However, we must also include interaction coming from images lying in the surrounding boxes. That is an innite number of terms and it is not feasible to include interactions with all images. For short-range interaction this problem is resolved by invoking what is called minimum image convention. In this case we construct a simulation box of same size as the original box with the red atom at its center. Now minimum image convention says that the red atom interact with those atoms which lie in this region, that is with the closest periodic images of the other N 1 atoms (see Fig.6.4).
192
B A B B
Figure 6.5: Number of pair computation can be reduced by introducing a spherical cut-o. Here the red(A) particle only interacts with particles within cuto sphere i.e. pink(B) colored particles. 1. The cuto distance should be smaller than L/2, where L is the minimum of the box dimension in any periodic direction. 2. Thermodynamic properties are dierent for a truncated potential compared to non-truncated case. However, we can apply long range correction to get back approximately the non-truncated properties. 3. Cuto introduces discontinuity in the force and energy computation. This has serious consequence on the energy conservation and stability of the simulation
This additional term is constant and does not aect the force and hence equation of motion of the system. However, the force is still discontinuous at or near cut-o. This can cause numerical instability in the numerical solution of
193
the dierential equations. To avoid this one can use shifted-force potential[18] where a linear term is added to the potential, so that its derivative is zero at the cuto. The form of the potential with the shifted force is given as follows usf (r) = u(r) u(rc )
du dr r=rc
(r rc ) r rc r > rc
(6.68)
The discontinuity now appears in the gradient, not in the force itself. The force goes smoothly to zero at cuto. But this shift makes the potential deviate from the true potential, so the calculated thermodynamics properties will be changed. However, the thermodynamic properties of the system interacting with un-shifted potential can be recovered from shifted-force potential simulation result using perturbation scheme. But it is dicult to do so and is rarely used in simulation.
The switching function S(r) gradually tapers the potential between two cutos: it smoothly changes its value from 1 at a cut-o distance rl (lower cuto) to a value of 0 at cuto distance ru (upper cuto) and satises the following criteria Sr=rl = 1.0 Sr=ru = 0
dS dr r=rl dS dr r=ru
=0
d2 S dr 2 d S dr 2
2
r=rl r=ru
=0 =0 (6.70)
=0
Zero rst derivatives ensure that the force approaches to zero smoothly at the cutos. A continuous second derivative ensures the stability of the integration algorithm. The lowest order polynomial which satises the above criteria is the third order polynomial (or cubic spline) and given by
2 2 2 2 S(r) = (ru r2 )2 23 [ 3 (ru r2 )] where = ru rl 2 2 2 2 2 2 2 2 3 = c0 (ru r ) + c1 (ru r ) + c2 (ru r )
(6.71)
For better accuracy higher order polynomial can also be used for switching function. For example we can use 5th order polynomial of the following form S(r) = c0 + c1 +c4
rrl ru rl rrl ru rl
+ c2
4
rrl ru rl rrl ru rl
+ c3
5
rrl ru rl
(6.72)
+ c5
194
i rl ru
Figure 6.6: Scheme for constructing Verlet list for non-bond computation.
The coecients can be obtained using the boundary condition given in equation (6.70) and are given as follows c0 = 1, c1 = 0, c2 = 0, c3 = 10, c4 = 15, c5 = 6. We can use higher order polynomial also.
195
196
Figure 6.8: (a) the original charge distribution, (b) A Gaussian charge distribution of opposite sign is added at every charge site, (c) another Gaussian charge distribution of Scheme for Ewald summation. 2. Potential due to Gaussian screening charge cloud with charge qi 3. Potential due to compensating charge cloud with charge qi The contribution to the electrostatic potential at point ri due to a set of screened charges can easily be computed by direct summation because the electrostatic potential due to screened charge is a rapidly decaying function of r. A Gaussian charge distribution of following form is commonly used (r) = qi 3 exp(2 r2 ) 3/2 (6.73)
controls the width of the distribution. Using the solution of Poissons equation either in real space or in Fourier space potential for the above three cases can be computed. One thing should be noted here, in the above mentioned recipe interaction of cancelling distribution centered at a site with itself is also included. Hence it should be subtracted from the above mentioned sum. This is termed as the self correction. Another correction comes from smeared nature of charge. Thus total potential energy due to Long range Coulomb interaction is given by the expression
197
Uc = Uq () Uself () + U ()
(6.74)
The rst term Uq () represents the reciprocal sum due to compensating charge cloud with charge qi and is given by Uq = = =
1 2 1 2 1 2 i qi (ri ) qi qj ik(ri rj ) 4V k2 /4 i,j V 2 e k=0 k2 e 2 4V k /4 |(k)|2 k=0 k2 e
(6.75)
(6.76)
And nally U () is the potential in the real space due to the interaction between the original charges qi plus the Gaussian screening charge cloud with charge qi and is given by U = =
1 2 1 2 n n
(6.77)
exp(u2 ) du
(6.78)
And the complementary error function erf c is dened as follows erfc(x) = 1 erf(x) 2 = x exp(u2 ) du (6.79)
Larger the value of , sharper the distribution hence large number of K summation has to be included for better accuracy. On the other hand, large value of reduces range of screened potential hence we can use smaller cuto radius. Hence value of is optimized between these two factors to give better accuracy and eciency. Note that Ewald summation as presented above scales as N 2 only. However, with suitable choice of and k-space summation cut-o K, Finchman[21] was able to optimize the summation which scales as N 3/2 . The optimized values for and k-space summation cut-o K is given by = /R = K (6.80)
Where R = R ; = L; K = KL and R and L are the real space cut-o L 2 (same as rc as discussed above) and the simulation box length respectively. Ewald summation can further be optimized through the use of Fast Fourier transform (FFT) in evaluating the reciprocal summation. This goes in the name
198
of Particle Mesh Ewald (PME)[22]. Particle Mesh based approaches rely on the use of xed cuto on the direct space sum together with an FFT based approximation for the reciprocal space sum that scales as N log(N).
1 2
2 N 2 mi vi i=1 3 N kB
At each time step the velocities are multiply by and T(t) is calculated from the KE at time t. This drastic approach, however, implies rapid energy transfer to, from and among the various degrees of freedom in the system. In particular it can be shown that velocity rescaling leads to an articial pumping of energy into low
199
frequency modes. This does not represent any statistical ensemble. However this algorithm is simple and easy to implement.
(6.84)
1 2 mi vi (t)}/t] 2
(6.85)
N being the total number of particles Form equation (6.83) we have the the change vi of the velocity over a time interval of t = 0 t is give by (assuming t small) dierence in velocity in two time step given by vi = vi (t + t) vi (t) = 1 mi
t+t t
[Fi (t ) mi vi (t ) + Ri (t )]dt
(6.86)
dt
t+t t
dt Ri (t )Ri (t ) = 6N mkT0 t
(6.87)
vi Fi + 2
i=1
3N kT0 Ek 2
(6.88)
First term on the right hand side is equal to minus of the time derivative of potential energy and is related to the eect of sytematic force which will have no eect of the imposed thermostat. Hence the coupling to the heat bath is
200
represented by the second term and can be associated with the time dependence of the system temperature dT = 2 (T0 T ) dt (6.89)
where the time constant for heat bath coupling T is equal to 2 1 . So the temperature deviation decays exponentially with time with time constant T and the equation of motion can be written as mi vi = Fi mi T0 1 vi T (6.90)
In the Berendsen weak coupling scheme coupling to the heat bath and subsequent temperature controll is achieved by appropiate rescaling of the velocities by a time dependent scaling factor during the integration of equation of motion. In the following we show the scaling factor in the context of leap-frog integrator equ. (6.26). t 2 t 2 t t 2
ri t +
= =
(t)r i t + (t) ri
+ m1 Fi (t)t i
(6.91)
To acheive temperature variation consistent with equ.(6.89) the scaling factor can be found by imposing [42] t 2 t 2 t 2
t+
=T
1 + T t T0 T
(6.92)
t 2
(6.93)
T0 t+
1/2
t 2
(6.94)
So the temperature is controlled by scaling the velocities of the particle as each time step with a time dependent constant given by,
201
= 1+
t T
T0 1 T
1/2
(6.95)
If T is large, then the coupling will be weak. If T is small, the coupling will be strong and when the coupling parameter equal to integration time step (T = t) then this algorithm is equivalent to simple velocity rescaling method. A good value of T is 0.5 1ps when integration time step t = 1fs. Advantage 1. Strength of the coupling can be varied and adapted to the use requirement 2. Very easy to code 3. Very ecient to bring the system to a desired temperature. Disadvantage 1. Does not represent a true canonical ensemble. Velocity rescaling articially prolongs any temperature dierence among components of the system, which can lead to the phenomena of hot solvent and cold solute, even though the temperature of the system is at its desired value. This can be avoided by having separate temperature coupling to the solute and solvent, but this leads to the unequal distribution of energy among various components.
(6.96)
(6.97)
(6.98)
202
is the internal virial for pair additive potential. Fij is the force on particle i, due to particle j and V is the volume of the simulation cell. During the simulation pressure can be changed by changing the virial which can be accomplished by scaling inter-particle distances. Constant pressure simulation is achieved by scaling the coordinates along with appropriate volume scaling. At each time step particle coordinate x is scaled to x and box length l to lwhere is given by = [1 t (P P0 )]1/3 p (6.99)
For further details of the Berendsen weak coupling method the reader is referred to the original paper.
203
Partial positive atomic charge on hydrogen (+0.4170e) is balanced by negative charge on oxygen (0.8340e). Intermolecular force between two water molecules is described by Lennard-Jones potential with a single interaction site centered on oxygen atom with potential depth of = 0.6364 KJ/mol and range of = 3.1506 . In CNT C-C bond distance length r0 = 1.43and C-C-C A angel is xed to 0 = 2/3. Interaction of water molecules with carbon atom is described by Lennard-Jones potential having potential depth = 0.1143 Kcal/mol and range = 3.275. The solvated structures were subjected to A 1000 steps of steepest descent minimization of potential energy, followed by another 2000 steps of conjugate gradient minimization. During this minimization the nanotube structure was kept xed in their starting conformations using a harmonic constraint with a force constant of 500 kcal/mol/. This allowed A the reorganization of the water molecules to eliminate bad contacts with the nanotube. The minimized structure was then subjected to 45 ps of MD, with 2 fs time step. During the dynamics, the system was gradually heated from 0 to 300 K with harmonic constraints on the solute using the SHAKE method. This was followed by 200 ps constant volume constant temperature (NVT) dynamics with a temperature-coupling constant of 0.5 1.0ps on the solute. Finally, 20-30 ns NPT production dynamics was carried out with a time constant for heat bath coupling of 1ps. The electrostatics interactions were evaluated with the Particle Mesh Ewald [20] (PME) method, using a real space cut o of 9. Again during this production runs nanotube was held xed to its starting A conguration using a force constant of 500 kcal/mol/2 . The trajectory was A saved forevery 1 ps interval and was used for data analysis. ResultsDespite of hydrophobic nature of carbon nanotubes, water spontaneously goes inside CNT in quantitative agreement with earlier works[29]. Water enters from one end of the tube and leave from other end. Water molecules inside the tube are arranged in a single le with all dipole either pointing up or down as shown in Fig.6.9(b) & (c). Figure 6.10 shows the typical trajectories of four neighbouring, conned water molecules in the single le chain. The extreme correlation among the trajectories is apparent. To quantify the positional ordering of the water molecules inside the nanotube we have calculated the pair-correlation function for the water molecules inside the nanotube using g(z) = 1 N
N N
i=1 j=1,j=i
(z zij
(6.100)
where zij is the axial separation between the ith and the j th water molecules,N is the number of water molecules inside the nanotube and the angular brackets indicate an average over time. We have also looked at the reorientational dynamics of the conned water. The water molecules are tightly packed inside a nanotube due to hydrogen bonding, with the average density almost four times that of bulk water. The
204
(a)
(b)
(c)
Figure 6.9: (a) TIP3P model of water (b) Water inside carbon nanotube (c) typical ordering of water inside CNT
Figure 6.10: The z coordinates of four conned, neighbouring water molecules as a function of time. dipole moments of all the conned water molecules are almost always aligned, pointing either up (along +z direction, see Fig.6.9b & 6.9c) or down (along -z direction). Their orientation changes by cooperative ips that take them from one of these states to the other. Fig.6.11 shows the time-dependence of
205
Figure 6.11: (a) and (b) show the average dipole moment of the water molecules conned inside the 14and 28and nanotubes respectively. The 14and A A A 28tubes accommodate 5 and 10 water molecules on an average. The dipole A moments of the water molecules inside the nanotube are mostly aligned in either up or down states, with cooperative ips between these two states. The ips become rarer as the length of the water chain increases. The red curve in 11(b) shows the axial component of the average dipole moment of the water molecules in the bulk, where there is no orientational order. the axial component of the average dipole moment, Mz (t) of the conned water molecules, dened by 1 Mz (t) = N (t)
N (t)
i=1
pi (t) n
(6.101)
Where N (t) is the total number of water molecules inside the tube at time t and pi (t) is the dipole moment of ith water molecules inside tube at time t and n represents axis of the tube. Figures 6.11(a) and 6.11(b) show data for water conned inside 14and A 28long nanotubes, respectively. The average number of conned water A molecules is 5 and 10 inside the 14and 28nanotubes, respectively. The net A A dipole moment of the chain of water molecules makes collective ips between the up and down states. The mean time interval between successive ips increases with the length of the tube. The red curve in Fig.6.11(b) shows the time-dependence of the axial component of the average dipole moment of the water molecules in bulk water outside the nanotube. This clearly shows that connement leads to ordering of the dipole moments of the water molecules. Further details about these interesting reorientational dynamics can be found in our original works [26, 30].
206
207
Figure 6.12: Energy versus strain for BNNTs for single wall nanotube (SWNT) and double wall nanotube (DWNT) of dierent chiralities.
(a)
(b)
(c)
(d)
(e)
Figure 6.13: Snapshots of (7,7) single wall BNNT under compression, snapshots correspond to discontinuities in Fig.6.12. The elastic properties of the C/BN nanotubes arise due to the strength of in plane C-C/B-N bonds in comparison with the ease of out of plane deformation. The strength of the tube is reected in the high Young modulus. Young modulus is calculated from the second derivative of the strain energy, Y = 1 2 E() V0 2 (6.102)
=0
208
Figure 6.14: Youngs modulus Y as a function of tube radius and charge on B/N atoms in BNNT.
Where V0 is the equilibrium volume of the nanotube, given by Vo = 2 (R + h) R2 L; where L is the length, R is the radius and h is the shell thickness. The values of h used [38, 39] are 0.066 nm[38] and 0.34 nm[39]. Empirical force constant model [39] and non-orthogonal tight binding calculation [40] suggest that h should be 0.34 nm, a value used by us. The values of Young modulus of DWNTs show that there is not much signicance of inter-wall van der Waal interaction on the Young modulus. In Fig.6.14 we show the variation of Young modulus as a function of tube radius for dierent charge scheme on B/N atoms in BNNT. The following observations can be made. (i) BNNTs with charge 0.68e, 1.0e and 1.41e on B and N atoms have higher Young modulus than that of CNT whereas BNNTs with charge 0e and 0.41e on B and N atoms have smaller Young modulus as compared to CNT. This dierence has its origin in electrostatic interactions. (ii) Young modulus increases with the radius and asymptotically approaches the value of the at sheet, consistent with ref. 31 and 41, but in contrast to ref. 39. Young modulus of SWNTs and DWNTs are very similar which can be explained by the fact that the elastic properties of nanotubes are determined by the strength of B-N bonds in the bent BN sheet and C-C bonds in bent Graphene and also the vdW interaction is comparatively very weak in MWNTs.
209
Acknowledgement
This article is based on the set of lectures I gave in the SERC school at IIT, Guwahati in 2008. I have greatly benetted from the excellent set of lectures on the subject by Prof. David Kofke, University at Bualo, The State University of New York.
References
[1] P. S. d. Laplace, Theorie Analytique des Probabilities, 3 ed. (GauthierVillars, Paris). [2] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, Journal of Chemical Physics 21, 1087 (1953). [3] E. Fermi, J. Pasta, and S. Ulam, Document LA-1940 (1955). [4] B. J. Alder and T. E. Wainwright, J. Chem. Phys. 27, 1208 (1957). [5] J. B. Gibson, A. N. Goland, M. Milgram, and G. H. Vineyard, Physical Review 120, 1229 (1960). [6] A. Rahman, Physical Review a-General Physics 136 (2A), A405 (1964). [7] L. Verlet, Physical Review 165, 201 (1968). [8] L. Verlet, Physical Review 159, 98 (1967). [9] A. Rahman and F. H. Stillinger, J. Chem. Phys. 55, 3336 (1971). [10] G. D. Harp and B. J. Berne, Physical Review A 2, 975 (1970); G. D. Harp and B. J. Berne, Journal of Chemical Physics 49, 1249 (1968). [11] M. P. Allen and D. J. Tildesley, Computer Simulation of liquids. (Oxford: Clarendon, 1987). [12] D. Frenkel and B. Smit, Understanding Molecular Simulations. (Academic Press: San Diego 2002). [13] A. Leach, Molecular Modelling. (Prentice Hall 2001); T. Schlick, Molecular Modeling and Simulation - An Interdisciplinary Guide (Springer, New York, 2002). [14] M. Tuckerman, B. J. Berne, and G. J. Martyna, J. Chem. Phys. 97, 1990 (1992). [15] C. W. Gear, Report No. ANL-7126, 1966. [16] J. P. Ryckaert, G. Ciccotti, and H. J. C. Berendsen, J. Comput. Phys. 23, 327 (1977).
212
Chapter 6. References
[17] J. D. Honeycutt and H. C. Andersen, J. Phys. Chem. 90, 1585 (1986); J. D. Honeycutt and H. C. Andersen, Chem. Phys. Letts. 108, 535 (1984). [18] S. D. Stoddard, Phys Rev A8, 1504 (1973). [19] B. Quentrec and C. Brot, J. Comp. Phys 13, 430 (1973). [20] D. M. Heyes, Electrostatic potentials and elds in innite point charge lattice. (1981). [21] D. Fincham, Molecular Simulation 13, 1 (1994). [22] T. Darden, York D., Pedersen L., J Chem Phys 98, 10089 (1993); T. A. Darden, A. Toukmaji, and L. G. Pedersen, Journal De Chimie Physique Et De Physico-Chimie Biologique 94 (7-8), 1346 (1997). [23] S. Nose, J. Chem. Phys 81, 511 (1984); S. Nose, Molecular Physics 52, 255 (1984); W. G. Hoover, Physical Review A 31, 1695 (1985). [24] H. J. C. Berendsen, J. P. M. Postma, A. van Gunsteren, A. DiNola, and H. R. Haak, J. Chem. Phys. 81, 3684 (1984). [25] B. Mukherjee, P. K. Maiti, C. Dasgupta, and A. K. Sood, J. Chem. Phys. 126, 124704 (2007). [26] B. Mukherjee, P. K. Maiti, C. Dasgupta, and A. K. Sood, Acs Nano 2, 1189 (2008); B. Mukherjee, P. K. Maiti, C. Dasgupta, and A. K. Sood, J. Phys. Chem. B 113, 10322 (2009). [27] D. A. Case, Darden, T.A., Cheatham III, T.E., Simmerling, C.L.,Wang, J., Duke, R.E., Luo, R., Merz, K.M.,Wang, B.,Pearlman, D.A., Cowley, M., Brozell, S., Tsui, V., Gohlke, H., Mongan, J., Hornak, V., Cui, G., Beroza, P., Schafmeister, C., Caldwell, J.W., Ross, W.S., Kollman, P.A., AMBER 8 (University of California, San Francisco., 2004). [28] W. D. Cornell, P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz, D. M. Ferguson, D. C. Spellmeyer, T. Fox, J. W. Caldwell, and P. A. Kollman, JACS 117, 5179 (1995). [29] G. Hummer, J. C. Rasaiah, and J. P. Noworyta, Nature 414 (6860), 188 (2001); B. Mukherjee, P. K. Maiti, C. Dasgupta, and A. K. Sood, J. Chem. Phys. 126(2007). [30] B. Mukherjee, P. K. Maiti, C. Dasgupta, and A. K. Sood, Journal of Nanoscience and Nanotechnology 9, 5303 (2009). [31] E. Hernandez, C. Goze, P. Bernier, and A. Rubio, Phys. Rev. Letts. 80, 4502 (1998). [32] K. N. Kudin, G. E. Scuseria, and B. I. Yakobson, Physical Review B 64, (2001). [33] A. K. Rappe, Goddard, W. A. III, J. Phys. Chem 95, 3358 (1991).
213 [34] S. L. Mayo, B. D. Olafson, and W. A. Goddard, J. Phys. Chem. 94, 8897 (1990). [35] W. H. Moon and H. J. Hwang, Physica E-Low-Dimensional Systems & Nanostructures 23 (1-2), 26 (2004). [36] M. Santosh, P. K. Maiti, and A. K. Sood, Journal of Nanoscience and Nanotechnology 9, 5425(2009). [37] J. F. Waters, P. R. Guduru, M. Jouzi, J. M. Xu, T. Hanlon, and S. Suresh, App. Phys. Letts. 87 (2005). [38] B. I. Yakobson, C. J. Brabec, and J. Bernholc, Phys. Rev. Letts. 76, 2511 (1996). [39] J. P. Lu, Phys. Rev. Letts., 79, 1297 (1997). [40] A. Rubio, J. L. Corkill, and M. L. Cohen, Phys. Rev. B 49, 5081 (1994). [41] C. Y. Li and T. W. Chou, Journal of Nanoscience and Nanotechnology 6, 54 (2006). [42] P. H. Hunenberger, Adv. Polym. Sci., 173, 105(2005).