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.