Multiple Time Step Methods and an Improved Potential Function for

13 Multiple Time Step Methods and an Improved Potential Function for Molecular Dynamics Simulations of Downloaded by UNIV OF CALIFORNIA SANTA BARBARA on November 24, 2015 | http://pubs.acs.org Publication Date: June 1, 1978 | doi: 10.1021/bk-1978-0086.ch013

Molecular Liquids W. B. STREETT and D. J. TILDESLEY 1

2

Science Research Laboratory, United States Military Academy, West Point,NY10996 G. SAVILLE Department of Chemical Engineering, Imperial College of Science, Prince Consort Road, London S.W.7., England Among the important recent advances i n computing methods for molecular dynamics s i m u l a t i o n s o f condensed matter are the s i n g u l a r i t y free a l g o r i t h m f o r rigid polyatomics, developed by Evans and Murad(1,2),and the m u l t i p l e time step (MTS) method, developed by S t r e e t t , et al. (3). In the method o f Evans and Murad, the equations o f r o t a t i o n a l motion are expressed i n terms o f four parameters, x, n, ξ, ζ, c a l l e d quaternions, def i n e d by X=Cos (6/2) Cos ((!/;+) /2) n=Sin(9/2)Cos((!ji-)/2) S=Sin(6/2)Sin((iJH>)/2) C=Cos(0/2)Sin((i|;+)/2) , where θ, Ψand ф are the E u l e r angles d e f i n e d by G o l d s t e i n (4). They report s i g n i f i c a n t improvements i n accuracy and e f f i c i e n c y over methods based on the equations o f E u l e r , Lagrange o r Ham ilton, using E u l e r angles. The MTS method of S t r e e t t , e t al., based on the use o f two o r more time steps o f d i f f e r e n t lengths to i n t e g r a t e the equations o f motion i n systems governed by continuous p o t e n t i a l f u n c t i o n s , has r e s u l t e d i n i n c r e a s e s i n computing speed by f a c t o r s o f three t o eight over conventional methods (3). The o r i g i n a l MTS paper d e s c r i b e s i n d e t a i l how the method i s a p p l i e d t o systems o f s p h e r i c a l molecules. In t h i s paper we give s i m i l a r d e t a i l s o f the a p p l i c a t i o n o f the MTS method t o systems o f rigid polyatomic molecules o f a r b i t r a r y shape, using the a l g o r i t h m o f Evans and Murad. We a l s o

Current address: Chemical Engineering Department, Cornell University, Ithaca, NY 14853. Current address: Physical Chemistry Laboratory, South Parks Road, Oxford OX1 3QZ, England. 1

2

0-8412-0463-2/78/47-086-144$05.00/0 © 1978 American Chemical Society In Computer Modeling of Matter; Lykos, P.; ACS Symposium Series; American Chemical Society: Washington, DC, 1978.

13.

Multiple Time Step Methods of Molecular Liquids 145

STREETT ET AL.

d e s c r i b e a modified s i t e - s i t e p o t e n t i a l , the s h i f t e d force po t e n t i a l , that gives b e t t e r accuracy and s t a b i l i t y than the usual truncated p o t e n t i a l .

Downloaded by UNIV OF CALIFORNIA SANTA BARBARA on November 24, 2015 | http://pubs.acs.org Publication Date: June 1, 1978 | doi: 10.1021/bk-1978-0086.ch013

S i t e - S i t e P o t e n t i a l s For Molecular

Liquids

Since s i t e - s i t e p o t e n t i a l s (also c a l l e d atom-atom poten t i a l s ) have been widely used i n computer s i m u l a t i o n and theo r e t i c a l s t u d i e s o f molecular l i q u i d s , we describe the a p p l i c a t i o n o f MTS methods to a model p o t e n t i a l of t h i s type. MTS methods are completely general, and can be a p p l i e d e q u a l l y w e l l to other types o f p o t e n t i a l f u n c t i o n s . In s i t e - s i t e models the p o t e n t i a l energy between a p a i r o f molecules i s the sum o f the p o t e n t i a l energies between p a i r s o f s i t e s on d i f f e r e n t molecules. These s i t e s are u s u a l l y atoms or groups of atoms, which i n t e r a c t with s i t e s on other molecules through s p h e r i c a l l y symmetric p o t e n t i a l f u n c t i o n s . Thus the p o t e n t i a l energy between mole cules 1 and 2 takes the form

u

(

r

12

h

)

i 2 ^ u

,

B

u

aB

(

r

aB ' )

where OK i s the set of angles d e f i n i n g the o r i e n t a t i o n o f mole cule i , g ( g ) Is the p o t e n t i a l energy between s i t e a on molecule 1 and s i t e 6 on molecule 2 (a f u n c t i o n only o f the s i t e - s i t e distance and the summation i s over a l l s i t e s i n each molecule. The methods d e s c r i b e d h e r e i n are a p p l i c a b l e to models governed by s i t e - s i t e p o t e n t i a l s that are continuous functions of r . Truncated P o t e n t i a l s . In the i n t e r e s t o f computing e f f i c iency, s i t e - s i t e p o t e n t i a l s are truncated at a c u t o f f d i s t a n c e , r , u s u a l l y taken to be about 2.5 times the e f f e c t i v e "diameter" o f the s i t e . I n t e r a c t i o n s beyond t h i s d i s t a n c e are neglected, and c o r r e c t i o n s to the c a l c u l a t e d thermodynamic p r o p e r t i e s , due to the neglected " t a i l " , are added using a simple p e r t u r b a t i o n scheme (_5) . E r r o r s inherent i n the use of truncated p o t e n t i a l s i n molecular dynamics simulations have u s u a l l y been overlooked, or at l e a s t neglected. These e r r o r s a r i s e because of the im pulse f e l t by two molecules whose s e p a r a t i o n r crosses the boundary r=r i n s u c c e s s i v e time s t e p s . There i s a gain or l o s s of p o t e n t i a l energy that i s not balanced by an equal and opposite change i n k i n e t i c energy; t o t a l energy i s not conser ved, and the c a l c u l a t e d t r a j e c t o r i e s depart s l i g h t l y from true Newtonian t r a j e c t o r i e s . (Numerical round o f f causes s i m i l a r errors.) In simulations of homogeneous f l u i d s there are, on average, as many c r o s s i n g s of the boundary i n one d i r e c t i o n as the other, so there i s no appreciable d r i f t i n the t o t a l energy and no s i g n i f i c a n t e r r o r i n the c a l c u l a t e d e q u i l i b r i u m p r o p e r t i e s ; however, i n simulations o f non-homogeneous systems (e.g., g a s - l i q u i d i n t e r f a c e ) these impulses can be a s i g n i f i u

r

a

a

r

c

In Computer Modeling of Matter; Lykos, P.; ACS Symposium Series; American Chemical Society: Washington, DC, 1978.

146

COMPUTER MODELING OF MATTER

cant source o f e r r o r and i n s t a b i l i t y . In a d d i t i o n , the accum u l a t i o n o f small e r r o r s i n p a r t i c l e t r a j e c t o r i e s i s l i k e l y t o r e s u l t i n s i g n i f i c a n t e r r o r s i n c a l c u l a t i o n s o f time-dependent phenomena that are based on c o r r e l a t i o n s o f p r o p e r t i e s a t long times, such as the "long time t a i l s " o f the v e l o c i t y autocor r e l a t i o n f u n c t i o n (6^) . For reasons t o be made c l e a r i n the next s e c t i o n , truncated p o t e n t i a l s have an adverse e f f e c t on energy c o n s e r v a t i o n i n some MTS methods. These problems can be overcome by simply i n c r e a s i n g the c u t o f f distance r t o a p o i n t where the i n t e r a c t i o n s omitted are t r u l y n e g l i g i b l e ; how ever, s i n c e the number o f p a i r i n t e r a c t i o n s c a l c u l a t e d i n the s i m u l a t i o n increases approximately as r^,, t h i s remedy i s ex pensive. As an a l t e r n a t i v e , we have developed a modified poten t i a l , c a l l e d the s h i f t e d force p o t e n t i a l , that provides accep t a b l e accuracy and energy conservation without increases i n r . The S h i f t e d Force P o t e n t i a l . When two s i t e s i n t e r a c t through a s p h e r i c a l l y symmetric p o t e n t i a l f u n c t i o n , they exert upon each other equal and opposite forces that a c t along the vector £ between them. (A wavy underline i n d i c a t e s a vector quantity.) The magnitude o f the force i s given by F=-(du/dr). In Figure 1(a) the p o t e n t i a l energy and force between two s i t e s are p l o t t e d against distance f o r a t y p i c a l truncated p o t e n t i a l . I t i s the d i s c o n t i n u i t y i n the force at r = r that i s the source of e r r o r i n molecular dynamics s i m u l a t i o n s . To reduce the e r r o r we use a m o d i f i e d f o r c e , F , formed by s u b t r a c t i n g from the usual e x p r e s s i o n f o r the force a constant, H, equal to the mag nitude o f t h i s d i s c o n t i n u i t y . The modified force i s


c

c

c

m

F =-(du/dr)-H,

rr.

where H=-(du/dr) The

r=r m o d i f i e d p o t e n t i a l i s defined

i n terms o f t h i s f o r c e : (2)

This

gives u(r)-Hr+G,

r i , s i n c e the i - j and j - i i n t e r a c t i o n s are i d e n t i c a l . ) Separate c o n t r i b u t i o n s t o the thermodynamic p r o p e r t i e s (pressure, energy, etc.) i n t o primary and secondary components, and c a l c u l a t e the time d e r i v a t i v e s o f the l a t t e r , using the expressions given i n reference 3. (b) Subsequent Time Step. C a l c u l a t e the primary force on each s i t e and the primary-neighbor c o n t r i b u t i o n s t o the thermodynamic p r o p e r t i e s by a summation over a l l p a i r s i n c l u d e d i n the l i s t o f primary neighbors compiled i n the f i r s t time step. C a l c u l a t e secondary f o r c e s from equation (3) and secondary-neighbor c o n t r i b u t i o n s t o the thermodynamic p r o p e r t i e s from s i m i l a r equations (see reference 3 ) . Add the primary and secondary components o f the c o n t r i b u t i o n s 2 + r

2

2

2

2

2

2

2

c

s

m

a

s


13.

Multiple Time Step Methods of Molecular Liquids

STREETT ET AL.

157

to the thermodynamic p r o p e r t i e s . Add the primary and secon dary forces f o r each s i t e , and c a l c u l a t e the exact f o r c e and torque, i n the s p a c e - f i x e d frame, on each molecule. Trans form s p a c e - f i x e d torques, T, i n t o torques i n the p r i n c i p a l axis system, T , using Jp=g~l.T, where E " i s the i n v e r s e o f the matrix i n equation (8). (The matrix i s orthogonal.) C a l c u l a t e the exact t r a n s l a t i o n a l and r o t a t i o n a l a c c e l e r a t i o n s , R' and w', from the d i f f e r e n t i a l equations o f motion: 1

p


1

R' '=F/M (A2) u)'=T / I p

y

y

y=x,y,z,

,

P

y

where M i s the molecular mass and Ip , Ip , I p the p r i n c i p a l moments o f i n e r t i a (4). C a l c u l a t e exact Values o f the f i r s t d e r i v a t i v e s o f the quaternions from equation (10). 4. C o r r e c t o r S e c t i o n . Correct the p r e d i c t e d values o f R and , and t h e i r d e r i v a t i v e s , and the p r e d i c t e d values o f the qua ternions (but not t h e i r d e r i v a t i v e s ) , according to z

l(j> =R(J) + -corr ~pred ( . =(J

» C

°

r

. .. ° ^+

r

P

=

r

e

A t )

j

j=0,I,---,*,

'

/ j 8

1

A. *Ao) ,

j = 0 , l , ( A 3 )

(At)Vj!

d

^corr ^pred

A • • AR' ' ~ ,

+

e

'

t

c

"

At/1! 1 1

1

where AR , Aa) and A£' are the d i f f e r e n c e s between the pre d i c t e d values from s e c t i o n 2 and the exact values from sec t i o n 3. The constants Aj f o r the p r e d i c t o r - c o r r e c t o r method can be found i n t a b l e 9.1 o f reference 7. 5. Return to P r e d i c t o r S e c t i o n . To use the l i n e a r e x t r a p o l a t i o n MTS method, the same pro cedure i s followed, except that i n the second time step o f each block o f n steps, the primary and secondary f o r c e s are again c a l c u l a t e d s e p a r a t e l y , using the l i s t o f primary neighbors com p i l e d i n the f i r s t time step, and values o f F from equation (4) are used i n a f i r s t order T a y l o r s e r i e s (equation 3) during the next n-2 time steps. I f the torques a r e separated i n t o p r i mary and secondary components, and t r e a t e d i n the same manner as the f o r c e s , equations (3) and (4) can be a p p l i e d t o the t o t a l f o r c e on each molecule, r a t h e r than t o the force on each s i t e , r e s u l t i n g i n a s u b s t a n t i a l savings i n storage. (A s i m i l a r savings i n storage can be r e a l i z e d i n the T a y l o r s e r i e s MTS methods i f the torques are separated i n t o primary and secondary s


158

COMPUTER MODELING OF MATTER

components, T and T , and d e r i v a t i v e s of T c a l c u l a t e d from d e r i v a t i v e s o f equation (5a); however, the savings i n storage w i l l be o f f s e t by the i n c r e a s e i n computing time r e q u i r e d to evaluate these d e r i v a t i v e s . In t h i s case i t i s more e f f i c i e n t to sum the forces on the s i t e s i n a double loop over a l l r e l e vant molecular p a i r s , and then to c a l c u l a t e the torques from equation (5a) i n a s i n g l e loop over a l l molecules.) g

g


Literature Cited 1. 2. 3.

Evans, D. J . , Molec. Phys. (1977) 34, 317. Evans, D. J . & Murad, S., Molec. Phys. (1977) 34, 327. S t r e e t t , W. B., T i l d e s l e y , D. J . and S a v i l l e , G., Molec. Phys. (1978) in p r e s s . 4. G o l d s t e i n , H., " C l a s s i c a l Mechanics", Chaps. 4 & 5, AddisonWesley, Reading, Massachusetts, 1971. 5. McDonald, I. R. & Singer, K., D i s c . Faraday Soc. (1967) 43, 40. 6. A l d e r , B. J . & Wainwright, T. E., Phys. Rev. (1970) Al, 18. 7. Gear, C. W., "Numerical Initial Value Problems i n Ordinary Differential Equations," pp. 148-54, P r e n t i c e - H a l l , Englewood Cliffs, New J e r s e y , 1971.

RECEIVED September 7, 1978.


Multiple Time Step Methods and an Improved Potential Function for

Recommend Documents