Molecular Dynamics Simulations of Simple Fluids with Three-Body

with long range, three body interactions explicitly included. The technique reported ... energy (12), and the radial distribution function (10, _13, 1...
0 downloads 0 Views 1MB Size
15 Molecular Dynamics Simulations of Simple Fluids with

Downloaded by CALIFORNIA INST OF TECHNOLOGY on November 27, 2017 | http://pubs.acs.org Publication Date: June 1, 1978 | doi: 10.1021/bk-1978-0086.ch015

Three-Body Interactions Included J. M. HAILE Chemical Engineering Department, Clemson University, Clemson, SC 29631

A primary goal i n the study of the fluid s t a t e s of matter i s an ability to p r e d i c t all static and dynamic p r o p e r t i e s using only information about the c o n s t i t u e n t molecules. The f l u i d systems of i n t e r e s t here are those i n which the temperature and molecular masses are sufficiently high that quantum e f f e c t s are n e g l i g i b l e and f o r which the classical Hamiltonian can be assumed separable i n t o independent t r a n s l a t i o n a l , r o t a t i o n a l , i n t e r n a l , and c o n f i g u r a t i o n a l c o n t r i b u t i o n s . The p r i n c i p a l remaining diffi­ culty i n r e a l i z i n g the goal of complete property p r e d i c t i o n f o r such systems i s q u a n t i t a t i v e d e s c r i p t i o n of the i n t e r m o l e c u l a r f o r c e s w i t h i n the fluid. These f o r c e s are commonly discussed i n terms of an i n t e r m o l e c u l a r p o t e n t i a l energy U which, f o r a f l u i d of N s p h e r i c a l molecules, depends only on the p o s i t i o n s of the molecular mass centers r : i

— l

U = U(r

r

r , . . . r^)

(1)

2

This i n t e r m o l e c u l a r p o t e n t i a l i s u s u a l l y w r i t t e n as a s e r i e s of terms which i n d i v i d u a l l y account f o r 2-body, 3-body,... N-body interactions: U - 11 i r i j * ik' jk L c # 1

T

(11)

and the t r u n c a t i o n o f eq. (6) a p p l i e s . T h i s a l g o r i t h m i n which a neighbor l i s t i s used f o r both two and three body f o r c e e v a l u ­ a t i o n and i n which three body f o r c e s are c a l c u l a t e d e x p l i c i t l y at each time step s h a l l be r e f e r r e d to as the c o n v e n t i o n a l mole­ c u l a r dynamics (CMD) method. U n f o r t u n a t e l y , even with the neighbor l i s t a p p l i e d to the three body f o r c e s , the CMD program f o r 108 Lennard-Jones plus A x i l r o d - T e l l e r p a r t i c l e s executes about 15 times slower than the same program f o r 108 Lennard-Jones p a r t i c l e s . In the work r e ­ ported here, source programs were w r i t t e n i n FORTRAN IV u s i n g s i n g l e p r e c i s i o n a r i t h m e t i c , compiled on an IBM FORTRAN IV-G compiler, and executed on the IBM 370/165 a t Clemson U n i v e r s i t y . 3.2 M u l t i p l e Time Step (MTS) Method A p p l i e d to Three Body Forces. In order to improve the execution speed of simula­ tion programs with three body i n t e r a c t i o n s i n c l u d e d , the mul­ t i p l e time step method of S t r e e t t , e t a l . (5) has been a p p l i e d to the e v a l u a t i o n of the three body f o r c e s . The m u l t i p l e time step method attempts to take advantage of the f a c t that molecular motions i n a f l u i d may be reduced to components which operate on very d i f f e r e n t time s c a l e s . More p r e c i s e l y , one can o f t e n i d e n ­ t i f y components o f the f o r c e on a molecule which have r e l a t i v e l y l a r g e d i f f e r e n c e s i n t h e i r r a t e s of change with time. I t i s the q u i c k l y v a r y i n g component of the f o r c e which l i m i t s the s i z e of the time step At which must be used to o b t a i n s t a b l e s o l u t i o n s

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

Downloaded by CALIFORNIA INST OF TECHNOLOGY on November 27, 2017 | http://pubs.acs.org Publication Date: June 1, 1978 | doi: 10.1021/bk-1978-0086.ch015

180

COMPUTER

MODELING OF MATTER

to the d i f f e r e n t i a l equations of motion. In the MTS method these q u i c k l y v a r y i n g components of the f o r c e are e x p l i c i t l y evaluated at each time step i n the u s u a l way; however, the slowly v a r y i n g components are only e x p l i c i t l y evaluated at longer i n t e r v a l s , say, every n time steps. At time steps between e x p l i ­ c i t determination, the slowly v a r y i n g components of the f o r c e are estimated by e x t r a p o l a t i o n from the previous e x p l i c i t e v a l u ­ ation. The e x t r a p o l a t i o n may be l i n e a r (LMTS method) or, i f one i s more ambitious, based on a T a y l o r s e r i e s expansion. The degree to which the MTS method improves execution speed of the program depends on the amount of computation avoided i n e x t r a ­ p o l a t i n g the slowly v a r y i n g f o r c e components r a t h e r than e x p l i ­ c i t l y c a l c u l a t i n g them. I f we consider a p p l y i n g t h i s MTS method to s i m u l a t i o n s with three body i n t e r a c t i o n s , eq. (8) already represents the f o r c e on molecule i d i v i d e d i n t o two p a r t s . Barker and Henderson have noted that the three body A x i l r o d - T e l l e r i n t e r a c t i o n i s slowly v a r y i n g (21). To t e s t t h i s , we performed a molecular dynamics s i m u l a t i o n f o r 108 p a r t i c l e s i n t e r a c t i n g with Lennard-Jones plus A x i l r o d - T e l l e r p o t e n t i a l s u s i n g the conventional molecular dyna­ mics algorithm. Figure 3 shows a t y p i c a l comparison of a compo­ nent of the two and three body f o r c e s on one p a r t i c l e during a segment of the s i m u l a t i o n . The f i g u r e i n d i c a t e s that the three body component of the f o r c e does indeed have a slower r a t e of change than the two body f o r c e component. Encouraged by F i g u r e 3, the MTS method was then a p p l i e d to the s i m u l a t i o n i n the f o l l o w i n g manner. The e n t i r e three b o d y ^ f o r c e i s considered to be slowly v a r y i n g , t h e r e f o r e a l l of _F. i s e x p l i c i t l y evaluated only at p e r i o d i c i n t e r v a l s i n ^ t h e simu­ lation. We choose to use l i n e a r e x t r a p o l a t i o n of F. for, although S t r e e t t , et a l . f i n d that a t h i r d order T a y l o r s e r i e s e x t r a p o l a t i o n i s more accurate, the a n a l y t i c e v a l u a t i o n s of the time d e r i v a t i v e s o^ the A x i l r o d - T e l l e r f o r c e are h o p e l e s s l y tedious. Thus, F\ i s e x p l i c i t l y evaluated at two s u c c e s s i v e times, t and (t + A t ) , and then l i n e a r l y e x t r a p o l a t e d over the next n time steps. T h i s gives the three body f o r c e on p a r t i c l e i at any intermediate time step ( t + k At) as: fi

Q

Q

F. — i

3 B

(t + k At) = F . o —1

3 B

3 B

( t ) + k [ F . ( t + At) - F . o — i o —x

3 B

( t )] o

(12)

In a s i m u l a t i o n program using a p r e d i c t o r - c o r r e c t o r a l g o ­ rithm, eq. (12) would appear i n the f o r c e e v a l u a t i o n step. The only a d d i t i o n a l storage r e q u i r e d f o r the MTS method over that f o r the CMD method i s ^ s p a c e f o r 6N numbers - 3N l o c a t i o n s f o r the components of _F. ( t ^ ) * 3N l o c a t i o n s f o r the components of the d i f f e r e n c e term i n brackets i n eq. (12). We experimented with v a r i o u s values of the time step At and number of e x t r a p o l a t i o n steps n. A compromise among s t a b l e s o l u t i o n s of the d i f f e r e n t i a l equations, conservation of system energy, and execution speed was obtained using n = 8 and a n c

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

Downloaded by CALIFORNIA INST OF TECHNOLOGY on November 27, 2017 | http://pubs.acs.org Publication Date: June 1, 1978 | doi: 10.1021/bk-1978-0086.ch015

15.

HAILE

MD Simulations of Fluids with Three-Body Interactions

181

At = 2.15 (10 ) s e c . T h i s value of the time step i s about onef i f t h that used by V e r l e t (31) and Rahman (34) i n s i m u l a t i o n s of the p a i r w i s e a d d i t i v e Lennard-Jones f l u i d . Note t h a t , i n the MTS method, any property which i s being c a l c u l a t e d as a running ensemble average a t each time step during the s i m u l a t i o n must a l s o be d i v i d e d i n t o slowly (three body) and q u i c k l y (two body) v a r y i n g c o n t r i b u t i o n s . The MTS method i s a p p l i e d to the slowly v a r y i n g components of the p r o p e r t i e s v i a equations analogous to eq. (12). T h i s comment a p p l i e s , f o r example, to the i n t e r n a l energy and pressure i n the work reported here. Two p a r a l l e l s i m u l a t i o n s f o r 108 p a r t i c l e s i n t e r a c t i n g with Lennard-Jones p l u s A x i l r o d - T e l l e r p o t e n t i a l s have been performed. The f i r s t c a l c u l a t i o n u t i l i z e d the CMD method i n which the f o r c e s were e x p l i c i t l y evaluated a t each time step. In the second run the two body f o r c e s were determined i n the standard way and the LMTS method d e s c r i b e d above was a p p l i e d to the three body f o r c e s . Both runs were s t a r t e d from the same i n i t i a l par­ t i c l e p o s i t i o n s and v e l o c i t i e s and both were continued f o r 1650 time steps. A comparison of the p r o p e r t i e s obtained from the two c a l c u l a t i o n s i s given i n Table I . In a d d i t i o n to the pro­ p e r t i e s l i s t e d i n Table I , r a d i a l d i s t r i b u t i o n f u n c t i o n s , v e l o ­ c i t y , speed, and f o r c e a u t o c o r r e l a t i o n f u n c t i o n s , and atomic mean squared displacements (from which d i f f u s i o n c o e f f i c i e n t s were obtained) were c a l c u l a t e d . For a l l of these p r o p e r t i e s , the LMTS values were w i t h i n 0.1% of the values obtained by the CMD method. Figure 4 shows the per cent d e v i a t i o n i n the i n s t a n ­ taneous t o t a l energy o f the two c a l c u l a t i o n s . Since, i n the LMTS method d e s c r i b e d above, the three body f o r c e s a r e e x p l i c i t l y evaluated only two of every ten time s t e p s , we might expect the LMTS method to be about 5 times as f a s t as the CMD method. The r e s u l t s i n Table I confirm t h i s ; i . e . , u s i n g the LMTS method, a s i m u l a t i o n of 108 Lennard-Jones p l u s A x i l r o d - T e l l e r p a r t i c l e s r e q u i r e s about 20 CPU minutes f o r 2000 time steps on an IBM 370/165. T h i s execution speed remains about 2.5 times slower than a CMD method, p a i r w i s e a d d i t i v e LennardJones s i m u l a t i o n . The r e s u l t s from the LMTS method and the CMD method a r e i n good agreement f o r those p r o p e r t i e s which a r e averaged over the s i m u l a t i o n . However, a f t e r about 1000 time steps, the i n s t a n ­ taneous components of the three body f o r c e s on i n d i v i d u a l par­ t i c l e s i n the LMTS c a l c u l a t i o n begin to d e v i a t e from those i n the CMD method, as i n d i c a t e d i n Table I . These d e v i a t i o n s a r e due to gradual accumulation o f i n a c c u r a c i e s from the l i n e a r e x t r a p o l a t i o n of the three body f o r c e . I t should be emphasized that while these d e v i a t i o n s occur, the t o t a l energy and momentum of the LMTS system remain c o n s e r v a t i v e and agree c l o s e l y with those of the CMD method (see F i g u r e 4 ) . This disagreement i n the three body f o r c e s can be i n t e r p r e t e d as a d r i f t i n g o f the

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

COMPUTER MODELING OF MATTER

Downloaded by CALIFORNIA INST OF TECHNOLOGY on November 27, 2017 | http://pubs.acs.org Publication Date: June 1, 1978 | doi: 10.1021/bk-1978-0086.ch015

182

Figure 3. Comparison of z-component of two-(2B) and three-body (3B) forces on particle #27 in conventional molecular dynamics simulation of Lennard-Jones plus Axilrod-Teller fluid. At = time step = 2.15 (10' ) sec, a = 0.65, kTV = 1.033. Note that the scale for three-body force is an order of magnitude smaller than that for the two-body force. 15

P

1

3

%

-0.1 L

1 1400

1

1 1600

I

I 1800

I

1 2000

At

Figure 4. Percent deviation in instantaneous total system energy calculated by MTS method from that determined by CMD method for 108 particles interacting with Lennard-Jones plus Axilrod-Teller forces. System parameters same as those in Figure 3. Percent deviation calculated as: (E MD — EUTS)100/E UD' C

C

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

15. HAILE

MD Simulations of Fluids with Three-Body Interactions

183

Table I Test of the L i n e a r M u l t i p l e Time Step Method A p p l i e d to the A x i l r o d - T e l l e r P o t e n t i a l i n Molecular Dynamics Simulations of Lennard-Jones p l u s A x i l r o d - T e l l e r I n t e r a c t i o n s .

Downloaded by CALIFORNIA INST OF TECHNOLOGY on November 27, 2017 | http://pubs.acs.org Publication Date: June 1, 1978 | doi: 10.1021/bk-1978-0086.ch015

CMD Method

Number P a r t i c l e s , N

108

LMTS Method

108

3 0.65

0.65

2.15

2.15

Number d e n s i t y , pa Time step, At (10"^)sec 71.2

17.1

CPU time, mins. 1.0325 Temperature,

1.0326

kTe ^ -4.356

-4.356

0.104

0.103

-2.795

-2.797

3B 3 -1 a e

-0.363

-0.373

3B 3 -1 a e

0.125

0.120

3B 3 -1 ae

0.0465

0.0519

xa

1.457

1.462

v e c t o r f o r p a r t i c l e 27

ya

4.981

4.987

at time step 1650

za

1.358

1.364

Average Conf. I n t . Energy, U(Ne) Average Pressure, P(pkT) ^ Instantaneous T o t a l Energy E(Ne) at time step 1650 Instantaneous Components of 3-body Force

p

on p a r t i c l e 27 a t time

r

F

x

F y

step 1650

Components of p o s i t i o n

r

F

z

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

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

P(pkT)

-1

-1

-1

U(Nc)

kTe

Number Time Steps

Time Step ( 1 0 ) sec

1 5

Number P a r t i c l e s

3

= 0.65

-0.11

-4.52

1.036

1200

9.6

864

-0.05

-4.54

1.052

2000

2.15

108

Verlet This ( r e f . (31)) Work

Lennard-Jones

pa

+0.10

-4.36

1.033

1650

2.15

108

MTS Method

L J + AT

pa

-0.20

-5.90

+0.38

-5.64

0.746

1730

2.15

108

MTS Method

LJ + AT

= 0.817

0.746

2000

2.15

108

This Work

Lennard-Jones

P r e l i m i n a r y R e s u l t s Showing E f f e c t of Three Body A x i l r o d - T e l l e r Force on I n t e r n a l Energy and Pressure

Table I I

Downloaded by CALIFORNIA INST OF TECHNOLOGY on November 27, 2017 | http://pubs.acs.org Publication Date: June 1, 1978 | doi: 10.1021/bk-1978-0086.ch015

Downloaded by CALIFORNIA INST OF TECHNOLOGY on November 27, 2017 | http://pubs.acs.org Publication Date: June 1, 1978 | doi: 10.1021/bk-1978-0086.ch015

15.

HAILE

MD Simulations of Fluids with Three-Body Interactions

185

LMTS system onto a d i f f e r e n t t r a j e c t o r y i n phase space from that followed by the system evolved by the CMD method. Hence, there i s no s t r i c t l y Newtonian connection between widely separated phase p o i n t s ; but then, there i s no true Newtonian connection between widely separated phase p o i n t s i n the system generated by the CMD method due to round o f f e r r o r s and the truncated poten­ tial. The c o n c l u s i o n i s that the phase space t r a j e c t o r i e s become d i f f e r e n t i n the two c a l c u l a t i o n s , but one i s not n e c e s s a r i l y more wrong (or r i g h t ) than the other. Hoover and Ashurst have discussed t h i s p o i n t i n another context with s i m i l a r c o n c l u s i o n s (35). The r e s t r i c t i o n which t h i s imposes i s that the LMTS method cannot be used to study very l o n g - l i v e d phenomena; such as the long t a i l of the v e l o c i t y a u t o c o r r e l a t i o n f u n c t i o n . 4.

Preliminary

Results

for Equilibrium

Properties

The l i n e a r m u l t i p l e time step method described i n S e c t i o n 3 has been used to simulate 108 p a r t i c l e s i n t e r a c t i n g w i t h LennardJones {jlus A x i l r o d - T e l l e r p o t e n t i a l s at two s t a t e c o n d i t i o n s : (a) p a =0.65, kTe~ = 1.036, which i s one of the c o n d i t i o n s of the Lennard-Jones f l u i d s t u d i e d by V e r l e t (31) and by Singh, et a l . (14, 15), (b) p a = 0.817, kTe = 0.746, which i s c l o s e to the c o n d i t i o n of the Lennard-Jones f l u i d s t u d i e d by Rahman (34). From these two c a l c u l a t i o n s the c o n f i g u r a t i o n a l i n t e r n a l energy, pressure, and r a d i a l d i s t r i b u t i o n f u n c t i o n have been determined. For comparison, s i m u l a t i o n runs f o r a p u r e l y Lennard-Jones f l u i d have a l s o been performed at s t a t e c o n d i t i o n s c l o s e to the above. E q u i l i b r i u m property r e s u l t s at both s t a t e c o n d i t i o n s f o r the Lennard-Jones and Lennard-Jones plus A x i l r o d - T e l l e r f l u i d s are compared i n Table I I . Long range c o r r e c t i o n s f o r the truncated two body p o t e n t i a l have been added to the i n t e r n a l energy and pressure; however, long range c o r r e c t i o n s f o r the truncated three body p o t e n t i a l have been neglected i n the v a l u e s given i n Table I I . These r e s u l t s confirm the e a r l i e r work by Barker, et a l . i n that the three body e f f e c t i s only a few percent on the i n t e r n a l energy, but i s s i g n i f i c a n t l y l a r g e r on the pressure. Note t h a t , at the s t a t e c o n d i t i o n s s t u d i e d , the net e f f e c t of the A x i l r o d - T e l l e r i n t e r a c t i o n i s r e p u l s i v e , as evidenced by p o s i t i v e c o n t r i b u t i o n s to the energy and pressure. In Figure 5 the r a d i a l d i s t r i b u t i o n f u n c t i o n s f o r the LennardJones and Lennard-Jones plus A x i l r o d - T e l l e r f l u i d s are compared at the high d e n s i t y s t a t e c o n d i t i o n . I t appears t h a t the three body i n t e r a c t i o n s have no e f f e c t except i n the f i r s t peak r e g i o n where the t r i p l e t i n t e r a c t i o n s lower g(r) by about 3%. However, we estimate the s t a t i s t i c a l e r r o r i n g(r) f o r these s i m u l a t i o n s i s about +3%, so i t i s unclear to what extent the observed effect i s real. We note that Schommers has reported s m a l l three body e f f e c t s on the d i s t r i b u t i o n f u n c t i o n f o r the two dimensional f l u i d , as w e l l (16). However, the p e r t u r b a t i o n theory c a l c u l a 3

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

Downloaded by CALIFORNIA INST OF TECHNOLOGY on November 27, 2017 | http://pubs.acs.org Publication Date: June 1, 1978 | doi: 10.1021/bk-1978-0086.ch015

COMPUTER MODELING OF MATTER

Figure 5. Effect of Axilrod-Teller potential on radial distribution function at pa = 0.817. For Lennard-Jones fluid, kTY = 0.746; for Lennard-Jones plus Axilrod-Teller fluid, kTY = 0.740. 3

1

1

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

Downloaded by CALIFORNIA INST OF TECHNOLOGY on November 27, 2017 | http://pubs.acs.org Publication Date: June 1, 1978 | doi: 10.1021/bk-1978-0086.ch015

15. HAILE

MD Simulations of Fluids with Three-Body Interactions

Figure 6. Comparison of trajectories of particle #27 in Lennard-Jones and Lennard-Jones plus Axilrod-Teller fluids at pa = 0.65. Both runs were started from the same point in phase space and trajectories shown are from time steps 5002000 in each simulation. For this calculation the Axilrod-Teller strength constant v was assigned a value of three times that for argon given in Section 2. Note that the circles represent positions of the center of mass of the atom, not the atomic diameter. 3

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

187

188

COMPUTER

MODELING OF MATTER

Downloaded by CALIFORNIA INST OF TECHNOLOGY on November 27, 2017 | http://pubs.acs.org Publication Date: June 1, 1978 | doi: 10.1021/bk-1978-0086.ch015

1 1

t i o n s of Singh and Ram show " s i g n i f i c a n t e f f e c t of the t r i p l e t i n t e r a c t i o n s on g(r) at low temperatures and high d e n s i t i e s (14). The g(r) f o r the low d e n s i t y c a l c u l a t i o n shows e s s e n t i a l l y no e f f e c t due to the three body f o r c e s . Singh and Ram a l s o f i n d l i t t l e e f f e c t at low d e n s i t i e s . Although some dynamic p r o p e r t i e s have been c a l c u l a t e d from these s i m u l a t i o n s , we f i n d , i n g e n e r a l , that 2000 time steps f o r 108 p a r t i c l e s do not provide s u f f i c i e n t data to o b t a i n s t a ­ t i s t i c a l l y meaningful dynamic p r o p e r t i e s . Work i s now i n pro­ gress aimed at i n c r e a s i n g the number of p a r t i c l e s i n the system and the l e n g t h of the run to improve the s t a t i s t i c a l p r e c i s i o n of dynamic property c a l c u l a t i o n s . In an e f f o r t to g a i n f u r t h e r enlightenment as to the e f f e c t s of the A x i l r o d - T e l l e r p o t e n t i a l , we have compared i n d i v i d u a l p a r t i c l e t r a j e c t o r i e s f o r the same p a r t i c l e over the same time segment from two d i f f e r e n t s i m u l a t i o n s . In one s i m u l a t i o n the i n t e r m o l e c u l a r p o t e n t i a l was Lennard-Jones while i n the other, the p o t e n t i a l was Lennard-Jones plus A x i l r o d - T e l l e r . Both c a l ­ c u l a t i o n s were done at the same d e n s i t y and were s t a r t e d from the same c o n f i g u r a t i o n and p a r t i c l e v e l o c i t i e s . F i g u r e 6 com­ pares the t r a j e c t o r i e s obtained from such a study. In the Lennard-Jones f l u i d the p a r t i c l e shown i n F i g u r e 6 i s e x h i b i t i n g v i b r a t o r y motion due to c o l l i s i o n s with neighboring p a r t i c l e s . In the Lennard-Jones plus A x i l r o d - T e l l e r f l u i d , the same p a r t i c l e during the same time segment i s undergoing an extended t r a j e c t o r y of l a r g e l y d i f f u s i o n a l motion. I t may be that study of three body f o r c e s on dynamic p r o p e r t i e s w i l l r e v e a l f a i r l y strong e f f e c t s on i n d i v i d u a l p a r t i c l e motion but that these e f f e c t s tend to cancel when averaged to o b t a i n bulk dynamic p r o p e r t i e s . Hence, Schommers two dimensional s i m u l a t i o n s show strong three body e f f e c t s on the v e l o c i t y a u t o c o r r e l a t i o n f u n c t i o n but i n s i g n i f i c a n t e f f e c t on the s e l f d i f f u s i o n c o e f f i c i e n t (16). 1

5.

Conclusions

A m u l t i p l e time step method has been s u c c e s s f u l l y developed f o r performing molecular dynamics simulations i n which three body i n t e r a c t i o n s are e x p l i c i t l y i n c l u d e d . The method has been t e s t e d f o r 108 p a r t i c l e s i n t e r a c t i n g with a Lennard-Jones plus A x i l r o d - T e l l e r p o t e n t i a l . When compared with a conventional molecular dynamics program u s i n g the same p o t e n t i a l model, the MTS program executes - 4 . 5 times f a s t e r while g i v i n g e s s e n t i a l l y the same values f o r e q u i l i b r i u m and dynamic p r o p e r t i e s , f o r runs of up to 2000 time steps. At t h i s p o i n t , the MTS method i s of d o u b t f u l r e l i a b i l i t y f o r studying l o n g - l i v e d phenomena. Based on t h i s study, i t seems that the m u l t i p l e time step methods could provide a mechanism f o r extending molecular dynamics simu­ l a t i o n s to a v a r i e t y of phenomena which i n v o l v e components that i n h e r e n t l y operate on d i f f e r e n t time s c a l e s ; f o r example, study of phase t r a n s i t i o n s or long-chain molecules.

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

Downloaded by CALIFORNIA INST OF TECHNOLOGY on November 27, 2017 | http://pubs.acs.org Publication Date: June 1, 1978 | doi: 10.1021/bk-1978-0086.ch015

15.

HAILE

MD Simulations of Fluids with Three-Body Interactions

189

In p r e l i m i n a r y s t u d i e s of three body e f f e c t s w i t h the MTS program, we have confirmed the work of Barker, et_ a l . (11, 12) that the A x i l r o d - T e l l e r c o n t r i b u t i o n i s only a few per cent f o r the i n t e r n a l energy but s i g n i f i c a n t l y l a r g e r f o r the p r e s s u r e . The e f f e c t of the A x i l r o d - T e l l e r i n t e r a c t i o n on the r a d i a l d i s t i r b u t i o n f u n c t i o n i s found to be s m a l l . This r e s u l t i s i n agreement w i t h Schommers' two dimensional s i m u l a t i o n s (16); but disagrees w i t h the t h e o r e t i c a l c a l c u l a t i o n s of Singh, et_ a l . at high d e n s i t y (14, 15). F i n a l l y , we note that the work reported here begins a study of long range, three body i n t e r a c t i o n s but does not i n c l u d e p o s s i b l e short range, three body e f f e c t s . P. A. E g e l s t a f f i s studying such e f f e c t s u s i n g neutron s c a t t e r i n g and Monte C a r l o s i m u l a t i o n (36). Acknowledgments The author i s g r a t e f u l t o : W. B. S t r e e t t and D. J . T i l d e s l e y f o r v a l u a b l e d i s c u s s i o n s on the m u l t i p l e time step method; H. W. Graben f o r d i s c u s s i o n s on three body f o r c e s ; P. A. E g e l s t a f f and J . A. Barker f o r i n s t r u c t i v e correspondence. T h i s work was supported, i n p a r t , by a grant from the F a c u l t y Research Committee, Clemson U n i v e r s i t y . The Clemson U n i v e r s i t y Computer Center generously provided the computer time used i n t h i s work.

Literature Cited 1. Hansen, J . P. and McDonald, I . R., "Theory of Simple L i q u i d s , " Academic P r e s s , New York, 1976. 2. Watts, R. 0. and McGee, I. J., " L i q u i d State Chemical P h y s i c s , " J . Wiley and Sons, New York, 1976. 3. M e t r o p o l i s , M., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. N., and T e l l e r , E., J. Chem. Phys. (1953) 21, 1087. 4. A l d e r , B. J . and Wainwright, T. E., J. Chem. Phys. (1957) 27, 1208. 5. S t r e e t t , W. B., T i l d e s l e y , D. J., and Saville, G., Molec. Phys. (1978) in p r e s s . 6. Copeland, D. A. and Kestner,N.R.,J.Chem. Phys. (1968) 49, 5214. 7. Casanova, G., D u l l a , R. J., Jonah, D. A., Rowlinson, J . S., and S a v i l l e , G., Molec. Phys. (1970) 18, 589. 8. D u l l a , R. J., Rowlinson, J. S., and Smith, W. R., Molec. Phys. (1971) 21, 299. 9. Sherwood, A. E. and P r a u s n i t z , J .M.,J.Chem. Phys. (1964) 41, 413. 10. Fowler, R. and Graben, H. W., J. Chem.Phys. (1972) 56, 1917. 11. Barker, J . A., F i s h e r , R. A., and Watts, R. O., Molec.Phys. (1971) 21, 657.

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

190

COMPUTER MODELING OF MATTER

Downloaded by CALIFORNIA INST OF TECHNOLOGY on November 27, 2017 | http://pubs.acs.org Publication Date: June 1, 1978 | doi: 10.1021/bk-1978-0086.ch015

12.

Barker, J . A., Henderson, D., and Smith, W. R., Molec. Phys. (1969) 17, 579. 13. Y a r n e l l , J . L., Katz, M. J . , Wenzel, R. G., and Koenig, S. H., Phys. Rev. A (1973) 7, 2130. 14. Ram, J . and Singh, Y., J. Chem. Phys. (1977) 66, 924. 15. Sinha, S. K., Ram, J . , and Singh, Y., J . Chem. Phys. (1977) 66, 5013. 16. Schommers, W., Phys. Rev. A (1977) 16, 327. 17. B e l l , R. J., J. Phys. B (1970) 3, 751. 18. A x i l r o d , B. M. and Teller, E., J. Chem. Phys.(1943) 11, 299. 19. Muto, Y., Proc. Phys. Math. Soc. Japan (1943) 17, 629. 20. Jansen, L. and Lombardi, E., Faraday D i s c . Chem. Soc. (London) (1965) 40, 78. 21. Barker, J . A. and Henderson, D., Rev. Mod. Phys. (1976) 48, 587. 22. V e r l e t , L., Phys. Rev. (1968) 165, 201. 23. F i s h e r , R. A. and Watts, R. O., A u s t r . J . Phys. (1972) 25, 529. 24. Singh, Y., Molec. Phys. (1975) 29, 155. 25. Shukla, K. P., Ram, J . , and Singh, Y., Molec. Phys. (1976) 31, 873. 26. Stogryn, D. E., J . Chem. Phys. (1970) 52, 3671. 27. Gear, C. W., "Numerical Initial Value Problems in Ordinary Differential Equations," Prentice-Hall, Englewood Cliffs, New Jersey, 1971. 28. Cheung, P. S. Y. and Powles, J . G., Molec. Phys. (1974) 30, 921. 29. H a i l e , J . M., "Surface Tension and Computer Simulation of Polyatomic F l u i d s , " Ph.D. D i s s e r t a t i o n , U n i v e r s i t y of F l o r i d a , G a i n e s v i l l e , 1976. 30. Wood, W. W., i n "Physics o f Simple L i q u i d s , " H. N. V. Temperley, J . S. Rowlinson, and G. S. Rushbrooke, eds, North-Holland, Amsterdam, 1968. 31. V e r l e t , L., Phys. Rev. (1967) 159, 98. 32. S c h o f i e l d , P., Comput. Phys. Comm. (1973) 5, 17. 33. Quentrec, B. and Brot, C., J. Comput. Phys. (1973) 13, 430. 34. Rahman, A., Phys. Rev. (1964) 136, 405. 35. Hoover, W. G. and Ashurst, W. T. in " T h e o r e t i c a l Chemistry," vol. 1, H. E y r i n g and D. Henderson, eds., Academic Press, New York, 1975. 36. E g e l s t a f f , P. A., p r i v a t e communication (1978). RECEIVED September 7, 1978.

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