New Interatomic Potentials For Silica

Ersan Demiralp, Tahir Çagin, Norman T. Huff*, William A. Goddard III

Materials and Process Simulation Center, Beckman Institute (139-74)

Division of Chemistry and Chemical Engineering

California Institute of Technology, Pasadena, California 91125

* Owens Corning Science & Technology Center,

2790 Columbus Road, Granville, OH 43023-1200

 

Abstract

A new interatomic potential for the description of various forms of silica has been developed. This new potential uses Morse type functions to describe short range interactions while the long range electrostatic interactions are taken as purely Coulombic. The atomic charges are calculated using the Charge Equilibration method of Rappe and Goddard. Thus, the charges vary depending upon the other atoms in the vicinity. The accuracy of the crystalline cell parameters of 7 polymorphs of crystalline silica calculated using this potential as well as previously published potentials are presented. This potential appears to qualitatively predict high pressure phase transitions in crystalline silica. In addition, when used in NPT molecular dynamics simulations at atmospheric pressures, it appears to predict the structural features of amorphous silica.

 

1. Introduction

Silica can exist in several crystalline polymorphs as well as in the amorphous state. A number of interatomic potential functions have been developed which can, with reasonable accuracy, reproduce the structural characteristics of one or more of the silica polymorphs. For example, the two body interatomic potential of Tsuneyuki et. al. (TTAM potential) predicts with a fair degree of accuracy the crystal cell parameters of 7 forms of silica. ' However, the predictions of silica melts and amorphous silica are not in good agreement with experiment. Three body interatomic potentials of Garofalini et. al. do a good job of predicting amorphous silica structures (using NVT molecular dynamics (MD) simulations) but predicted crystal cell parameters of various polymorphs is not as accurate.

In this work, a new type of two body interatomic potential was developed which predicts with fair accuracy, structural characteristics of both crystalline and amorphous polymorphs of silica. This new interatomic potential is different from previous potentials in that a Morse type functionality is assumed for short range forces. Long range interactions are represented by purely Coulombic interactions. However, the charges used to calculate the Coulombic interactions are variable. The charge on individual atoms vary according to the location and type of surrounding atoms according to the Charge Equilibration (QEq) method of Rappe and Goddard.1

 

This work uses a two body interatomic potential function with the form:

 

 

 

Uij (Rij) = + Do (1)

 

 

In this equation, qi and qj represent the charges on the ith and jth atoms. Rij is the distance between atoms i and j. Do, Ro and g are adjustable parameters. Do is generally associated with a bond strength, Ro is associated with an equilibrium interatomic separation distance while g is related to curvature at the minimum of the potential. The parameters for each of the interatomic potentials are shown in Table 1.

 

Table 1. Morse Potential Parameters.

 

Interaction

Ro (Å)

Do (kcal/mol)

g

O-O

3.7910

0.5363

10.4112

Si-Si

3.7598

0.17733

15.3744

Si-O

1.6280

45.9970

8.6342

 

The values of Do, Ro, and g were calculated based mainly upon experimental data for b -cristobalite structures and slightly modified to obtain better descriptions for the other silica polymorphs.

The charges for optimized silica polymorphs varied between 1.22 and 1.4 for Si and -0.61 to -0.70 for O. This compares with potential derived charges and Mulliken charges obtained from a Si(OH)4 cluster ab initio Hartree Fock calculation using a 6-31G** basis set of 1.528 and 1.493 respectively. Thus, the QEq charges appear reasonable.

 

2.Results

The crystal cell parameters calculated using this new interatomic potential are compared with experimental values as well as with the cell parameters calculated using the TTAM potential and the potential in the Glass 2.01 force field (gff 2.01). Gff 2.01 potential in the Cerius2 molecular modeling software developed by Molecular Simulations Incorporated (MSI) is very similar to the three body potential function of Garofalini et al.4 Some results of calculated crystalline structures are presented in Table 2. Predicted densities of both crystalline and amorphous polymorphs are presented in Table 3. Calculated bulk moduli of the crystalline and amorphous polymorphs are presented in Table 4.

 

 

 

 

 

 

Table 2. Calculated and experimental structural parameters.

 

a -quartz

A

B

C

a

b

g

Experimental

4.9160

4.9160

5.4050

90

90

120

This Work

5.0973

5.0973

5.5777

90

90

120

TTAM

4.9893

4.9893

5.5269

90

90

120

gff 2.01

5.1508

5.1508

5.6341

90

90

120

a -cristobalite

           

Experimental

4.9570

4.9570

6.8903

90

90

90

This Work

5.0588

5.0588

6.5601

90

90

90

TTAM

4.9640

4.9640

6.6365

90

90

90

gff 2.01

5.2836

5.2836

7.4721

90

90

90

b -cristobalite

           

Experimental

7.1600

7.1600

7.1600

90

90

90

This Work

7.1339

7.1339

7.1339

90

90

90

TTAM

7.0537

7.0537

7.0537

90

90

90

gff 2.01

7.4721

7.4721

7.4721

90

90

90

Stishovite

           

Experimental

4.1801

4.1801

2.6150

90

90

90

This Work

4.2610

4.2610

2.7873

90

90

90

TTAM

4.2584

4.2584

2.7445

90

90

90

gff 2.01

4.1759

4.1759

3.0261

90

90

90

 

Table 3. Experimental and calculated densities (g/cm3) for silica polymorphs.

 

Polymorph

This Work

TTAM

gff 2.01

Experimental

Glass

2.23

2.30

2.16

2.19

a -quartz

2.3849

2.5121

2.3123

2.6458

a -cristobalite

2.3772

2.4405

1.9133

2.3572

b -cristobalite

2.1985

2.2744

1.9133

2.1746

Stishovite

3.9432

4.0096

3.7815

4.2808

b -tridymite

2.0338

2.2085

1.9136

2.1833

Keatite

2.3533

2.5668

2.2525

2.4987

Coesite

2.6036

2.7232

2.7319

2.9215

 

 

 

 

 

 

 

 

 

Table 4. Experimental (300 K) and calculated (0 K) bulk modulus (Gpa) for silica polymorphs.

 

Polymorph

This Work

TTAM

gff 2.01

Experimental

a -quartz

30.05

34.60

208.52

37.2 ± 0.2

a -cristobalite

19.82

18.67

197.86

16.8

b -cristobalite

20.64

19.38

197.86

 

Stishovite

254.66

327.15

340.34

295 ± 2

b -tridymite

14.61

16.47

196.53

 

Keatite

25.96

38.88

86.21

 

Coesite

79.77

99.50

130.06

94 ± 1

 

 

Amorphous structures were generated by performing NVT, followed NPT MD simulations with 648 atoms. The starting structure for these simulations was a b -cristobalite structure containing 27 unit cells (3 unit cells in x, y and z directions). The MD simulation consisted of 20 ps of NVT dynamics at 8000 K followed by 60 ps at 4000 K, 4 ps at 100 K increments from 3900 K to 400 K and then 10 ps at 300 K. This was followed by 30 ps of NPT dynamics at 300 K. This final trajectory was used to generate densities, RDF, and angle distributions. The angle distributions for O-Si-O and Si-O-Si for these potentials are shown in Figure 1.

 

Figure 1: Bond Angle Distributions

 

 

 

 

 

 

 

The peak values for O-Si-O are 108° ,108° , 106° and for Si-O-Si are 145° , 154° , 147° from the glass simulations by using the potential described in this paper, the gff 2.01 and the TTAM potentials respectively. The experimental values for O-Si-O and Si-O-Si angles were found to be 109° 47' and 144° by Mozzi and Warren. It is interesting to note that for the O-Si-O angle all these potentials have almost the same distributions. However, for Si-O-Si angle distributions, our potential had the peak at 145° , but gff 2.01 had the peak at 154° . The widths of the curves at the half of maximum value (FWHM) are 37° for all the potentials.

 

Discussion and Conclusions

This new two body potential appears to predict structural parameters of both crystalline and amorphous silica fairly well. The prediction of the density and angle distributions for amorphous silica using NPT dynamics is encouraging. This suggests that the main forces present in crystalline and amorphous silica appear to be ionic in nature as opposed to covalent. However, the broadening of the Si-O-Si and O-SI-O bond angle distributions suggest that in amorphous materials, the covalent bonding plays a significant role.

 

References