Molecular Dynamics Simulation

Chapter 3

Molecular Dynamics Simulation

Molecular dynamics calculates the biological system in liquid state; this property is the same as the liquid NMR. Besides, the time average property and kinetic property of molecular dynamics simulation is in good accordance with experiments.

In order to understand molecular motion, Lipari and Szabo ^{24; 25} developed the Model-Free approach connecting NMR relaxation and molecular motion. This approach also made NMR relaxation and MD simulation more comparable. Because MD simulates the system usually in ns-ps time scale, the Model-Free approach efficiently calculated magnitude of motion in the same time regime. Thus, this approach increased the connection of MD and NMR. And MD simulation is a good tool to analyze dynamics features in atomic view. That is the reason why we conducted MD simulation to assist us to understand the dynamics in catalytic pathway.

3.1 The Force Field Formalism

Ιn recent standard molecular dynamics simulation, molecular mechanics is invariably used to perform calculations on systems containing significant numbers of atoms. Molecular mechanics is based upon a rather simple model of the interaction within a system with contributions from processes such as the stretching of bonds, the opening and closing of angles and the rotations about single bonds.

Many of the molecular modeling force field in use today for molecular systems can be interpreted in terms of a relatively simple four-component picture of the intra- and inter- molecular force within the system. Energetic penalties are associated with

the deviation of bonds and angles away from their “reference” or “equilibrium” values, and there is also a function that describes how the energy changes as a bond rotates, and force field also contains terms that described the interaction between non-bonded parts of the system. An attractive feature of this representation is that the various terms can be ascribed to changes in specific internal coordinates such as bond lengths, angles, and the rotation of bonds or movements of atoms relative to each other. This makes it easier to understand how changes in the force field parameters affect its performance, and also helps in the parameterizations process. One functional form for such a force field that can be used to model single molecules or assemblies of atoms or molecule is³⁶:

( ) ( ) ( )

∑ ∑

∑

∑

∑

−

= =+ ⎟⎟

⎟

⎠

⎞

⎜⎜

⎜

⎝

⎛

⎥+

⎥

⎦

⎤

⎢⎢

⎣

⎡

⎟⎟

⎠

⎞

⎜⎜

⎝

−⎛

⎟⎟

⎠

⎞

⎜⎜

⎝ + ⎛

− +

+

− +

−

=

1

1 1 0

6 12

2 0 , 2

0 ,

4 4

) cos(

2 1 2

) 2 (

N

i N

i

j ij

j i ij

ij ij

ij ij

torsions n angles

i i i i

i bonds N i

r q q r

r

V n l k

k l r

U

πε σ

ε σ

γ ω θ

r

θ

(3.1)

) (r^N

U r denotes the potential energy which is a function of the set of the position of N particles (usually atoms). The various contributions are schematically represented in Figure 3.1. The first term in equation [3.1] describes the interaction between pairs of bonded atoms, modeled here by a harmonic potential that gives the increase in energy as the bond length l

rrs

i deviates from the reference value li.0. The second term is a summation over all valence angles in the molecule, also modeled by a harmonic potential. The third term is torsional potential that models how the energy changes as a bond rotates. The fourth term contribution is the non-bonded term. This calculation between all pairs of atoms (i and j) that are in different molecules or that are in the same molecule but separated by a least three bonds. In a simple force field the non-bonded term is usually modeled using a Coulomb potential term for electrostatic and a Lennard-Jones potential for van der Waals interactions.

Figure 3.1

Schematic representation of the four key contributions to a molecular mechanics force field: bond stretching, angle bending and torsional terms and non-bonded interactions.

3.2 Integration of the Newtonian Equation of Motion

For any arrangement of the atoms in the system, the force acting on each atom due to interactions with other atoms can be calculated by differentiating the energy function. From the force on each atom it is possible to determine its acceleration via Newton’s second law. Integration of the equations of motion should then yield a trajectory that describes how the positions, velocities and accelerations of the particles vary with time, and from which the ensemble average values of properties can be determined via the ergodic theorem:

∫

=

∞

=

→

=

^τ

τ

τ

⁰

^{( ( ), ( ))}

lim 1 )

(

t

N

N

t r t dt

p A t

A

A

(3.2)

The instantaneous value of the property A can be written as A(p^N(t),r^N(t)), where p^N(t) and r^N(t) represent the N momenta and positions, respectively, at time t.

In molecular dynamics, successive configurations of the system are generated by integrating Newton’s laws of motion. The result is a trajectory that specifies how positions and velocities of the particles in the system vary with time. For a system in the micro-canonical ensemble, a trajectory is obtained by solving the differential equations embodied in Newton’s second law (F=ma):

i i x

m F dt

x

d²₂ = _i (3.3)

This equation describes the motion of a particle of mass mi along one coordinate (xi) with Fxi being the force on the particle in that direction. The essential idea is that the integration is broken down into many small stages, each separated in time by a fixed time δ t. The total force on each particle in the configuration at a time t is calculated as the vector sum of its interactions with other particles. From the force we can determine the accelerations of the particles, which are then combined with the

positions and velocities at a time t to calculate the positions and velocities at a time t+ δ t . The force is assumed to be constant during the time step. The forces on the particles in their new positions are then determined, leading to new positions and velocities at time t+2δ t, and so on.

It is assumed that the positions and dynamic properties (velocities, accelerations, etc.) can be approximated as Taylor series expansions:

+ L +

+ +

=

+ t r t tv t t a t t b t

r

^{2 3}

6 ) 1 2 (

) 1 ( )

( )

( δ δ δ δ

where v is the velocity (the first derivative of the positions with respect to time), a is the acceleration (the second derivative), b is the third derivative, and so on. The Verlet algorithm ³⁷ is probably the most widely method for integrating the equations of motion in a molecular dynamics simulation. The Verlet algorithm uses the position and accelerations at time t, and the position from the previous step, r(t-δt), to calculate the new positions at t+δt, r(t+δt). We can write down the following relationships between these quantities and the velocities at time t:

+ L +

+

=

+ ( )

2 ) 1 ( )

( )

( t t r t tv t t

²

a t

r δ δ δ

(3.4)

− L +

−

=

− ( )

2 ) 1 ( )

( )

( t t r t tv t t

²

a t

r δ δ δ

(3.5)

Adding these two equations gives

r ( t + δ t ) = 2 r ( t ) − r ( t − δ t ) + δ t

²

a ( t )

(3.6) The velocities do not explicitly appear in the Verlet integration algorithm. The velocities can be calculated in a variety of ways; a simple approach is to divide the difference in positions at times t+δt and t-δt by 2δt:

v ( t ) = [ r ( t + δ t ) − r ( t − δ t )] / 2 δ t

(3.7)

Alternatively, the velocities can be estimated at the half-step, ^t δ^t 2 +1 :

v t

δ

t) [r(t

δ

t) r(t)]/

δ

t 2

( + 1 = + − (3.8)

Implementation of the Verlet algorithm is straightforward and the storage requirements are modest, comprising two sets of positions r(t-δt) and r(t) and the accelerations a(t).

3.3 Temperature and Pressure Coupling

Molecular dynamics is traditionally performed in the constant NVE (or NVEP) ensemble. Although thermodynamic results can be transformed between ensembles, this is strictly only possible in the limit of infinite system size (‘the thermodynamic limit’). It may therefore be desired to perform the simulation in a different ensemble.

The two most common alternative ensembles are the constant NVT and the constant NPT ensembles. In our system, we used the NPT ensembles.

In order to control the temperature of molecular system, we coupled the system to an external heat bath that is fixed at the desired temperature³⁸. The bath acted as a source of thermal energy, supplying or removing heat from the system as appropriate.

The velocities are scaled at each step, such that the rate of change of temperature is proportional to the difference in temperature between the bath and the system:

( ( ) 1

)

( T T t

dt t dT

bath

−

= τ )

(3.9)

where τ is a coupling parameter whose magnitude determines how tightly the bath and the system are coupled together. This method gives an exponential decay of the system towards the desired temperature. The change in temperature between

successive steps is:

(

T T(t) T = t _bath−

∆

^δ τ )

(3.10) If τ is large, then the coupling will weak; and if τ is small, the coupling will be strong.

This approach permits the system to fluctuate about the desired temperature.

A simulation in the isothermal-isobaric ensemble maintains constant pressure by changing the volume of the simulation cell. The amount of volume fluctuation is related to the isothermal compressibility, κ:

P T

V

V ⎟

⎠

⎜ ⎞

⎝

⎛

∂

− ∂

κ = 1 (3.9)

An easily compressible substance has a larger value of κ, so large volume fluctuations occur at a given pressure than in a more incompressible substance. Many of the methods used for pressure control are analogous to those used for temperature control.

Thus, the pressure can be maintained at a constant value by simply scaling the volume.

An alternative is to couple the system to a ‘pressure bath’, analogous to a temperature bath ³⁸. The rate of change of pressure is given by:

)) ( 1 (

)

( P P t

dt t dP

bath p

−

=

τ

(3.10)

τ_p is the coupling constant, Pbathis the pressure of the ‘bath’, and P(t) is the actual pressure at time t. The volume of the simulation box is scaled by a factor λ, which is equivalent to scaling the atomic coordinates by a factor λ^1/3. Thus:

) (

1

bath

p

P t P −

−

= τ

κ δ

λ

(3.11)

The new positions are given by:

r

i

′ = λ

^1/3

r

i (3.12) The constant κ can be combined with the relaxation constant τp as a single constant.

This expression can be applied isotropically or anisotropically. In general, it is best to use the anisotropic approach as this enables the box dimensions to change

independently.

In the extended pressure-coupling system methods, first introduced by Anderson

39 an extra degree of freedom, corresponding to the volume of the box, is added to the system. The kinetic energy associated with this degree of freedom is ( / )²

2

1Q dV dt ,

where Q is the ‘mass’ of the piston. The piston also has potential energy PV, where P

is the desired pressure and V is the volume of the system. The volume can vary during the simulation, with the average volume being determined by the balance between the internal pressure of the system and the desired external pressure. The extended-system temperature-scaling method uses a scaled time; in the extend pressure method the coordinates of the extended system are related to the ‘real’ coordinates by

(3.13)

i

i

V r

r ′ =

^1/3

3.4 Constrained Dynamics

When the simulation involves conformational flexible molecules then the motion is inevitably desired in terms of the atomic Cartesian coordinates. The high frequency motions (e.g. the vibration of bonds between hydrogen atoms and heavy atoms) are usually of less interest than the lower frequency modes, which often correspond to major conformational changes. Unfortunately, the time step of a molecular dynamics simulation is dictated by the highest frequency motion present in the system. In order to increase the time step without affecting the other internal degrees of freedom, the constrained dynamics is applied in many biological systems.

The most commonly used method for applying constraints, particularly in molecular dynamics, is the SHAKE procedure of Ryckaert, Ciccotti and Berendsen ⁴⁰. In constraint dynamics the equations of motion are solved while simultaneously satisfying the imposed constraints. Constrained systems have been much studied in classical mechanics. Constraints are often categorized as holonomic or non-holonomic.

Holonomic constraints can be expressed in the form

0 ) , , , ,

( q

₁

q

₂

q

₃

t =

f K

(3.14)

q₁, q2, etc., are the coordinates of the particles. Non-holonomic constraints can not be expressed in this way. SHAKE uses holonomic constraints. In a constrained system the coordinates of the particles are not independent and the equations of motion in each of the coordinates of the particles are not independent and the equations of motion in each of the coordinate directions are connected. As we know, the motion of a system of N particles can be described in terms of 3N independent coordinates or degree of freedom. If there are k holonomic constraints then the number of degrees of freedom is reduced to 3N-k.

In the general case, the equations of motion for a constrained system involve two type of force: the ‘normal’ forces arising from the intra- and intermolecular interactions, and the forces due to the constraints. The constraint σk requires the bond between atoms i and j to remain fixed. The constraints influences the Cartesian coordinates of atom i and j. The force due to this constraint can be written as follows:

F

_{ckx k}

x

^k

∂

= λ ∂ σ

(3.15)

where λk is the Lagrange multiplier and x represents one of the Cartesian coordinates of the two atoms. If an atom is involved in a number of constraints then the local constraint force equals the sum of all such terms. The nature of the constraint for a bond between atoms i and j is:

(3.16)

( − )

²

−

²

= 0

=

_{i j ij}

ij

r r d

σ

The constraint force lies along the bond at all times. For each constraint bond, there is an equal and opposite force on the two atoms that comprise the bond. The overall effect is that the constraint forces do no work. Suppose the constraint k corresponds to the bond length with respect to the coordinates of atoms i and j and multiplying by an undetermined multiplier,

) (

and ) (

2

) (

so ) (

2

j i cj

j i j

k

j i ci

j i i

k

r r F

r r r

r r F

r r r

−

−

=

−

−

∂ =

∂

−

=

−

∂ =

∂

σ λ σ λ

(3.17)

The above expression for the forces can be incorporated into the Verlet algorithm as follows:

( )

2 ( ) ( ) ² ( ) ²r (t) m

t t m F t t t r t r t t

r _ij

k i

k i

i i

i

i ⁺

^δ

^{= − −}

^δ

⁺

^δ

⁺

∑ ^{λ δ}

(3.18)

The SHAKE method uses an alternative approach in which constraint is

considered in turn and solved. Satisfying one constraint may cause another constraint to be violated, and so it is necessary to iterate around the constraints until they are all satisfied to within some tolerance. The tolerance should be tight enough to ensure that the fluctuations in the simulation due to the SHAKE algorithm are much smaller than the fluctuations due to other sources, such as the use of cutoffs. Another important requirement is that the constrained degree of freedom should be only weakly coupled to the remaining degree of freedom, so that the motion of the molecule is not affected by the application of the constraints. The sampling of unconstrained degree of freedom should also be unaffected. Therefore, angle constraints can be easily accommodated in the SHAKE scheme by recognizing that an angle constraint simply corresponds to an additional distance constraint. The SHAKE method has been extended by Tobias and Brooks ⁴¹ to enable constraints to be applied to arbitrary internal coordinates. This enables the torsion angle of a rotatable bond to be constrained to a particular value during a molecular dynamics simulation.