Integration Schemes

Appendix E

Integration Schemes E.1 Higher-Order Schemes The basic idea behind the predictor-corrector algorithms is to use information about the position and its first rt derivatives at time t to arrive at a prediction for the position and its first rt derivatives at time t + At. We then compute the forces (and thereby the accelerations) at the predicted positions. And then we find that these accelerations are not equal to the values that we had predicted. So we adjust our predictions for the accelerations to match the facts. But we do more than that. On the basis of the observed discrepancy between the predicted and observed accelerations, we also try to improve our estimate of the positions and the remaining r t - 1 derivatives. This is the "corrector" part of the predictor-corrector algorithm. The precise "recipe" used in applying this correction is a compromise between accuracy and stability. Here, we shall simply show a specific example of a predictorcorrector algorithm, without attempting to justify the form of the corrector part. Consider the Taylor expansion of the coordinate of a given particle at time t + At:

~r

r(t + At) -- r ( t ) + Atxv + 0t

At 2 02r At 3 03r + 2! 0t 2 3! 0t 3

~

~

~

~

Appendix E. Integration Schemes

534 Using the notation x0(t)

=

x~ ( t )

-

r(t) Or at--

Ot

x2(t)

-=

At 2 02r 2! 0t 2

x3(t)

--

At 3 03r 3! Ot3'

we can write the following predictions for x0(t + At) through x3(t + At):

x0(t+At)-

x0(t)+

xl(t+At)

=

x2(t+At)

=

x3(t + At)

=

Xl(t)if-

xl(t)

+

x2(t)

q-

x3(t)

2x2(t)

+

3x3(t)

x2(t)

+

3x3(t) x3 (t).

Now that we have x0 (t + At), we can compute the forces at the predicted position, and thus compute the corrected value for x2 (t + At). We denote the by Ax2: difference between X~~ and x predicted 2 AX2 D X~~

D xPredicted.

We now estimate "corrected" values for x0 through x3, as follows: xcorrected predicted rt -- xrt + CrLAx2,

(E.1.1)

where the C~ are constants fixed for a given order algorithm. As indicated, the values for C~ are such that they yield an optimal compromise between the accuracy and the stability of the algorithm. For instance, for a fifth-order predictor-corrector algorithm (i.e., one that uses x0 through x4), the values for C~ are Co

=

C1

--

C2

-

19 120 3 4 (of course)

1

1 C3

=

C4

--

2

1 ~. 12

E.2 Nos)-Hoover Algorithms

535

One may iterate the predictor and corrector steps to self-consistency. However, there is little point in doing so because (1) every iteration requires a force calculation. One would be better off spending the same computer time to run with a shorter time step and only one iteration because (2) even if we iterate the predictor-corrector algorithm to convergence, we still do not get the exact trajectory: the error is still of order At" for an rtth-order algorithm. This is w h y we gain more accuracy by going to a shorter time step than by iterating to convergence at a fixed value of At.

E.2 Nos4-Hoover Algorithms As discussed in section 6.1.2, it is advantageous to implement the Nos6 thermostat using the formulation of Hoover, equations (6.1.24)-(6.1.27). Since the velocity also appears on the right-hand side of equation (6.1.25), this scheme cannot be implemented directly into the velocity Verlet algorithm (see also section 4.3). To see this, consider a standard constant-N,V,E simulation, for which the velocity Verlet algorithm is of the form r(t + At) v(t+At)

--= r(t) + v ( t ) A t + f(t)At2/(2m) f(t + At) + f(t) = v(t)+ At. 2m

When we use this scheme for the Nos4-Hoover equations of motion, we obtain for the positions and velocities ri(t + At)

=

q ( t ) + v(t)at + [fi(t)/m{ - t(t)vi(t)] At2/2

v t ( t + At)

---- vt(t) + [fi(t + A t ) / m i -- t ( t + At)v~(t + At) + f t ( t ) / m i -- t(t)v~(t)] A t / 2 .

(E.2.1) (E.2.2)

The first step of the velocity Verlet algorithm can be carried out without difficulty. In the second step, we first update the velocity, using the old "forces" to the intermediate value v(t + At/2) = v'. And then we must use the new "forces" to update v': v t ( t + At) = v't + [ft(t + A t ) / m t -- t ( t + At)v~(t + At)] A t / 2 .

(E.2.3)

In these equations vi (t + At) appears on the right- and left-hand sides; therefore, these equations cannot be integrated exactly.1 For this reason the Nos6Hoover method is usually implemented using a predictor-corrector scheme or solved iteratively [138]. This has a disadvantage that the solution is no longer time reversible. Martyna et al. [85] have developed a set of explicit reversible integrators using the Liouville approach (see section 4.3.3) for this type of extended systems. 1For the harmonic oscillator it is possible to find an analytic solution (see Case Study 12).

536

E.2.1

Appendix E. Integration Schemes

Canonical E n s e m b l e

For M chains, the Nos6-Hoover equations of motion are given by (see also section 6.1.3) ti

=

pi/mi

f,~ -

F ~ - - P~~ ,, ,~ _

Lk

=

P~,k

15~,,

-

k-1,...,M

Qk

p~/mi-LkB-1-

[Qk-i

19~.k

kBT

LQG:'

--Q~2~,~,,

PLk+I. --

~

P

~

Qk+l

k

.

The Liouville operator for the equation of motions is defined as (see section 4.3.3) 0 iL - fl On with 11 = (r TM, pN, LM p ~ ) . Using the equations of motion, pi - mivi, and P ~.k -- QkV~.k, we obtain as Liouville operator for the Nos6-Hoover chains iLNHc

-

E vi" Vr~ + E

9Vv~ -

i=1

i=1

1TI.i

0

M-1

+ E

va, vi" Vv~ + i=1

Yak 0Lk k--1

0

(Gk -- vL~v~+,) 0 v ~ + G M 0vL~

k=l

with

G1 Gk

=

Q1

-

Qk

-

i=1

1

As explained in section 4.3.3, the Liouville equation combined with the Trotter formula is a powerful technique for deriving a time-reversible algorithm for solving the equations of motion numerically. Here we will use this technique to derive such a scheme for the Nos6-Hoover thermostats. We use a simplified version; a more complete description can be found in ref. [85].

537

E.2 Nos~-Hoover Algorithms

We have to make an intelligent separation of the Liouville operator. The first step is to separate the part of the Liouville operator that only involves the positions (iG) and the velocities (iL~) from the parts that involve the Nos6-Hoover thermostats (iLc): iLNHc - iL,. + iL,, + iLc with N

iLr

--

5 Vi " t=l

Vri

N

iLv

E Fi (ri) 9 Vvi 1Tti i=1

iLc

M

0 Y - vL~ 0Lk

N

~ - va, vi.

k=l

i=1

Vv~

M--1 c3 c3 E -I-GM + (Gk VLkVLk+,) 0VLk 0VLM k=l

There are several ways to factorize iLNHc USing the Trotter formula; we follow the one used by Martyna et al. [85]" e (iLAt) - - e ( i L c A t / 2 ) e ( i L ~ A t / 2 ) e ( i L ~ A t ) e ( i L ~ A t / 2 ) e ( i L c A t / 2 )

+ (.9 ( A t 3) .

(E.2.4)

The Nos6-Hoover chain part Lc has to be further factorized. Here, we will do this for a chain of length M = 2; the more general case is discussed in ref. [85]. The Nos6-Hoover part of the Liouville operator for this chain length can be separated into five terms: iLc - iL L + iLc~ + iLcl + iLvL, + iLG2, where the terms are defined as 2

iL L =

0 Y- vL~ OLk

k=l

N

iLc~

-

-Eva'v~'~Tv~ i=1

O

iLcl

=

G1 Ova,

iLvL,

=

--(VL, VL2)

iLG2

=

G28VL2.

O

8VL~

538

Appendix E. Integration Schemes

The factorization for the Trotter equation that we use is 2 e(iLc 2At/4)e(iLva, At/4+iLc, At/4)

e(iLcAt/2)

X e (iLa At/2) e(iLcv At/2) e(iEc, At/4+iLv a 1At/4) e(iLc 2At/4) e(iLc2At/4) [e(iL~.,At/8)e(iLc,At/4)e(iL~,At/8)] X e (iLLAt/2) e(iEc~ At/2) X [e(iLva, At/8)e(iLG'At/4)e(iL~a, At/8)] e (iLc2At/4). (E.2.5) L ]

Our numerical algorithm is now fully defined by equations (E.2.4) and (E.2.5). This seemingly complicated set of equations is actually relatively easy to implement in a simulation. To see how the implementation works, we need to know how each operator works on our coordinates rl = (r TM, v TM, ~,1 ,v~.,, ~2~V~2). If we start at t = 0 with initial condition rl, the position at time t = At follows from eiLNHcAtf[rN,v TM L1 VL, ~2,V~. 2] ~

~

~

Because of the Trotter expansion, we can apply each term in iLNHc, sequentially. For example, if we let the first term of the Liouville operator, iLG2, act on the initial state 11, exp

At 0 ) f [r N , pN , ~l,Vs163 -~--G2 0v&2 OO

=

=

E

n:0

]

11.

(G2at/4)

~!

a~

OV---~2f[rN'

f [ r N pN ~.1 V~., ~2

V ~ 2 "3L

pN ~1,VL,, ~2, VE.2] G2At/4]

This shows that the effect of iLc2 is to shift v~2 without affecting the other coordinates. This gives as transformation rule for this operator: e (iLc2At/4) ;

Vs 2 --4 Va2 + G2At/4.

(E.2.6)

The operators (iLva, and iLcv) are of the form exp (ax0/0x); such operators give a scaling of the x coordinate: exp ( ax O / f(x)

-

0 /f{exp[ln(x)]} exp ( a Oln(x)

=

f{exp[ln(x) + a ] } - f[xexp(a)]

2The second factorization, indicated by [... ], is used to avoid a hyperbolic sine function, which has a possible singularity. See ref. [85] for details.

E.2 Nos&HooverAlgorithms

539

If we apply this result 3 to iLia,, we obtain for this operator exp

( At

i3 ) f [ r N pN ~,,,V~,,,~,2,V~,2 ] ----8 va2va' Ova, ' '

= f[rN,pN,~,I,exp(---~VL2) Va,,&i,Va2], giving the transformation rule e(iLv~, At/S) .9

va, --+ exp

[-vazAt/8] va,.

(E.2.7)

In a similar w a y we can derive for the other terms e (iLc'at/4) : e (iLLAt/2)

e (iLcvAt/2)

v,., :

:

--+ va, + G1At/4

(E.2.8)

~1

-+

},,, --V~,, At/2

(E.2.9)

~,2

--+ }-2 - vaz A t / 2 --+ exp [-v~,, At/2] vi.

(E.2.10) (E.2.11)

vi

Finally, the transformation rules that are associated to iLv and iL~ are similar to the velocity Verlet algorithm, i.e., e (iLvAt/2)

:

e (iL~At) :

vt

--+ v~ + F~At/(2rrt)

(E.2.12)

r~

--+ rt + v~At.

(E.2.13)

With these transformation rules (E.2.6)-(E.2.13) we can write d o w n our numerical algorithm by subsequently applying the transformation rules according to the order defined by equations (E.2.4) and (E.2.5). If we start with initial coordinate 11(0) = ( r N , v N , ~,1, V~,,, ~2, VL2), we have to apply first e iLc r I. Since this operator is further factorized according to equation (E.2.5), the first step in our algorithm is to apply e ( i L G 2 A t / 4 ) . According to transformation rule (E.2.6) applying this operator on 11 gives as new state v~,2 (At/4) - v,,2 + G2At/4. The o u t p u t of this rule is the new state on which we apply the next operator in equation (E.2.5), iLv,,,, with transformation rule (E.2.9)"

v~, (At/8) - exp [-vL2 ( A t / 4 ) A t / 8 ] v~,. 3Thiscanbegeneralized,givingthe identity

exp( ~

~

e,p(o 0g(-I 0 )r a

[g-1 (y

540

Appendix E. Integration Schemes

Algorithm 30 (Equations of Motion: Nos4-Hoover)

subroutine

integrate

integrate equations of motion Nos6-Hoover thermostat

call chain(uk) cal i pos_vel (uk) call chain(uk) return end

Comments to this algorithm: 1. This subroutine solves the equations of motion for a single time step At using the Trotter equations (E.2.4)and (E.2.5). 2. In the subroutine c h a i n we apply e iLcat/4 to the current state (see Algorithm 31). 3. In the subroutine p o s _ v e l we apply e (iL~+iLv)At to the current state (see Algorithm 32). 4. uk is the total kinetic energy.

The next step is to apply iLG1, followed by again il_va,, etc. In this way we continue to apply all operators on the output of the previous step. Applying the Nos6-Hoover part of the Liouville operator changes 6k, v~.k, and vi. The other two Liouville operators change vi and r~. This makes it convenient to separate the algorithm into two parts in which the positions and velocities of the particles and the Nos6-Hoover chains are considered separately. An example of a possible implementation is shown in Algorithm 30.

E.2.2

The Isothermal-Isobaric Ensemble

Similar to the canonical ensemble we can derive a time-reversible integration scheme for simulation in the NPT ensemble. The equations of motions are given by expressions (6.2.1)-(6.2.8): t~

=

Pi + PE ~ri

V

=

dVpJw

541

E.2 Nos~-Hoover Algorithms

Algorithm 31 (Propagating the chain)

subroutine chain (uk)

G2=(QI*vxiI*vxiI-T) vxi2=vxi2+G2*delt4 vxil=vxil*exp(-vxi2*delt8) GI=(2*uk-L*T)/QI vxil=vxil+Gl*delt4 vxil=vxil*exp(-vxi2*delt8) xil=xil+vxil*delt2 xi2=xi2+vxi2*delt2 s=exp(-vxil*delt2) do i=l,npart v(i)=s*v(i) enddo uk=uk*s*s vxil=vxil*exp(-vxi2*delt8) GI=(2*uk-L*T)/QI vxil=vxil+Gl*delt4 vxil=vxil*exp(-vxi2*delt8) G2=(Ql*vxil*vxil-T)/Q2 vxi2=vxi2+G2*delt4 return end

apply equation (E.2.5) to the current position Update va2 using equation (E.2.6) Update va, using equation (E.2.7) Update va, using equation (E.2.8) Update va, using equation (E.2.7) Update t l using equation (E.2.9) Update &2 using equation (E.2.10) Scale factor in equation (E.2.11) update v~ using equation (E.2.11) update kinetic energy Update va, using equation (E.2.7) Update va 1 using equation (E.2.8) Update va, using equation (E.2.7) Update va2 using equation (E.2.6)

Comments to this algorithm: 1. In this subroutine T is the imposed temperature, d e l r d e l t 4 = At/4, and d e l t 8 = A t / 8 . 2. u k

is the total kinetic energy.

At, d e l r

= At/2,

542

Appendix E. Integration Schemes

A l g o r i t h m 32 (Propagating the Positions and Velocities)

subroutine

pos_vel (uk)

uk=O do i=l, npart x (i) =x (i) +v(i) *delt2 enddo call force do i=l, npart v(i) =v(i) +f (i) * d e l t / m x(i) =x(i)+v(i) *delt2 u k = u k + m * v (i ) *v (i ) / 2 enddo return end

apply equation (E.2.4) to the current position

update x~ using equation (E.2.13) calculate the force

update vi using equation (E.2.12) update xi using equation (E.2.13) update kinetic energy

Comments to this algorithm: 1. In this subroutine c l e l t - At and d e l t 2 -- At/2. 2. The subroutine f o r c e calculates the force on the particles.

1~~

-

~,k

=

19a,

-

l~ak

=

l~a~

=

1 ~ - p2 dV (Pint -- Pext) -4- ~- i=1 rrti Pak Qk

fork-1,...,M

p2 i=1

p~,,pr Q1

mi + ~

P~k-1 Qk-1

P~M-1 QM-1

PL2p

- - ( d N + 1 ) k B T - -~2 L,

-- kBT - Pa~+------~'pLk Qk+l

for k -- 2 , . . . , M - 1

kBT.

To derive a time-reversible numerical integration scheme to solve the equations of motion we use again the Liouville approach. A state is characterized by the variables 11 - (r TM, pN, e, pc, tM, p~A). The Liouville operator is defined by 0 iLNPT -- fl 0q"

E.2 Nos~-Hoover Algorithms

543

For a chain of length M - 2, using pi - mivt(# mill), Pak - Qkvak, r - (ln V)/d, and Pn - Wvn, the Liouville operator for these equations of motion can be written as iLNPT -- iL~ + iL~ + iLce, in which we define the operators N iLr

--

iLv

=

iLcp

0

E (Vi + V c r i ) 9 ~Tr i + Vr 0 r i=1

N rF /r/] Vv ~L

mi N

M

i=1

k=l

O

=

(| +

+GM 0v~. M0

M-1

Ovak

k--1

~

i:lL vcvi 9Vvr + (Gr - vcva, ) 0re

with

G1

--

I__Q[ ]L miv2 t + Wv2 = _ l(Nf + 1)kBT

Gk

=

Qk (Qk-,V~ , - k B T )

GC

1

1 w

1+

g

7- ~,q i=1

aU(r,V)

+ y - r ~ - F ~ ( r / ~ - d V ~0 V

- dPextV

]

i=l

An appropriate Trotter equation for the equations of motion is [85] e (iLNpTAt) = e ( i L c p A t / 2 ) e ( i L v A t / 2 ) e ( i L ~ A t ) e ( i L ~ A t / 2 ) e (ikcPAt/2) + O ( A t 3) .

(E.2.14) The operator iLcp has to be further factorized" iLce - iL L + ikcv + ikcc + iL~c + iLcl + iL~L, + iLc2,

9

544

Appendix E. Integration Schemes

where the terms are defined as

2

iL L =

0 5" Yak O~k

k=l iLc~

-

-f-[va,

+(l+~N)]Vi.Vv

i

i=l 0

iLcr

=

Gr162

iLve

=

- ( v a , vr

iLG1

=

G1 Ova,

iLva,

-

- ( v a , va2)

ovc

0

iLG2

=

ava,

0

G20v~2.

The Trotter expansion of the term iLc is

e(iLcpAt/2)

e(ikG 2At/4+iLv a, At/4) e(iLG, At/4) e(iLG r At/4+iLv e At/4) X e (iLaAt/2) e(ikc~ At/2) xe(iL~r162 , At/4)e(iLc 2At/4) e(iLc2At/4) [e(iL~,At/8)e(iLc,At/4)e(iLva, At/8)] X [e(iL~ At/8) e(iLc ~At/4) e(iL~ At/8) ] X e(iLa At/2)e(iLc~At/2)[e(iL~ ~At/8) e(iLc ~At/4) e(iLv ~At/8) ] X [e(ik~a, At/8)e(iLc'At/4)e(iL~a, At/8)] e(iCc2At/4).(E.2.15)

Similar to the NVT version the transformation rules of the various operators can be derived and translated into an algorithm. Such an algorithm is presented in ref. [85].