Molecular dynamics simulations of fluid carbon dioxide using the model potential based on ab initio MO calculation

Computational Materials Science 14 (1999) 220±226

Molecular dynamics simulations of ¯uid carbon dioxide using the model potential based on ab initio MO calculation Seiji Tsuzuki *, Kazutoshi Tanabe National Institute of Materials and Chemical Research, Tsukuba, Ibaraki 305, Japan

Abstract The nonbonding interaction energy of carbon dioxide dimer was calculated using several basis sets to evaluate the basis set e€ect. Large basis sets including di€use polarized functions were necessary to evaluate the interaction energy correctly. Small basis sets considerably underestimate the interaction energy. Nonbonding parameters of a model potential of carbon dioxide were re®ned based on the interaction energies of 40 geometrical con®guration dimers calculated by an MP2/6-311 + G(2df)-level ab initio method with the BSSE correction. The molar volume and heat of evaporation of liquid carbon dioxide obtained from molecular dynamics simulations with the model potential reproduced the experimental values with only a few percentage errors. On the other hand the model potential based on the BSSE uncorrected interaction energies overestimated the attractive interaction and failed to reproduce the experimental values. The same model potential also reproduced the experimental pressure and self-di€usion coecient of super critical ¯uid satisfactorily. Ó 1999 Published by Elsevier Science B.V. All rights reserved. Keywords: Ab initio molecular orbital calculation; Model potential; Molecular dynamics simulation; Carbon dioxide; Supercritical ¯uid; Nonbonding interaction; van der Waals parameters

1. Introduction Molecular dynamics simulation has turned out to be an important methodology in the study of condensed matter [1]. Starting from the knowledge of molecular interactions, the simulation provides structural information as well as thermodynamics and transport properties of liquid and solid states. The reliance of the simulation depends on the interaction potential used in the simulation [2]. However it is a dicult task to re®ne the nonbonding parameters accurately only from experimental data. A number of nonbonding interaction * Corresponding author. Tel.: +81 298 544522; fax: +81 298 544487; e-mail: [email protected].

potentials have been proposed for carbon dioxide. However these potentials have several de®ciencies [2±5]. In this work we re®ne the parameters of nonbonding interaction potential of carbon dioxide based on the interaction energies obtained by ab initio molecular orbital calculations.

2. Method The Gaussian 92 program [6] was used for ab initio calculations. The structure of a single carbon dioxide molecule was optimized at the MP2/6311G(3d) level. The optimized C±O bond distance  The optimized structure was used for is 1.169 A. the calculation of the dimers. The intermolecular

0927-0256/99/$ ± see front matter Ó 1999 Published by Elsevier Science B.V. All rights reserved. PII: S 0 9 2 7 - 0 2 5 6 ( 9 8 ) 0 0 1 1 1 - 6

S. Tsuzuki, K. Tanabe / Computational Materials Science 14 (1999) 220±226

interaction energies were calculated with a 6311 + G(2df) basis set and the electron correlation correction by the Mùller±Plesset method [7,8] and the basis set superposition error (BSSE) [9] correction by the counterpoise method [10]. Molecular dynamics simulation of 216 rigid carbon dioxide molecules was carried out in the NPT ensemble for liquid. The equations of motion of rigid bodies were integrated by the leap-frog algorithm with a time step of 2 fs. Periodic boundary conditions with cubic unit cells were employed. Constant temperature and pressure conditions were maintained by using the method of Nose and that of Andersen [11±13]. The system was equilibrated initially by 10-ps simulation. The density and nonbonding potential energy were obtained from the following 100-ps simulation. The heat of evaporation DHv was calculated from the nonbonding potential per molecule, Ei , by DHv ˆ jEi j ‡ RT [14]. Molecular dynamics simulation of 216 rigid carbon dioxide molecules was carried out in the NVT ensemble for supercritical ¯uid. The system was equilibrated initially by 20-ps simulation. The pressure was obtained from the following 50-ps simulation. The self-di€usion coecient was obtained by 100-ps simulation after equilibration. 3. Result and discussion 3.1. Electron correlation energy To evaluate the e€ect of electron correlation, the intermolecular interaction energies of the T shape dimer (Fig. 1, dimer C) were calculated by

221

the Hartree±Fock and by the second, third and fourth-order Mùller±Plesset method using the 6311G(2d) basis set. The calculated interaction energies at the HF, MP2, MP3 and MP4 (SDTQ) levels were ÿ0.29, ÿ0.78, ÿ0.67 and ÿ0.78 kcal  molÿ1 , respectively, at the C/C distance of 4.2 A. The correction of electron correlation is important to evaluate the interaction. The dispersion energy is not negligible. The HF level calculation underestimates the attractive interaction, due to the lack of the evaluation of the dispersion energy. The calculated interaction energy at the MP2 level is close to that at the MP4 level. Due to the good performance of the MP2 level calculation, we decided to correct the electron correlation energies by the MP2 method in all the following calculations. 3.2. The e€ect of basis set The intermolecular interaction energy of the T shape dimer (dimer C) was calculated using several basis sets to evaluate the basis set e€ect. The calculated interaction energies greatly depend on the used basis set as shown in Table 1. Small basis sets (6-31G* and 6-311G*) greatly underestimate the attractive interaction. 3.3. Interaction potential The interaction potentials of the four orientations of carbon dioxide dimers (Fig. 1) were calculated at the MP2/6-311 + G(2df) level. The calculated potentials were compared with those obtained by Williams' force ®eld [15] and by Kuchta's force ®eld [4] as shown in Figs. 2±5. The

Fig. 1. Four types of orientation of carbon dioxide dimers considered in this work.

222

S. Tsuzuki, K. Tanabe / Computational Materials Science 14 (1999) 220±226

Table 1 Calculated interaction energies of T shape dimer using several basis set

a

Scale factor of exponents of gaussians

Basis set

fd1

fd2

fd3

fd4

b

ff1

A

6-31G*

B

6-311G*

1.0

C

6-311G(2d)

0.5

2.0

D

6-311G(2df)

0.5

2.0

E

6-311G(3d)

0.25

1.0

4.0

F

6-311G(3df)

0.25

1.0

4.0

G

6-311G(4d)

0.125

0.5

2.0

H

6-311G(2d)

0.25

1.0

I

6-311 + G(3d)

0.25

1.0

J

6-311G(2df)

0.25

1.0

0.5

K

6-311G(2df)

0.25

1.0

0.25

L

6-311G(2df)

0.25

1.0

0.125

M

6-311G(2d2f)

0.25

1.0

0.125

N

6-311 + G(2df)

0.25

1.0

0.25

O

6-311G(3df)

0.0625

0.25

Interaction energies ff2

1.0

1.0 8.0

4.0

1.0

0.25

0.5

MP2

c

ÿ0.67 (ÿ1.35) ÿ0.63 (ÿ1.15) ÿ0.78 (ÿ1.37) ÿ0.84 (ÿ1.35) ÿ0.96 (ÿ1.71) ÿ0.99 (ÿ1.67) ÿ0.95 (ÿ2.02) ÿ0.97 (ÿ1.70) ÿ0.97 (ÿ1.41) ÿ1.04 (ÿ1.71) ÿ1.05 (ÿ2.19) ÿ1.01 (ÿ2.95) ÿ1.07 (ÿ2.66) ÿ1.07 (ÿ1.88) ÿ1.07 (ÿ2.87)

BSSE

HF

disp

0.68

ÿ0.47

ÿ0.20

0.52

ÿ0.54

ÿ0.09

0.59

ÿ0.29

ÿ0.49

0.52

ÿ0.37

ÿ0.47

0.75

ÿ0.30

ÿ0.66

0.67

ÿ0.34

ÿ0.66

1.07

ÿ0.23

ÿ0.71

0.73

ÿ0.32

ÿ0.65

0.44

ÿ0.28

ÿ0.69

0.67

ÿ0.35

ÿ0.69

1.14

ÿ0.30

ÿ0.75

1.94

ÿ0.28

ÿ0.73

1.59

ÿ0.32

ÿ0.75

0.81

ÿ0.30

ÿ0.77

1.80

ÿ0.30

ÿ0.77

 Intermolecular distance is 4.2 A. Exponents of the d functions of carbon atom adi (C) are given by 0:626  fdi . Exponents of the d functions of oxygen atom adi (O) are given by 1:292  fdi . Exponents of the f functions of carbon atom afi (C) are given by 0:8  ffi . Exponents of the f functions of oxygen atom afi (O) are given by 1:4  ffi . c The calculated interaction energies with the BSSE correction. The values in parentheses correspond to the interaction energies not corrected for the BSSE. a

b

Williams' force ®eld overestimates the repulsion. The Kuchta's potential agrees with the potential obtained by the ab initio calculations better than the Williams' potential. On the other hand the Kuchta's force ®eld underestimates the repulsion in the linear orientation. 3.4. Re®nement of nonbonding parameters The parameters of the Buckingham type potential [16] Eq. (1) were re®ned to ®t the interac-

tion energies of the 40 dimers calculated by the ab initio method. Eij ˆ Bij exp…ÿCij rij † ÿ Aij rijÿ6 ‡ qi qj rijÿ1 :

…1†

The charges ÿ0.34 and 0.68 e (obtained by the ®tting of electrostatic potential from the MP2/6311 + G(2df) level wave function with the Merz± Singh±Kollman scheme [17]) were put on oxygen and carbon atoms, respectively. The nonbonding interaction center of oxygen atom was relocated toward the connected carbon atom by 9% of the

S. Tsuzuki, K. Tanabe / Computational Materials Science 14 (1999) 220±226

223

Fig. 2. Comparison of various intermolecular interaction potential of dimer A.

Fig. 3. Comparison of various intermolecular interaction potential of dimer B.

C±O bond distance to incorporate the anisotropy of the nonbonding interaction of oxygen atoms as shown in Fig. 6. The optimized Buckingham parameters are shown in Table 2 (Param 1). The incorporation of the anisotropy is important. The relocation decreases the RMS error considerably. The ®tted potential agrees with the interaction energies obtained by the ab initio method better than those from the Williams' and Kuchta's potentials.

bon dioxide with constant temperature and pressure conditions. The calculated molar volume and heat of evaporation at three di€erent temperatures agree well with the experimental values [18] as shown in Table 3. The calculated molar volumes are 1±2% larger than the experimental values. The calculated heats of evaporation are 6±9% larger than the experimental values. The same parameters were employed for the simulation of supercritical ¯uid with constant temperature and volume conditions. The calculated pressures and self-di€usion coecients reproduced the experimental values [18,19] satisfactorily as shown in Tables 4 and 5. The good agreement shows that the used model potential

3.5. MD simulation of ¯uid carbon dioxide The newly re®ned parameters were used for the molecular dynamics simulation of liquid car-

224

S. Tsuzuki, K. Tanabe / Computational Materials Science 14 (1999) 220±226

Fig. 4. Comparison of various intermolecular interaction potential of dimer C.

Fig. 5. Comparison of various intermolecular interaction potential of dimer D.

from the ab initio molecular orbital calculation is reliable. 3.6. Basis set superposition error The BSSE's of the dimer C using several basis sets are shown in Table 1. The BSSE's are large, even if large basis sets are used. The accuracy of the counterpoise method has been a controversial issue [20]. To compare the accuracy of the calculated intermolecular interaction energies with and without the BSSE correction, we also carried out the molecular dynamics simulation using the model potential re®ned based on the BSSE uncorrected interaction energies calculated at the

Fig. 6. Anisotropic oxygen nonbonding interaction model used in this work.

S. Tsuzuki, K. Tanabe / Computational Materials Science 14 (1999) 220±226 Table 2 Optimized parameters of the Buckingham equation Parameter set Param 1 C/C O/O C/O

b

Param 2 C/C O/O C/O

c

B (kcal molÿ1 )

C ÿ1 ) (A

356 389 372

217 354000 8760

2.27 4.40 3.16

451 802 601

108 452000 6990

2.00 4.44 2.98

A 6 (kcal A molÿ1 )

225

MP2/6-311 + G(2df) level. The parameters of the model potential are shown in Table 2 (Param 2). The parameters A of the Buckingham equation correspond to the strength of the attractive interaction. The parameters A of the parameter set 2 are much larger than those of the parameter set 1 due to the deeper potential depth of the BSSE uncorrected potential. The calculated heat of evaporation using the parameter set 2 is about three times as large as the experimental value as shown in Table 3. Clearly this model potential overestimates the attractive interaction. The poor performance of this model potential shows that the BSSE uncorrected interaction energies are not reliable.

a

RMS error (kcal molÿ1 ) 0.141

0.174

a The charges ÿ0.34 and 0.68 e (obtained by the ®tting of electrostatic potential from the MP2/6-311 + G(2df) level wave function with the Merz±Singh±Kollman scheme) were put on oxygen and carbon atoms, respectively. b Parameters were re®ned based on the BSSE corrected interaction energies. van der Waals center of oxygen atom is relocated toward the connected carbon atom 9% of the C±O bond distance. c Parameters were re®ned based on the BSSE uncorrected interaction energies. van der Waals center of oxygen atom is relocated toward the connected carbon atom 10% of the C±O bond distance.

4. Conclusion We have evaluated the basis set and electron correlation e€ects on the calculated intermolecular interaction energies of carbon dioxide dimers and have found that large basis set including di€use

Table 3 Molar volume and heat of evaporation of liquid carbon dioxide obtained by molecular dynamics simulation Param set

Temp

a

Press

a

Molar volume Calc

1 1 1 2

220.0 230.0 240.0 220.0

5.2 8.6 14.0 8.3

38.4 39.8 40.9 25.5

b

Exp (+1.9%) (+2.1%) (+1.2%) (ÿ32.4%)

b

Heat of evap c

Calc

37.7 39.0 40.4 37.7

3.8 3.7 3.6 10.9

Exp (+5.6%) (+8.8) (+9.1) (+202.8)

c

3.6 3.4 3.3 3.6

a

The average temperature in K and average pressure in 105 Pa in the simulation. Molar volume in cm3 molÿ1 and heat of evaporation in kcal molÿ1 . The di€erences from the experimental values are shown in parentheses. c Ref. [18]. b

Table 4 Calculated pressure of ¯uid carbon dioxide obtained by molecular dynamics simulation Temperature (K) 320 330 350 400 a

a

ÿ3

Density (mol cm ) 0.004

0.006

0.010

0.015

6.60 (6.94) 7.54 (7.43) 8.38 (8.40) 10.6 (10.7)

8.15 (8.45) 9.29 (9.29) 10.9 (10.9) 14.6 (14.8)

8.83 (9.94) 10.0 (11.6) 13.9 (14.8) 21.3 (22.8)

10.8 12.9 18.8 37.3

(12.6) (15.8) (22.2) (38.6)

Simulations were carried out under constant temperature and volume conditions. Pressure in 106 Pa. Experimental values are in parentheses. Ref. [18].

226

S. Tsuzuki, K. Tanabe / Computational Materials Science 14 (1999) 220±226

Table 5 Calculated self-di€usion coecient

References

a

Density (mol cmÿ3 )

2 fsÿ1 ) Self-di€usion coecient (A calc

obsb

0.00417 0.00573 0.00940 0.01409 0.01937

1.16 0.764 0.457 0.267 0.151

1.198 0.875 0.523 0.287 0.189

a b

Temperature is at 323.15 K. Ref. [19].

polarization functions are important to evaluated the attractive interaction. We have re®ned Buckingham type nonbonding interaction potential parameters based on ab initio calculations using reasonably large basis set with electron correlation and BSSE correction. The calculated properties of ¯uid carbon dioxide with this model potential reproduced the experimental values, suggesting that the model potential is accurate. Acknowledgements We thank the Research Information Processing System, Tsukuba, Japan for the support of computational facilities.

[1] M.P. Allen, D.L. Tildesley, Computer Simulation of Liquids, Clarendon Press, New York, 1987. [2] C.S. Murthy, S.F. O'Shea, I.R. McDonald, Mol. Phys. 50 (1983) 531. [3] B. Kuchta, R.D. Etters, Phys. Rev. B 38 (1988) 6265. [4] R.D. Etters, B. Kuchta, J. Chem. Phys. 90 (1989) 4537. [5] K.B. Domanski, O. Kitao, K. Nakanishi, Mol. Simul. 12 (1994) 343. [6] M.J. Frisch, G.W. Trucks, M. Head-Gordon, P.M.W. Gill, M.W. Wong, J.B. Foresman, B.G. Johnson, H.B. Schlegel, M.A. Robbs, E.S. Replogle, R. Gomperts, J.L. Andres, K. Raghavachari, J.S. Binkley, C. Gonzalez, R.L. Martin, D.J. Fox, D.J. Defrees, J. Baker, J.J.P. Stewart, J.A. Pople, Gaussian 92, Gaussian, Pittsburgh, PA, 1992. [7] C. Mùller, M.S. Plesset, Phys. Rev. 46 (1934) 618. [8] J.A. Pople, R. Seeger, R. Krishnan, Int. J. Quantum. Chem. Symp. 11 (1977) 149. [9] B.J. Ransil, J. Chem. Phys. 34 (1961) 2109. [10] S.F. Boys, F. Bernardi, Mol. Phys. 19 (1970) 553. [11] S. Nose, Mol. Phys. 52 (1984) 255. [12] S. Nose, J. Chem. Phys. 81 (1984) 511. [13] H.C. Andersen, J. Chem. Phys. 72 (1980) 2384. [14] W.L. Jorgensen, J. Am. Chem. Soc. 103 (1981) 335. [15] S.R. Cox, L.Y. Hsu, D.E. Williams, Acta Crystallogr. A 37 (1981) 293. [16] A.D. Buckingham, B.D. Utting, Annu. Rev. Phys. Chem. 21 (1970) 287. [17] B.H. Besler, K.M. Merz, Jr., P.A. Kollman, J. Comput. Chem. 11 (1990) 431. [18] S. Angus, B. Armstrong, K.M. deReuck, Carbon Dioxide, International Thermodynamic Tables of the Fluid State 3, Pergamon Press, Oxford, 1973. [19] P. Etesse, J.A. Zega, R. Kobayashi, J. Chem. Phys. 97 (1992) 2022. [20] F.-M. Tao, Y.-K. Pan, J. Phys. Chem. 95 (1991) 3582.