Computer Modeling of Matter - American Chemical Society

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.


Computer Modeling of Matter - American Chemical Society

Recommend Documents