Schlick T. Molecular Modeling and Simulation: An Interdisciplinary Guide

Подождите немного. Документ загружается.

13.4. The Verlet Algorithm 451

Though the position update formula of eq. (13.10) is fourth-order accurate, the

Verlet method is overall second-order accurate. This is because eq. (13.10)di-

vided by Δt

represents an O(Δt

) approximation to

X =



F (X). Though the

same trajectory positions are followed in theory, differences are realized in prac-

tice among Verlet variants due to ﬁnite computer arithmetic. Local accuracy

together with stability considerations affect global measurements, like the kinetic

temperature (as computed from a simulation by the relation of eq. (13.4)). This is

because such global properties must be interpreted statistically.

13.4.2 Leapfrog, Velocity Verlet, and Position Verlet

Three popular variants of Verlet are the leapfrog scheme [553], velocity Ver-

let [1244], and position Verlet [95, 99, 1277]. All originate from the same

leapfrog/St¨ormer/Verlet method [1227, 1302] and have favorable numerical

properties like symplecticness (e.g., [1199]) and time-reversibility [1277].

The leapfrog scheme is so called because the velocity is deﬁned at half steps

n+1/2

) while positions are deﬁned at whole steps (X

)ofΔt. The propaga-

tion formula of leapfrog thus transforms one point in phase space to the next as

n−1/2

} =⇒{V

n+1/2

n+1

}. The velocity Verlet scheme transforms

instead {X

} =⇒{X

n+1

We ﬁrst present these three methods using a symmetric propagation triplet.

Both leapfrog and velocity Verlet can be written as:

n+1/2

= V

Δt



n+1

= X

+ΔtV

n+1/2

(13.14)

n+1

= V

n+1/2

Δt



n+1

The position Verlet scheme can be written as:

n+1/2

= X

Δt

n+1

= V

+Δt



n+1/2

(13.15)

n+1

= X

n+1/2

Δt

n+1

Note that the force is evaluated at the endpoints in leapfrog and velocity Verlet

but at the midpoint of position Verlet.

Leapfrog

The leapfrog Verlet method as described above can also be written as

n+1/2

= V

n−1/2

+Δt



, (13.16)

n+1

= X

+ΔtV

n+1/2

. (13.17)

452 13. Molecular Dynamics: Basics

It can be derived from the position propagation formula of eq. (13.10)in

combination with the velocity formula (13.13) approximated at half steps:

n−1/2

=(X

− X

n−1

)/Δt. (13.18)

It is transparent to see the equivalence of the position update formula of

leapfrog (eq. (13.17)) to that used in the triplet scheme (eq. (13.14)). To see that

the velocity updates are equivalent as well, replace the n +1superscript by n in

the third equation of the triplet (13.14)toget:

= V

n−1/2

Δt



, (13.19)

and add this equation to the ﬁrst equation (for V

n+1/2

) of the triplet (13.14).

The result V

n+1/2

= V

n−1/2

+Δt



is the velocity formula of leapfrog

(eq. (13.16)).

Velocity Verlet

Besides the triplet equation (13.14), the velocity Verlet scheme can also be

written as:

n+1

= X

+ΔtV

Δt



, (13.20)

n+1

= V

Δt

(



n+1

) . (13.21)

The ﬁrst equation can be derived by using eq. (13.13) to express X

n−1

n+1

− 2ΔtV

) and plugging this substitution for X

n−1

in eq. (13.10). The

second equation of velocity Verlet is obtained by applying the velocity equa-

tion (13.13) to express V

n+1

in terms of positions and then using eq. (13.20)

to write X

n+2

in terms of X

n+1

, V

n+1

,and



n+1

,andX

in terms of X

n+1

,and



. That is, we write V

n+1

n+2

− X

2Δt



n+1

+ΔtV

n+1

Δt



n+1

2Δt



−



n+1

− ΔtV

−

Δt



2Δt



n+1

+ V

Δt

(



n+1

) . (13.22)

Collecting the two V

n+1

terms yields

n+1

Δt

(



n+1

) , (13.23)

leading immediately to the position update formula of eq. (13.21).

13.5. Constrained Dynamics 453

The equivalence of the triplet version (eq. (13.14)) of velocity Verlet to the form

in eqs. (13.20), (13.21) is straightforward.

Position Verlet

Besides the triplet (13.15), we can express position Verlet as:

n+1

= X

+ΔtV

Δt



n+1/2

, (13.24)

n+1

= V

+Δt



n+1/2

. (13.25)

MTS Preference

The velocity Verlet scheme has been customarily used in MTS schemes because

the forces are evaluated at the ends of the interval and this is more amenable to

force splitting. Moreover, it only requires one constrained dynamics application

per inner timestep cycle (when position is updated) when constrained dynam-

ics techniques are used; two such constrained-dynamics iterations are required

for position Verlet. However, position Verlet — ﬁrst suggested to be preferable

in large-timestep methods [95] — was recently shown to have advantages in

connection with moderate to large timesteps [97,99].

13.5 Constrained Dynamics

In constrained dynamics, the equations of motion are augmented by algebraic

constraints to freeze the highest-frequency interactions. To illustrate, consider a

typical modeling of bond-length potentials as the harmonic form

bond

− r

)

where r

is an interatomic distance with equilibrium value r

,andS

is a force

constant. The constraint implemented in the context of MD simulations is thus

= r

− r

=0.

Using the formalism of Lagrange multipliers, we have then in place of eq. (13.3)

the following system to solve:

V (t)=−∇E(X(t)) − g



(X(t))

Λ ,

X(t)=V (t) , (13.26)

g(X(t)) = 0 .

Here g(X(t)) is a vector with entries g

containing the individual constraints, and

Λ is the vector of Lagrange multipliers.

454 13. Molecular Dynamics: Basics

SHAKE

In 1977, the SHAKE algorithm was introduced on the basis of the leapfrog Verlet

scheme of eqs. (13.16, 13.17)[1079]:

n+1/2

= V

n−1/2

− Δt∇E(X

) − Δtg



)

n+1

= X

+ΔtV

n+1/2

, (13.27)

g(X

n+1

)=0.

Like the Verlet method, the forces −∇E must be computed only once per step.

However, the m constraints require at each timestep, in addition, the solution of

a nonlinear system of m equations in the m unknowns {λ

}. Thus, the SHAKE

method is semi-explicit.

RATTLE

A variant termed RATTLE was later introduced [44]. Similar techniques have

been proposed for constraining other internal degrees of freedom [1264], and

direct methods have been developed for the special case of rigid water mole-

cules [868]. Leimkuhler and Skeel [733] have shown that the RATTLE method

is symplectic and that SHAKE, while not technically symplectic, yields solu-

tions identical to those of RATTLE for the positions and only slightly perturbed

velocities.

Computational Advantage

With the fastest vibrations removed from the model, the integration timestep can

be lengthened, resulting in a typical force calculation per time ratio of 0.5 per

fs. This computational advantage must be balanced against the additional work

per step required to solve the nonlinear equations of constraints. Each proposed

constrained model is thus accompanied by a practical numerical algorithm.

For the time discretization of eq. (13.27), an iterative scheme for solution of the

nonlinear equations was presented in [1079], where the individual constraints are

imposed sequentially to reset the coordinates that result from an unconstrained

Verlet step. This iterative scheme is widely used due to its simplicity and frugal

consumption of computer memory. However, the SHAKE iteration can converge

very slowly, or even fail to converge,for complex bond geometries [92]. Improved

convergence was later reported by Barth et al. [92], who proposed enhancements

based on the relationship between the SHAKE iteration process and the iterative

Gauss-Seidel/Successive-Over-Relaxation techniques for solving nonlinear sys-

tems. This was also independently developed by Xie and Scott, as reported later

by Xie et al. [1401].

Limitations

The computational advantage of constrained models is clear, but the verity of

the constrained model for reproducing the dynamics of unconstrained systems

is a separate issue. Van Gunsteren and Karplus [1296] showed through simula-

tions of the protein BPTI in vacuum that the use of ﬁxed bond lengths does not

13.6. Various MD Ensembles 455

signiﬁcantly alter the dynamical properties of the system, whereas ﬁxing bond

angles does. Similar conclusions were reported by Toxvaerd [1269] for decane

molecules.

Though it is possible that the former study was also inﬂuenced by poor conver-

gence of SHAKE for bond-angle constraints (which can lead to overdetermination

and singularities in the constraints equations), the intricate vibrational coupling in

biomolecules argues against the general use of angle constraints. As shown in

[93], from the point of view of increasing the timestep in biomolecular simula-

tions, only constraints of light-atom angle terms in addition to all bond lengths

might be beneﬁcial. Furthermore, it is common practice not to constrain motion

beyond that associated with light-atom bonds because the overall motion can be

disturbed due to vibrational coupling.

13.6 Various MD Ensembles

13.6.1 Need for Other Ensembles

The algorithmic framework discussed thus far is appropriate for the microcanoni-

cal (constant NVE) ensemble (see the equilibration illustration in Figure 13.2,for

example), where the total energy E is a constant of the motion. This assumes that

the time averages are equivalent to the ensemble averages.

The study of molecular properties as a function of temperature and pressure —

rather than volume and energy — is of general importance. Thus, microcanonical

ensembles are inappropriate for simulating certain systems which require constant

pressure and/or temperature conditions and allow the energy and volume to ﬂuctu-

ate. This is the case for a micelle or lipid bilayer complex (proteins and lipids in an

ionic solution) under constant pressure (isotropic pressure or the pressure tensor),

and for crystalline solids under constant stress (i.e., constant pressure tensor). For

such systems, other ensembles may be more appropriate, such as canonical (con-

stant temperature and volume, NVT) isothermal-isobaric (constant temperature

and pressure, NPT), or constant pressure and enthalpy (NPH). Special techniques

have been developed for these ensembles, and most rely on the Verlet framework

just described, including the MTS variants discussed in the next chapter [836].

For pioneering works, see the articles by Andersen [43] and the extension by

Parrinello and Rahman [964] for NPH, and that of Nos´e[925] and Hoover [566]

for NVT and NPT ensembles (started in 1984). For details on algorithms for these

extended ensembles, see the texts [22,428,731] and current literature.

In this section we only give a ﬂavor of some of these methods. Methods for han-

dling these thermodynamic constraints or generating the various ensembles can be

grouped into simple constrained formulations, including approaches that involve

stochastic-motivated models, and more sophisticated techniques that involve ad-

ditional degrees of freedom (extended system methods). The former group, while

easy to implement, does not generally generate the desired ensemble in rigorous

terms.

456 13. Molecular Dynamics: Basics

13.6.2 Simple Algorithms

Constraint-based methods can involve Lagrange-multiplier coupling, or simple

scaling of phase-space variables. For example, the velocity vector might be simply

scaled at each timestep to ﬁx the desired kinetic temperature T at the target value

[22]. This can be achieved by the scaling

new

← c

old

, (13.28)

where



/T . (13.29)

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 artifactual pumping of energy into low-frequency

modes [524]; occasional resetting of velocities can mitigate this problem.

Weak Coupling Thermostat for Constant T

Similar algorithms that are motivated by stochastic approaches are also easy to

implement but can steer the system toward the desired constant (e.g., target tem-

perature or pressure) more gently. This can be accomplished by mimicking a

diffusive process governed by a ﬁctitious frictional coefﬁcient γ, whose value

controls the relaxation rate of this coupling [123].

For example, consider the weak coupling thermostat proposed by Berendsen

and coworkers [123], widely used for constant-temperature MD. In this approach,

the equations of motion (eq. (13.3)) are modiﬁed as:

X(t)=V (t) ,

V (t)=−∇E(X(t)) −γ

MV (t) , (13.30)

2τ



1 −



where γ

has units of inverse time, T is the instantaneous kinetic temperature, and

τ is the time constant of the coupling to the heat bath. (Compare to the Langevin

formalism described in the next chapter, Section 14.4).

This augmentation effectively scales the velocity vector via eq. (13.28)bythe

factor:



1 −

Δt



1 −



, (13.31)

where Δt is the (Verlet) timestep. This scaling produces a least-squares lo-

cal disturbance to the velocity vector, so as to satisfy the global temperature

constraint.

13.6. Various MD Ensembles 457

The advantage of this weak bath coupling approach is the introduction of added

thermostat parameters (not dynamic variables, as below). The factor τ controls the

characteristic decay time — and hence the strength — of the coupling to the heat

bath. When τ is large, γ

is small, and the scaling factor c

approaches unity (the

microcanonical ensemble is approached); when τ is small, γ

is large — that is,

the coupling to the heat bath is strong — and the energy exchange between the

molecular system and the thermal reservoir is signiﬁcant.

In particular, when τ =Δt, c



/T, producing the simplest scaling

approach for constant-T simulations, but also the most artiﬁcial, as mentioned

above. The advantage of the ﬂexible coupling τ parameter (or γ

) is the control of

the rate at which the target temperature is reached.

Still, this approach, while convenient, does not generate the true canonical en-

semble; it can also lead to artifacts as simple velocity scaling [524]. See below for

a more rigorous approach that reproduces the desired thermal ensemble properly.

Weak Coupling Barostat for Constant P

In the same spirit, for the purpose of constant pressure simulations for isotropic

systems modeled in a ﬁnite volume (generally via periodic boundary conditions),

Berendsen and coworkers suggest a simple scaling to the Cartesian positions,

new

← c

old

, (13.32)

as well as a volume scaling (see below), where the scale factor is:



1 −

βΔt

− P )



1/3

. (13.33)

Here β is the isothermal compressibility, P is the instantaneous(or internal) pres-

sure, P

is the target pressure (or external pressure P

), and τ is the pressure

coupling time (see Box 13.2). The goal of this formulation is to allow the vol-

ume of the macromolecular system (ν

) to ﬂuctuate as the instantaneous (internal)

pressure (P ) approaches the applied (external) pressure, P

This procedure is applied to the coordinates of particles under periodic bound-

ary conditions, and to the box length L for an isotropic system in a cubic box

new

← c

old

). From the point of view of the guiding dynamic equations, this

coordinate and volume scaling protocol effectively solves the augmented equation

of motion for the time derivative of position (

X(t)=V (t)), where the additional

term is proportional to X as follows:

X(t)=V (t) −



βX(t) ,

(13.34)



β =

β (P

− P (t))

3τ

458 13. Molecular Dynamics: Basics

This equation can also be written as a pair of differential equations describing the

evolution of the position (X)andvolume(ν

)byusingeq.(13.36)inBox13.2 to

produce the following equations (we suppress the time dependency for clarity):

X = V +

˙ν

(13.35)

˙ν

β (P − P

)

Thus, this approach to constant pressure simulations allows the volume of the

system to ﬂuctuate as the pressure is held constant by changing the cell size uni-

formly but not its shape, as appropriate for an isotropic system. For an anisotropic

system, instead of the scalar pressure variable, the 3 × 3 pressure matrix (see

Box 13.2) is relevant [853].

As for the weak coupling thermostat method, the above approach does not

generally produce trajectories from any known ensemble. A rigorous approach us-

ing an extended system method which includes dynamic thermostat and barostat

variables is outlined below. See also [382] for an improved constant-pressure sim-

ulation protocol using a Langevin piston that is straightforward to implement. The

Langevin piston method attempts to eliminate the nonphysical vibrations of the

volume associated with the approach above and to generate in theory the correct

NPT ensemble.

Box 13.2: Pressure Variables Deﬁnitions

The isothermal compressibility describes the change in volume (ν

) of a substance under

external pressure P via:

−1

βν



dν



(13.36)

P = − ˙ν

/ν

. (13.37)

The instantaneous pressure P is calculated from the thermodynamic relation between the

pressure, temperature, volume, and internal virial (vir)ofasystemby[22]:

P =

3 ν

− vir) . (13.38)

The internal virial is proportional to the inner product of the each atom’s position vector

) with the corresponding force component acting on atom i due to all particles (F

vir = −



) . (13.39)

(Intramolecular forces are thus considered in addition to intermolecular forces). Note that

in the expression for the pressure (eq. (13.38)) we have used the relation between the

temperature T and the kinetic energy E

(E

 =

T), where the number of degrees

of freedom, N

, is three times the number of atoms in the system.

13.6. Various MD Ensembles 459

For an anisotropic system, the pressure becomes a 3 × 3 matrix called the pressure tensor

P [853]. This tensor is deﬁned as:

P =





)



, (13.40)

where m

and v

are the scalar mass and 3-component velocity vector corresponding to

atom i, respectively; the result of an outer product between a 3 ×1 vector and its transpose

is a 3 × 3 matrix.

13.6.3 Extended System Methods

Extended-system methods introduce additional degrees of freedom to represent

the environment of the macromolecular system (e.g., pressure piston, thermostat,

etc.). These degrees of freedom have positions, momenta, and/or other associ-

ated variables. Examples include an external volume variable ν

associated with

a piston of given mass, whose potential energy is expressed in terms of the de-

sired pressure [924], or effective thermodynamic variables that are functions of

the positions and momenta of the added dynamic variables.

The equations of motion for the extended system are expressed as augmented

versions of the standard equations of motions, to represent the evolution of both

the internal and external variables in the desired ensemble. The solutions for the

time evolution of the extended system by standard numerical methods (for the

microcanonical ensemble) then produce the phase-space variables appropriate for

the desired ensemble.

Canonical Ensemble

For example, to produce the true canonical (NVT) ensemble — which the simple

scaling and thermostat approaches above do not accomplish — the method termed

Nos´e-Hoover [566,924] adds a ﬁctitious degree of freedom to the physical system

with ‘coordinate’ parameter x

(effectively a scaling parameter [428]), mass m

and thermodynamic friction coefﬁcient ζ

. (This friction coefﬁcient is related to

and the corresponding momentum ˙x

After appropriate scaling of variables (see [428], for example), the equations of

motion for the real and artiﬁcial system combine to yield:

X(t)=V (t) ,

V (t)=−∇E(X(t)) − ζ

MV (t) , (13.41)

(t)=2V

MV − gk

where g is the number of degrees of freedom in the system. The conserved quan-

tity under these augmented equations of motion has two energy terms in addition

to the internal kinetic and potential energies:



NVT

MV )+E(X)+

)+gk

. (13.42)

460 13. Molecular Dynamics: Basics

The choice of the ﬁctitious mass m

should ensure that the thermalization

process is efﬁcient. Too small a value implies large harmonic motion and rapid

thermal ﬂuctuations for the extended degree of freedom, and this in turn restricts

the timestep too severely. Too large a value, on the other hand, makes the ther-

malization process slow. In general, m

is chosen to be proportional to gk

and the integration of the above system is performed fairly accurately so that the

energy is well conserved. See the original literature and [428,835] for algorithmic

details.

Isothermal-Isobaric Ensemble

Controlling the pressure in MD simulations is more involved. Instead of ﬁxing the

volume, as in microcanonical simulations, the volume of the system is considered

as a dynamic variable, and it is allowed to ﬂuctuate while the pressure is held

constant. An approach analogous the Nos´e-Hoover NVT form above has been

described by Andersen [43] and Hoover [566] for the isothermal-isobaric, or NPT

ensembles.

Essentially, in addition to the effective coordinate, mass, and friction set

,ζ

) associated with the ﬁctitious thermostat variable, a set (x

,ζ

)

is introduced and associated with a virtual pressure piston (“barostat”).

The effective equations of motion for a 3-dimensional system become:

X(t)=V (t)+ζ

V (t)=F (X(t)) − M

V (t)





+ ζ



˙ν

=3ν

, (13.43)

(t)=3ν

− P

(2V

MV ) −m

(t)=2V

MV +

− (g +1)k

where g is the number of degrees of freedom in the system, P

is the external

applied pressure, and P

is the internal pressure, deﬁned as:

3 ν



− vir −



3 ν



∂E(X, ν

)

∂ν



; (13.44)

the virial vir is deﬁned in equation (13.39). The conserved quantity under these

augmented equations of motion is:



NPT

MV )+E(X, ν

)+(g +1)k

(13.45)

By this approach, the volume of the system ﬂuctuates under the applied thermostat

and barostats so that the system is driven to the steady state at which the average

internal pressure P is equal to the external applied force P