Molecular Dynamics Simulation of Methane Using a Singularity-Free

5 Molecular Dynamics Simulation of Methane Using a Singularity-Free Algorithm S. MURAD and K. E. GUBBINS

Downloaded by EMORY UNIV on March 8, 2016 | http://pubs.acs.org Publication Date: June 1, 1978 | doi: 10.1021/bk-1978-0086.ch005

School of Chemical Engineering, Cornell University, Ithaca, NY 14853

Rather extensive computer s i m u l a t i o n s t u d i e s have been r e ported f o r l i q u i d s o f linear molecules over the past few years (1). Much l e s s work has been reported f o r l i q u i d s of n o n l i n e a r mole c u l e s , water (2) being the only exception, although some simula t i o n work f o r ammonia (3), benzene (4,5) and n-butane (6) has been published r e c e n t l y . Simulations o f rigid polyatomics are complicated by the a d d i t i o n a l o r i e n t a t i o n v a r i a b l e s , and a l s o because the equations that must be solved f o r the r o t a t i o n a l motion possess a s i n g u l a r i t y at θ = 0 if the usual E u l e r angles are used to represent the o r i e n t a t i o n s . This l a t t e r difficulty was overcome by Evans and Murad (7) by using quaternions i n place of E u l e r angles. T h i s r e s u l t s i n a more efficient and accurate algorithm, and y i e l d s an order of magnitude improvement i n energy conservation. A f u r t h e r i n c r e a s e i n computing speed by a f a c t o r of three to eight may be obtained by using the m u l t i p l e time step method developed by S t r e e t t et al. (8). The a p p l i c a t i o n of t h i s method to r i g i d polyatomic molecules i s described i n another paper i n t h i s volume (9). In t h i s paper we use the quaternion method o f Evans and Murad to study methane, and report p r e l i m i n a r y r e s u l t s using a truncated and s h i f t e d s i t e - s i t e Lennard-Jones p o t e n t i a l . The m u l t i p l e time step method was not used i n t h i s work, p r i n c i p a l l y because o f storage l i m i t a t i o n s on the PDP 11/70 minicomputer used. The p o t e n t i a l model s t u d i e d i n v o l v e s f i v e Lennard-Jones s i t e s f o r each molecule, l o c a t e d on the C and H n u c l e i (C-H bond length 1.10 A). Thus the p a i r p o t e n t i a l i s the sum of 25 i n t e r a c t i o n s , 1 between C...C, 8 between C...H, and 16 between H...H sites. The p o t e n t i a l i s thus uCraUiU,) = I [ u ( r a 6

= 0

a ( J

) - u

o ( J

(r )] 0

r < T

Q

(1)

r > r

o where ri2 i s the v e c t o r from the center of molecule 1 to that o f

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

5.

MURAD

AND

63

Molecular Dynamics Simulation of Methane

GUBBINS

molecule 2, = 4>i0^^ i s the o r i e n t a t i o n o f molecule i , r g i s the distance between s i t e a of molecule 1 and s i t e 8 o f molecule 2, r = 10.025 A i s the c u t o f f d i s t a n c e , and a

Q

-

W

(2)

4 e

where e g and a g are the Lennard-Jones parameters (see Figure 1). Using t h i s p o t e n t i a l model, f i v e sets o f parameters (e 3>o* $) have been used to evaluate thermodynamic p r o p e r t i e s ( c o n f i g u r a t i o n a l energy, pressure and s p e c i f i c h e a t ) . The best o f these parameter sets i s used i n making a more d e t a i l e d study o f methane. This best set gives r e s u l t s s u p e r i o r to those o f a recent study (10) u t i l i z i n g the Williams (11) s i t e - s i t e exp-6 model, and i s also s u p e r i o r to the r e s u l t s obtained by Hanley and Watts (12) using an i s o t r o p i c m-6-8 p o t e n t i a l model. a

a


a

a

Method The quaternion method has been described i n d e t a i l i n recent p u b l i c a t i o n s CO so that we only give a b r i e f summary o f the method here. The quaternion parameters x> n, £ C can be defined i n terms of G o l d s t e i n E u l e r angles (13) by ?

X

= cos

(6/2)

n

= s i n (0/2)

c o s ( * + )/2, ^ c o s ( * - )/2,

(3)

> 5

= s i n (0/2)

s i n ( * - )/2,

C

= cos

sin(i|i + (J>)/2.

(0/2)

j

From (3) i t can be seen that X

2

+ n

2

+

C + C =1

(4)

2

The r o t a t i o n matrix A (13) defined by V = A • V -principal - -lab

(5)

i s given by

-S A =

2

I -2«n L

+

n

2

- C + x,

+ cx),

2(nc - 5x).

2Ux - 5n),

2

2

S - n 2

2

- c

-2(5c + nx).

2

+ x, 2

-S

2

2(nc + £x)

(6)

2(nx - SO

n

2

+ c

2

+ x

2

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

64

COMPUTER

MODELING

OF

MATTER

The p r i n c i p a l angular v e l o c i t y u>p i s r e l a t e d to the quaternions by the matrix equation W

' Px\ (7)


0 / The i n v e r s e of (7) i s simple to f i n d s i n c e the 4 x 4 matrix i s orthogonal. The equations of motion f o r these parameters are thus f r e e of s i n g u l a r i t i e s . A system o f 108 model methane molecules was s t u d i e d using standard p e r i o d i c boundary c o n d i t i o n s ( 1 ) . The i n t e r m o l e c u l a r forces and torques were c a l c u l a t e d i n the l a b o r a t o r y coordinate system. The torques were then converted to p r i n c i p a l torques Tp, using the r o t a t i o n matrix (equation ( 6 ) ) . P r i n c i p a l torques and f o r c e s F were used to evaluate the time d e r i v a t i v e s of the p r i n c i p a l angular v e l o c i t i e s and c a r t e s i a n p o s i t i o n s r of the molecules using the equations du) -77^ = T d t

/I Pa P

(a = x,y,z)

(8)

a

2

d r dP"

(9)

F/m

A t h i r d order p r e d i c t o r - c o r r e c t o r method (14) was used to i n t e grate the center of mass motion (equation ( 9 ) ) , while a second order method was used f o r o r i e n t a t i o n a l equations of motion of the molecules (equations (8) and ( 7 ) ) . A l l r e s u l t s were based on 2,000 time steps of 1.47 x 1 0 " seconds, a f t e r e q u i l i b r a t i o n . The energy and pressure were c o r r e c t e d f o r the long range p o t e n t i a l c o n t r i b u t i o n i n the usual way ( 1 ) , by p u t t i n g g g = 1 f o r r g > r . The s h i f t i n the p o t e n t i a l a l s o a f f e c t s the energy, but not the pressure or s p e c i f i c heat. The energy was c o r r e c t e d f o r t h i s p o t e n t i a l s h i f t by w r i t i n g U = + U h i f t > where Ujfl) i s the value corresponding to equations (1) and (2) and obtained i n the s i m u l a t i o n , and U h i f t i s the c o r r e c t i o n , given by 15

a

a

Q

s

s

where g ( r ) i s the s i t e - s i t e c o r r e l a t i o n f u n c t i o n . n

n


5.

MURAD

65


A N D GUBBINS

Energy conservation was w i t h i n about 0.75% per 10,000 time steps. One time step took approximately 25 seconds f o r a 108 molecule system on the PDP 11/70 minicomputer. Results The p o t e n t i a l o f equations (1) and (2) i n v o l v e s the para meter set ( e , a , ecc/ HH> Cc/°HH> ^CH* ^CH)» where ncH CcH g i the departure from the usual combining r u l e s f o r and CCH, e

c c

a

a

n

d

c c

v e

Q

CH

=

n

ke

CH Downloaded by EMORY UNIV on March 8, 2016 | http://pubs.acs.org Publication Date: June 1, 1978 | doi: 10.1021/bk-1978-0086.ch005

(a

* CH CC ^CH CC

+

a

HH

(10)

}

(11)

HH

In t h i s work we have s t u d i e d the e f f e c t o f v a r y i n g £QQ> °CC> e

e

a

n

d a

a

T

i e

v

a

u

e

s

a

n

a

r

e

e

CC^ HH CC^ HH* * l °f ^CH d CcH ^ P t fixed and are put equal to those given by Williams (11) i n h i s para meter set VII. Table I shows the Lennard-Jones parameter sets used i n t h i s work. In the f i r s t f i v e sets ZQQ and OQQ are kept f i x e d , and the r a t i o s ^QQ/O^H and ecc/ HH are v a r i e d from one set to the next. The parameters f o r set 1 were obtained by p u t t i n g the Lennard-Jones o* g and e g values equal to those used by Williams (11) i n h i s parameter set VII f o r the exp-6 model. Parameter sets 2-5 show increases i n £HH o f e i t h e r 10 o r 20%, and decreases i n of e i t h e r 0, 1, o r 2%. In a l l f i v e sets ncH n d were kept f i x e d at 0.972 and 1.132, r e s p e c t i v e l y , as given by Williams parameter set VII. For each o f the 5 parameter sets given i n Table I three s t a t e points were s t u d i e d and the mean squared d e v i a t i o n s between molecular dynamics and experimental r e s u l t s (15) were c a l c u l a t e d for the c o n f i g u r a t i o n a l energy, pressure and s p e c i f i c heat. The three s t a t e points used were: (a) T = 200 K, p = 10.00 mol fc" , (b) T = 225 K, p = 18.50 mol Z" , (c) T = 125 K, p = 25.27 mol ST . The f i r s t o f these c o n d i t i o n s i s i n the moderately dense super c r i t i c a l r e g i o n , the second i s i n the dense s u p e r c r i t i c a l r e g i o n , and the t h i r d i s i n the dense l i q u i d r e g i o n . Table I I gives the root mean squared d e v i a t i o n s found f o r each parameter s e t . Set 5 gave the best r e s u l t s . For t h i s parameter set s i x a d d i t i o n a l s t a t e c o n d i t i o n s were s t u d i e d . The r e s u l t i n g molecular dynamics r e s u l t s are compared with experiment i n Table I I I . The molecular dynamics r e s u l t s shown i n Table I I I may be used to r e s c a l e the values of SQQ and OQQ> keeping f i x e d CC/ HH> CC/ HH> ^CH ^CH* ^ t h o d o f c a r r y i n g out t h i s r e s c a l i n g i s i l l u s t r a t e d ( f o r two s t a t e c o n d i t i o n s ) i n Figure 2. For a s i n g l e s t a t e c o n d i t i o n (e.g. 1 i n Figure 2) i t i s p o s s i b l e to f i n d a set o f (OQQ CQQ) values that give exact agreement between molecular dynamics and experiment f o r the c o n f i g u r a t i o n a l e

a

a

a

1

1

e

G

a

a

a

n

d

T

h

1

e m

9



66

COMPUTER

MODELING

OF MATTER

Figure 1. Site-site and center-center distances for two CH, molecules at orientations o> and w t

t

2

Table I P o t e n t i a l Parameter Sets Studied

1 2 3 4 5 Rescaled 5

cc

CH (K)

HH (K)

"CH (X)

HH (A)

350 350 350 350 350

023 023 008 008 2.995

2.870 2.870 2.841 2.841 2.813

48.760 48.760 48.760 48.760 48.760

20.690 21.700 21.700 22.665 22.665

850 535 535 220 220

3.350

2.995

2.813

51.198

23.798

8.631

CC (A)

Parameter Set

e

/k

(K)

£

/ K

E

/ K

The b o n d l e n g t h i n parameter sets 1-5 i s kept f i x e d at 1.10 A. The parameter s e t " r e s c a l e d 5" a l s o has t h i s same value o f the bond l e n g t h , s i n c e the r e s c a l e d a's are unchanged. o


5.

MURAD

AND


GUBBINS

Table I I


Root Mean Squared Deviations For Models Studied

Root Mean Squared Deviation"*"

Parameter Set

C

U /NkT

P/pkT

Cy/Nk

1 2 3 4 5

0.45 0.41 0.26 0.18 0.22

0.46 0.41 0.34 0.21 0.14

0.43 0.35 0.44 0.37 0.37

0.45 0.40 0.30 0.21 0.21

Rescaled 5

0.08

0.07

0.22

0.10

V

= C V

- C

V

i d g

a

s

, where C

V

and C

i d g

T o t a l ft

a

s

are a t the

v

V

same T. We expect the e r r o r s i n our MD values to be c

no greater than 0.03 f o r U /NkT and 0.08 f o r P/pkT.

ft

For the t o t a l d e v i a t i o n , energy, pressure and s p e c i f i c heat were given the weights 4:2:1.

Table I I I Results f o r Parameter Set 5 Before R e s c a l i n g

M „ Temperature Density (K) (mol I ) n

C

U /NkT MD

Ex£

r

P/pkT MD

Ex£

C /Nk v MD Ex£

198 228 250

10.00 10.00 10.00

-1.49 -1.70 0.48 0.35 0.37 0.70 -1.28 -1.45 0.70 0.55 0.15 0.55 -1.13 -1.25 0.77 0.66 0.12 0.39

223 248 276

18.50 18.50 18.50

-2.45 -2.59 0.92 0.80 0.32 0.39 -2.15 -2.28 1.22 1.08 0.37 0.37 -1.92 -1.99 1.45 1.25 0.52 0.40

123 151 182

25.27 25.27 25.27

-6.68 -6.84 -0.10 0.0 0.55 0.85 -5.28 -5.41 1.23 1.15 0.73 0.98 -4.27 -4.33 2.19 1.91 0.72 0.79


67


68

COMPUTER

MODELING O F MATTER

Figure 2. Rescaling the potential parameters acc and e r, keeping the ratios VCC'-VCH'-VHH and tcc'>*cH-*nHfixed.Curves 1 and 2 correspond to different state conditions. C

Table IV Results f o r Parameter Set 5 A f t e r R e s c a l i n g

c

U /NkT

r

P/pkT

C /Nk v MD Exp

Temperature (K)

Density (mol J T )

MD

Exp

MD

Exp

208 239 263

10.00 10.00 10.00

-1.49 -1.28 -1.13

-1.59 -1.34 -1.12

0.48 0.70 0.77

0.41 0.60 0.71

0.37 0.15 0.12

0.70 0.44 0.35

234 260 290

18.50 18.50 18.50

-2.45 -2.15 -1.92

-2.45 -2.15 -1.92

0.92 1.22 1.45

0.90 1.13 1.34

0.32 0.37 0.52

0.37 0.37 0.43

129 158 191

25.27 25.27 25.27

-6.68 -5.28 -4.27

-6.52 -5.15 -4.18

-0.10 1.23 2.19

0.30 1.35 2.20

0.55 0.73 0.72

0.89 0.96 0.80

1



5.

MURAD


A N D GUBBINS

1

1

O-O

100

'

'

180

1

'

260

1

— 340

T(K)

Figure 3. Configurational internal en&J from molecular dynamics (points) and experiment (lines)

eT

1

1 00

i

180

1 T(K)

i

260

1

1—

3 40

Figure 4. Pressure from molecular dynamics (points) and experiment (lines)


69

70

COMPUTER

MODELING

OF MATTER

energy. This y i e l d s a l i n e on a p l o t of a^ vs. £QC* A second, but d i f f e r e n t , s e t o f (OQQ, cc) values i s obtained by equating molecular dynamics and experimental values f o r the pressure at the same s t a t e c o n d i i t o n . T h i s y i e l d s a second l i n e on the CC * CC P l « The i n t e r s e c t i o n o f these two l i n e s gives a unique s e t o f ( QQ, values f o r the s t a t e c o n d i t i o n s t u d i e d . By t r e a t i n g the molecular dynamics data f o r the other s t a t e c o n d i t i o n s i n a s i m i l a r way i t i s p o s s i b l e to estimate a best set o f (GQQ, cc) l from the r e s u l t i n g " i n t e r v a l o f c o n f u s i o n " on a OQQ V S . €QQ p l o t . T h i s best s e t of r e s c a l e d parameters i s given i n Table I , and the r e s u l t i n g mean squared d e v i a t i o n s are given i n Table I I . Table IV and Figures 3 and 4 compare the molecular dynamics values with experiment f o r t h i s rescaled p o t e n t i a l . The molecular dynamics r e s u l t s found from parameter s e t 1 are almost the same as those found f o r the Williams p o t e n t i a l using h i s s e t V I I , and reported elsewhere (10). The r e s u l t s f o r the r e s c a l e d Lennard-Jones p o t e n t i a l used here are sub s t a n t i a l l y b e t t e r , as seen from Table I I . We p l a n to study f u r t h e r refinements of t h i s p o t e n t i a l model, i n c l u d i n g the a d d i t i o n o f an octupole moment, and we are a l s o i n the process of i n v e s t i g a t i n g other p r o p e r t i e s , i n c l u d i n g the o r i e n t a t i o n a l c o r r e l a t i o n f u n c t i o n s , s e l f - d i f f u s i o n c o e f f i c i e n t s , and time correlation functions. c

e

a

v s

e

o t

0


e

v

a

u

e

s

Acknowledgments We thank the American Gas A s s o c i a t i o n , the N a t i o n a l Science Foundation, and the Petroleum Research Fund, as administered by the American Chemical S o c i e t y , f o r support o f t h i s work. L i t e r a t u r e Cited 1.

S t r e e t t , W. B. and Gubbins, K. E., Ann. Rev. Phys. Chem., (1977), 28, 373.

2.

S t i l l i n g e r , F. H., Adv. Chem. Phys., (1975), 12, 107.

3.

McDonald, I . R. and 64, 4790.

4.

Kushick, J . and Berne, B. J . , J . Chem. Phys., (1976), 64, 1362.

5.

Evans, D. J . and Watts, R. O., Mol. Phys., (1976), 32, 93.

6.

Ryckaert, J . P. and Bellemans, 30, 123.

Klein,

M. E., J . Chem. Phys.,

(1976),

A., Chem. Phys. L e t t . ,


(1975),


5.

MURAD

A N D GUBBINS


71

7.

Evans, D. J., Mol. Phys., (1977), 34, 317; Evans, D.J.and Murad, S., ibid., (1977), 34, 327.

8.

S t r e e t t , W. B., T i l d e s l e y , D. J . and (1978), 35, 639.

Saville,

G., Mol. Phys.,

9.

S t r e e t t , W. B., T i l d e s l e y , D. J. and t h i s volume.

Saville,

G., (1978),

10.

Murad, S., Evans, D. J., Gubbins, K. E., S t r e e t t , W. B. and T i l d e s l e y , D. J., Mol. Phys., (1978), submitted.

11.

W i l l i a m s , D. E., J. Chem. Phys.,

12.

Hanley, H. J. M. and Watts, R. O., Mol. Phys., 1907; Aust. J . Phys., (1975), 28, 315.

13.

G o l d s t e i n , H., " C l a s s i c a l Mechanics", Chapter 4, AddisonWesley, (1971).

14.

Gear, C. W., "Numerical Initial Value Problems i n Ordinary Differential Equations", P r e n t i c e - H a l l , Englewood Cliffs (1971).

15.

Goodwin, R. D., "The Thermodphysical P r o p e r t i e s o f Methane from 90 to 500 K at Pressures to 700 Bars", N.B.S. T e c h n i c a l Note No. 653 (National Bureau o f Standards, Washington, DC, 1974).

RECEIVED

(1967), 47, 4680. (1975), 29.,

August 15, 1978.


Molecular Dynamics Simulation of Methane Using a Singularity-Free

Recommend Documents