Bond-order reactive force fields for molecular dynamics simulations of crystalline silica

Computational Materials Science 111 (2016) 269–276

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Bond-order reactive force fields for molecular dynamics simulations of crystalline silica Benjamin J. Cowen a,b, Mohamed S. El-Genk a,b,c,d,⇑ a

Institute for Space and Nuclear Power Studies, University of New Mexico, Albuquerque, NM, USA Nuclear Engineering Department, University of New Mexico, Albuquerque, NM, USA c Mechanical Engineering Department, University of New Mexico, Albuquerque, NM, USA d Chemical and Biological Engineering Department, University of New Mexico, Albuquerque, NM, USA b

a r t i c l e

i n f o

Article history: Received 29 May 2015 Received in revised form 10 September 2015 Accepted 16 September 2015

Keywords: Bond-order Variable-charge reactive force fields MD simulation of silica polymorphs Alpha–beta transition Accuracy and transferability

a b s t r a c t This paper investigates the applicability of the bond-order, variable-charge (BOVC) force fields of the H2 O , and ReaxFFGSI Charge-Optimized Many-Body (COMB10), ReaxFFSiO SiO , for molecular dynamics (MD) simulations of crystalline SiO2. The calculated lattice constants and densities of the four SiO2 polymorphs, quartz, cristobalite, coesite, and stishovite, are compared to experimental values. Additionally, the calculated pair distribution and bond-angle distribution functions and the a–b transition for quartz, the most stable low-energy polymorph, are compared to experimental results. The simulations with the COMB10 force field accurately predict the properties of the SiO2 polymorphs, except the a-cristobalite, and the GSI 2O quartz a–b transition. The results with ReaxFFH SiO and ReaxFFSiO accurately predict the properties of 2O the SiO2 polymorphs, except the stishovite, but those with ReaxFFH SiO inaccurately predict the quartz a–b transition. Ó 2015 Elsevier B.V. All rights reserved.

1. Introduction Characterizing the interface of silica (SiO2) with other materials, such as silicon (Si), water (H2O), and oxygen (O), is important in a wide range of applications. Examples are the interaction of SiO2 with biomolecules, such as H2O, in chromatography, geochemistry, biomedical devices, and the encapsulation of biomolecules within silica gels [1–3]. The Si/SiO2 interface is encountered in the manufacturing of semiconductor transistors [4,5]. The O/SiO2 interface is encountered in the dry-etching and degradation of silica surfaces, commonly used in the aerospace industry, due to the collision by atomic oxygen [6,7]. The challenging characterization of these interfaces experimentally has stimulated a growing interest in using molecular dynamics (MD) simulations for that purpose. The MD simulations employ force fields for modeling atom– atom interactions in the materials of interest and the force fields are fitted to high-level ab initio or experimental results at specified temperatures and pressures [6,7]. Some force fields are rather simple, such as the Lennard-Jones pair potential [8], originally developed to model noble gases, while others are parameterized for modeling complex systems and interfaces. Examples are the twobody force fields with Coulomb interactions (e.g. BKS [6], Pedone ⇑ Corresponding author at: Institute for Space and Nuclear Power Studies, University of New Mexico, Albuquerque, NM, USA. http://dx.doi.org/10.1016/j.commatsci.2015.09.042 0927-0256/Ó 2015 Elsevier B.V. All rights reserved.

[7], TTAM [8], CHIK [9]), the three-body force fields with no Coulomb interactions (e.g. Munetoh [10]), and the bond-order, variable-charge (BOVC) force fields [11–13]. The BOVC force fields are the most suitable for modeling complex chemical reactions and interfaces. Two of the most popular BOVC force fields are the ChargeOptimized Many-Body (COMB) and ReaxFF. Since the introduction of the original COMB for simulating the Si/SiO2 interface [14], and the original ReaxFF for simulating hydrocarbons [15], several materials have been parameterized in forms of these powerful force fields. For instance, the capabilities with COMB have been successfully extended for modeling the interfaces of hafnium/ hafnium oxide [16] and Cu/SiO2 [17], as well as zirconium [18]. Similarly, the capabilities of ReaxFF have been extended for modeling the aluminum–molybdenum alloys [19], zinc oxide [20], and the silicon/silica system [21]. The COMB10 [11] force field is an extension of COMB07 [14], developed as an extension of the Tersoff-based potential by Yasukuawa [22] for semiconductors. The COMB10 force field, parameterized for modeling the Si/SiO2 interface and amorphous silica, accurately predicts the structural and elastic properties of crystalline SiO2 and the trend of defect formation in amorphous silica [11]. The ReaxFF is the foundation of two other parameterizations H2 O involving silica. The first is that of Fogarty et al. [12] (ReaxFFSiO ),

270

B.J. Cowen, M.S. El-Genk / Computational Materials Science 111 (2016) 269–276

who extended the parameterization by van Duin et al. [21] (ReaxFFSiO) to modeling the H2O/SiO2 interface. The ReaxFFSiO, based on the radical reactions described by quantum mechanical data, was unable to describe the proton transfer reactions at the interface of water with silica [12]. Fogarty et al. [12] replaced the O/H parameters in the ReaxFFSiO with a set fitted to the protontransfer reactions and the water clusters from two quantum mechanics datasets [24]. These datasets described the reaction energies for the Si(OH)4 polymerization, and the binding and dissociation of the water molecules from Si(OH)4 [23]. Thus, the Si/H, Si/ Si, and Si/O parameters were refitted [24]. Fogarty et al. [12] provide detailed information on the refitting schemes. The second parameterization of ReaxFF is that by Kulkarni et al. [16] ðReaxFFGSI SiO Þ. They extended ReaxFFSiO to modeling the gassurface interactions, specifically that of oxygen with silica. They carried out electronic structure calculations to obtain a highquality dataset for describing the geometry and the density functional theory energies, and used the new quantum–mechanical data to re-parameterize previous work [21]. The obtained parameters were optimized to reproduce the binding energy variations of the gas-surface interactions. Cowen and El-Genk [24] have reviewed and investigated four fixed-charge, two-body force fields (BKS [6], Pedone [7], TTAM [8], CHIK [9]) and one no-charge, three-body (Munetoh [10]) force field for modeling silica, and compared the results of the computational cost, accuracy, and transferability. The results of the MD simulations with the BKS force field accurately predicted the experimental values of the density, cell volume, and radial and bond-angle distribution functions for quartz, cristobalite, coesite, and stishovite, to within 2%. In addition, the results of the a–b transition and the quartz I–II transition were in good agreement with the experimental results [24]. The results of the simulations with both the CHIK and TTAM force fields, over the entire range of the silica phase diagram, were slightly less accurate than with the BKS potential. The simulations with the Munetoh force field were the cheapest in terms of the short computational time, but also the least accurate [24]. The predictions of the structural properties sometimes differed from experimental values by more than 10%. The simulations with the Munetoh force field did not produce the a–b and I–II phase transitions or the equation of state for stishovite silica. The results of the MD simulations with the Pedone potential of the material properties were within 5% of experimental values, except for the high-pressure polymorph of stishovite. The properties of stishovite were overestimated by 8% and the quartz I–II phase transition was overestimated by as much as 60% [24]. This work extends the previous effort by Cowen & El-Genk [24] to investigate the applicability of three variable-charge force fields H2 O ReaxFFSiO

ReaxFFGSI SiO

(COMB10 [11], [12], [13]) for modeling crystalline silica, developed specifically for characterizing complex interfaces and reactions. Investigated are the accuracy and transferability of these force fields for calculating the material properties, including the lattice constants, density, and the pair and bond-angle distribution functions. The transferability is examined by comparing the simulation results to the experimental values for the various silica polymorphs; namely, quartz, cristobalite, coeH2 O site, and stishovite. The COMB10 [11] and ReaxFFSiO [12] are also investigated for predicting the a–b phase transition for quartz. The MD simulations with the BOVC force fields require, as an input, the frequency of updating the atom charges. Generally, it is beneficial to update the atom charges as frequently as possible, without sacrificing either accuracy or computational efficiency. Updating the atom charges every time step, although computationally expensive, ensures that the charges are at their energy minima. Updating less frequently decreases the computational time,

but compromises the energy conservation. This work examined the effect of changing the frequency of updating the atom charges on the calculated results. Two updating options are considered: every time step, and every 10-time steps. The number that follows the abbreviation of the force field indicates the updating frequency (e.g., COMB10-1 refers to COMB10, with atom charges updated every time step). 2. Methodology The crystallographic information files (CIF) for all four polymorphs of silica were obtained from the Crystallography Open Database [25–27]. The silica polymorphs are replicated in all three spatial dimensions using symmetry operations to create a superGSI 2O cell. The MD simulations with COMB10, ReaxFFH SiO , and ReaxFFSiO ,

are carried out using the Large-Scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) code [28], with a time step of 0.2 fs. Aktulga et al. [29] implemented the ReaxFF parameterization into LAMMPS. The performed MD simulations track the positions and velocities of the atoms by time-integrating the Nose–Hoover style, non-Hamiltonian equations of motion, specifically by sampling from the isothermal–isobaric (NPT) ensemble [28]. The QEq [30] scheme is used to calculate the charge of each atom, based on the local bonding environment. The MD simulations of quartz, cristobalite, coesite, and stishovite used 5184, 768, 768, and 864 atoms, respectively. The energy minimization is carried out using the Polak–Ribiere version of the conjugate gradient (CG) algorithm [31], as discussed by Cowen and El-Genk [24]. The simulated systems are equilibrated for 100 ps and the results of the lattice constants, the average charge, and density, are collected over the following 100 ps. The ranges of charges are obtained by finding the minimum and maximum charges of the atoms following data collection. 2.1. Force fields The COMB10 [11] and ReaxFF [12,13] force fields are expressed somewhat differently. The COMB10 [11] force field is given as:

ET ¼ Eself þ ECoul þ EvdW þ Ebond þ Eothers

ð1Þ

The ReaxFF [12,13] is expressed as:

ET ¼ Eself þ ECoul þ EvdW þ Ebond þ Eangle þ Etorsion þ Econjugation þ EHbond þ Elonepair þ Eover þ Eunder þ Eothers

ð2Þ

The first four terms in Eq. (2) are basically the same for the COMB10 [11] force field, Eq. (1). Other terms: Eangle ; Etorsion ; Eover ; Eunder ; Econjugation

in Eq. (2) are incorporated

implicitly in the Ebond term of the COMB10 force field (1). The Eself term represents the energy difference in the various states of the atom charges, while Eothers consists of several terms that correct for allene-type terminal triple bond, and C2 species [32]. For more details on the physics of the various terms in Eqs. (1) and (2), the reader is referred to the original papers on the BOVC force fields [11–13]. 2.2. Quartz a–b transition In the performed MD simulations of the a–b transition for H2 O H2 O quartz using COMB10-1, ReaxFFSiO -1, and ReaxFFSiO -10, the temperature in the ranges of 300–700 K and 900–1500 K is increased in increments of 100 K. In the intermediate range of 700–900 K, the temperature increased in small increments of 10 K. The

271

B.J. Cowen, M.S. El-Genk / Computational Materials Science 111 (2016) 269–276

COMB10-10 force field required a slightly different approach, because the initial simulations predicted a structural transformation outside the 700–900 K range. Thus, the range for using the 10 K increments is shifted to higher temperatures to accurately capture the quartz a–b transition. In the simulations with the COMB10-10 force field, the temperature in the ranges of 300– 1100 K and 1400–1500 K is increased in increments of 100 K. Increments of 10 K are used in the intermediate range of 1100– 1400 K. The used heating time between temperatures is 1.0 ps, whereas equilibration times of 1, 2, 4, 8, and 16 ps and equal times for data collection are used. For example, a 16-ps of data collection time followed a 16-ps equilibration time at each temperature. 3. Results and discussion This section compares the calculated properties of the crystalline silica polymorphs of quartz, cristobalite, coesite, and stisho-

distribution functions. The presented results demonstrate the effect of changing the frequency of updating the atom charges every time step, and every 10-time steps on the calculated materials properties. Only the MD simulations with the COMB10 and H2 O ReaxFFSiO force fields are investigated for predicting the quartz a–b transition.

3.1. Quartz Quartz, the most abundant and stable SiO2 polymorph at ambient conditions, is used in many applications due to its superior optical and piezoelectric properties and the formation of SiO2 at the SiO2/Si interface [33–35]. Table 1 compares the calculated properties (standard deviation shown in parenthesis) and the atom charges for a-quartz, using MD simulations with the BOVC force fields to experimental values. The simulations results with the ReaxFFGSI SiO -1 accurately predict the lattice constants and density of

H2 O vite, using MD simulations with COMB10, ReaxFFSiO ; and

a-quartz, to within 2.1% of experimental values. The calculated

ReaxFFGSI SiO , with experimental results. These material properties are the lattice constants, density, and the pair and bond-angle

H2 O material properties in the simulations with the ReaxFFSiO -1 are within 3.6%, and those with COMB10 are within 7.3%, of the

Table 1 Comparison of calculated properties (standard deviation shown in parenthesis) and atom charges for a-quartz, using MD simulations with the BOVC force fields to experimental values. Parameter

a (Å) c (Å) q (g/cm3) qavg (O) qavg (Si) qrange (O) qrange (Si)

Experimental/force field Exp.[36]

COMB10-1

COMB10-10

H2 O -1 ReaxFFSiO

H2 O ReaxFFSiO -10

ReaxFFGSI SiO -1

4.916 5.405 2.646 – – – –

4.79 (4e3) 5.30 (6e3) 2.84 (3e3) 1.45 (1e4) 2.91 (3e4) 1.38 to 1.54 2.88–2.93

4.79 (3e3) 5.28 (5e3) 2.85 (2e3) 1.45 (5e5) 2.91 (1e4) 1.39 to 1.53 2.89–2.92

4.95 (3e3) 5.52 (3e3) 2.55 (1e3) 0.81 (2e3) 1.62 (4e3) 0.69 to 0.91 1.40–1.80

4.93 (1e2) 5.48 (3e2) 2.60 (3e3) 0.81 (1e3) 1.62 (4e3) 0.69 to 0.91 1.58–1.64

4.95 (3e3) 5.46 (3e3) 2.59 (1e3) 0.46 (9e4) 0.92 (2e3) 0.41 to 0.50 0.88–0.96

H2 O Fig. 1. Comparison of calculated pair distribution functions for a-quartz, at ambient conditions (300 K and 1 atm), in MD simulations with COMB10 and ReaxFFSiO to experimental results.

272

B.J. Cowen, M.S. El-Genk / Computational Materials Science 111 (2016) 269–276

experimental values. In the MD simulations of a-quartz with the COMB10 force field, the charge-equilibration frequency negligibly 2O affects the results. In the simulations with the ReaxFFH SiO , there is a larger disparity in the results, predicting a-quartz densities of 2.55 g/cm3 and 2.60 g/cm3 with a charge equilibration frequency of every time step and every 10-time steps, respectively. Fig. 1 compares the calculated pair distribution functions for a-quartz at ambient conditions, in the MD simulations with the

H2 O H2 O COMB10-1, COMB10-10, ReaxFFSiO -1, and ReaxFFSiO -10 force fields, with experimental values. The pair distribution functions are important for showing the structural order of crystalline materials. There are differences in the calculated distributions with updating the atom charges every time step and every 10-time steps. The distribution functions’ curves, produced by evaluating the atom charges every time step are smoother and more accurate, compared to experimental values. The most accurate curves for the O–Si and O–O pairs are those calculated in the simulations with the COMB10-1 force field. They are accurate to within 1% of the

H2 O the calculated curves with ReaxFFSiO -10, which disappear when updating the atom charges every time step. The MD simulations with the COMB10 force field underestimate the experimental results of the Si–O–Si BAD by 4%, while the simulations with H2 O ReaxFFSiO over-predict the experimental results by as much as 10% (Fig. 2). For judicious comparison to experimental results, Table 2 lists the ranges of the calculated bond-angles at the base and half-width, the average values of the bond-angles, and the bond-angles at the peaks of the BAD functions. The results in this table show close agreements with the experimental values of the average bond angles. The equilibration frequency of the atom charges clearly affects the simulations results, as demonstrated in Figs. 1 and 2. In addition to not conserving energy, when updating the atom charges each 10 time steps, the peaks of the calculated pair and

H2 O experimental results. The simulation results with ReaxFFSiO -1 for the Si–Si pairs are the most accurate; they are within 1% of experimental values. The MD simulations with the COMB10 force field underestimate the Si–Si distribution by 3%, and that with the H2 O ReaxFFSiO underestimate both the O–Si and O–O pairs by 3%, compared to the experimental values of the pair distribution functions (Fig. 1). For an accurate characterization, Table 2 lists the ranges of the calculated bond-lengths at the base and half-width, the average values of the bond-lengths, and the bond-lengths at the peaks of the pair distribution functions. Fig. 2 compares the calculated O–Si–O and Si–O–Si bond-angle distribution (BAD) functions for a-quartz at ambient conditions, using MD simulations with the COMB10-1, COMB10-10, H2 O H2 O ReaxFFSiO -1, and ReaxFFSiO -10 to the experimental values. The results with all BOVC force fields predict the O–Si–O BAD to within 1% of the experimental values. However, there are slight kinks in

Fig. 2. Comparison of calculated BAD functions for a-quartz, at ambient conditions H2 O (300 K and 1 atm), in MD simulations with COMB10 and ReaxFFSiO to experimental results.

Table 2 Comparison of the bond-lengths at the 1st peak of the pair distribution function, and the bond-angles for a-quartz, using MD simulations with the BOVC force fields, to experimental values. Force field

Atom Pairs

Bond lengths for 1st peak in RDF (Å) Min–Max @ Base

Min–Max @ Half-Width

@ Peak

Average value Calculated

Exp. [36]

COMB10-1

O–Si O–O Si–Si

1.46–1.86 2.29–2.87 2.76–3.41

1.57–1.64 2.54–2.72 2.95–3.05

1.61 2.63 2.98

1.62 2.63 3.01

1.61 2.61–2.64 3.06

COMB10-10

O–Si O–O Si–Si

1.51–1.80 2.33–2.82 2.85–3.38

1.59–1.64 2.56–2.70 2.95–3.02

1.61 2.61 2.99

1.62 2.62 3.02

1.61 2.61–2.64 3.06

H2 O ReaxFFSiO -1

O–Si O–O Si–Si

1.46–1.78 2.19–2.94 2.94–3.31

1.53–1.61 2.47–2.66 3.04–3.13

1.57 2.55 3.08

1.57 2.56 3.09

1.61 2.61–2.64 3.06

H2 O ReaxFFSiO -10

O–Si O–O Si–Si

1.52–1.65 2.16–2.90 3.02–3.23

1.56–1.59 2.44–2.67 3.06–3.10

1.57 2.52 3.08

1.57 2.57 3.09

1.61 2.61–2.64 3.06

Interface

Bond angles (°)

COMB10-1

O–Si–O Si–O–Si

89.9–130.5 99.2–153.5

103.6–115.2 133.6–141.5

108.9 136.7

109.3 137.7

109.5 143.7

COMB10-10

O–Si–O Si–O–Si

94.4–128.5 128.6–155.3

104.4–115.1 134.1–140.1

108.7 137.2

109.3 138.0

109.5 143.7

H2 O ReaxFFSiO -1

O–Si–O Si–O–Si

91.2–132.6 137.5–179.6

103.6–114.6 151.9–163.5

108.3 156.6

109.4 158.0

109.5 143.7

H2 O ReaxFFSiO -10

O–Si–O Si–O–Si

89.6–134.2 146.6–172.0

102.6–115.5 154.2–160.8

107.0 157.0

109.3 157.7

109.5 143.7

B.J. Cowen, M.S. El-Genk / Computational Materials Science 111 (2016) 269–276

bond-angle distribution functions are not smooth, unlike those with updating the atom charges each time step. Notwithstanding that, the average values are still remarkably close to the experimental results for both charge-equilibration frequencies. At ambient conditions, the impact of the charge-equilibration frequency on the MD simulation results is small; however, predicting the complex phase transition is a challenge. 3.1.1. a–b transition The a–b phase transition for quartz has been shown experimentally to occur at 846 K [37]. A good test of the transferability of a force field in MD simulations is to characterize this transition. This section compares the calculated cell volume and c/a ratio at the a– b transition in MD simulations with the different BOVC force fields, as a function of temperature, to reported experimental values. 3.1.1.1. Cell volume. Fig. 3 shows that in the performed MD simulations with the COMB10-1 force field for durations of 1, 2, and 4 ps, the calculated cell volume steadily increases with increasing temperature, with no structural transformation. The simulation results with the COMB10-1 force field and with equilibration and data collection times of 8 ps and 16 ps show a sharp drop in the calculated cell volume at the a–b transition. Conversely, the experimental results indicate a plateau in the cell volume expansion at the a–b transition, due to the rotation of the SiO4 tetrahedra. The MD simulations with the COMB10-10 force field (Fig. 3) and different equilibration times reveal a structural transformation, but at a much higher temperature than predicted with the COMB10-1 potential. As with the COMB10-1 potential, the not fully converged results show a continuous decrease in the cell volume at the a–b transition, rather than reaching a plateau, as indicated experimentally (Fig. 3). 2O ReaxFFH SiO -1

The MD simulation results with the converged for all cases (Fig. 3). However, a plateau is not predicted over the entire temperature range (300–1500 K) of interest, indicating the inabilH2 O ity of the ReaxFFSiO -1 to model the a–b transition for quartz. Only

273

the simulation results of the cell volume expansion using the H2 O ReaxFFSiO -10, with 16 ps of equilibration and data collection time, predict a structural transformation. However, the transformation does not predict the experimentally measured plateau at the a–b transition (Fig. 3). The simulations with the BOVC force fields and with updating the atom charges every 10 time steps give drastically different results than when updating the charges every time step. Indeed, the period of the vibrational modes for the high frequency asymmetric Si–O bond stretching is approximately an order of magnitude larger than the 2 fs time interval for the lower chargeequilibration frequency [38]. However, the charge-equilibration frequency is uncoupled to the vibrational modes. Thus, in the MD simulations with the force fields of COMB and ReaxFF, the atom charges need to be updated every time step in order to find the local energy minima. As indicated earlier, updating the atom charges less frequently no longer conserves energy. At ambient conditions, the effect of the charge-equilibration frequency on the calculated results is small, but more pronounced over the quartz a–b transition, as the results using chargeequilibration frequencies of each time step and each 10 time steps diverge. Therefore, for the complex a–b structural phase transition, it is important to update atom charges as often as possible to compensate for the change in the bonding environment. Perhaps there is a more optimal value between the two investigated in this work, namely each time step and 10-time steps. Even when the atom charges are updated every time step, the MD simulations results H2 O with the COMB10 and the ReaxFFSiO inaccurately predict the quartz a–b transition.

3.1.1.2. c/a. The c/a ratio at the a–b transition for quartz is difficult to model and has eluded many force fields. Herzbach et al. [39] indicated that the polarizable TS [40] is the only force field, among those they investigated, that is capable of modeling the sudden drop in c/a ratio at the a–b transition for quartz. The present MD simulation results with the COMB10-1 and COMB10-10 force fields

H2 O Fig. 3. Comparison of calculated cell volume expansion at the a–b transition of quartz in MD simulations with COMB10 and ReaxFFSiO to experimental results.

274

B.J. Cowen, M.S. El-Genk / Computational Materials Science 111 (2016) 269–276

H2 O Fig. 4. Comparison of calculated c/a ratio at the a–b transition of quartz, in MD simulations with COMB10 and ReaxFFSiO to experimental results.

inaccurately predict the c/a ratio at the a–b transition (Fig. 4). At the a–b transition, the calculated c/a ratio experiences a sudden increase, rather than a decrease, as obtained experimentally. The

a-quartz of 2.84 g/cm3. Thus, it is unclear how reported work

H2 O performed simulations with the ReaxFFSiO -10 show a similar

H2 O simulations with both ReaxFFSiO and ReaxFFGSI SiO are within 3% and 5% of the lattice constants and density, respectively, compared to experimental values (Table 3). These results demonstrate the transferability of these force fields to the a-cristobalite bonding environment.

increase in the c/a ratio near the a–b transition. The is the only force field that does not show any structural transformation in the calculated c/a ratio. H2 O ReaxFFSiO -1

[11] could have predicted values closer to experimental results. The calculated properties for a-cristobalite in the performed MD

3.2. Cristobalite

3.3. Coesite

Similar to quartz, cristobalite is formed at higher pressures, and is a four-coordinated polymorph. However, the density is lower, comparable to that of amorphous silica. Based on the previous results showing that evaluating the atom charges every 10 time steps is not quite accurate and does not conserve energy (Figs. 1–4), the results of the performed MD simulations discussed in remainder of this paper, are with updating the atom charges every time step. In Table 3, the calculated structural properties and charges of the oxygen and silicon atoms of a-cristobalite, in MD simulations with the three BOVC force fields investigated in this work, are compared to experimental values. Despite the fact that a-quartz and acristobalite have similar crystal structures, with the silicon atoms tetrahedrally coordinated, the results of the MD simulations with the COMB10 force field of the material properties are not as accurate for this polymorph, as they are for a-quartz (Table 3). Despite the extensive simulations with different time steps, frequencies of updating the atom charges of the atoms, equilibration and data collection times, and number of atoms, the calculated results for this polymorph are significantly different from those previously reported with the COMB10 force field [11]. Energy minimization of the force field produces results consistent with the experimental values, but equilibration at ambient conditions produces a 21% higher density. The calculated density in this work for the a-cristobalite of 2.81 g/cm3 is similar to that calculated for

Coesite is another four-coordinated silica polymorph, which forms at high pressures due to events such as meteor impacts [42]. The bonding environment is thus much different from those for quartz and cristobalite. The monoclinic crystal structure of coesite comprises four silicon dioxidetetrahedra in Si4O8 and Si8O16 rings [43]. Table 4 compares the calculated properties and the charges of the oxygen and silicon atoms for coesite with experimental values. The results of the performed MD simulations using all three BOVC force fields accurately predict the properties of coesite. The monoclinic cell of coesite is characterized as a – b – c and b – 90°. The b angle is shown experimentally to be 120.34° [43]. All three force fields predict this angle to within 1° of the experimental value. The MD simulation results of the coesite properties, H2 O with the COMB10, ReaxFFSiO , and ReaxFFGSI SiO force fields, are within 1.7%, 2.2%, and 2.3% of the experimental values, respectively. Such good accuracy attests to the transferability of these force fields to the coesite bonding environment.

3.4. Stishovite The crystal structure of stishovite resembles that of rutile (TiO2), with each of the octahedrally coordinated silicon atoms bound to six oxygen atoms [44]. Therefore, the transferability of

275

B.J. Cowen, M.S. El-Genk / Computational Materials Science 111 (2016) 269–276

Table 3 Comparison of calculated properties (standard deviation shown in parenthesis) and atom charges for a-cristobalite in MD simulations with the BOVC force fields, to experimental values. Parameter

a (Å) c (Å) q (g/cm3) qavg (O) qavg (Si) qrange (O) qrange (Si)

Experimental/force field Exp. [41]

COMB10 [11]

H2 O [12] ReaxFFSiO

ReaxFFGSI SiO [13]

4.978 6.948 2.318 – – – –

4.81 (1e2) 6.13 (1e2) 2.81 (7e3) 1.45 (1e4) 2.90 (3e4) 1.39 to 1.51 2.88–2.92

5.03 (7e3) 7.11 (1e2) 2.22 (3e3) 0.87 (4e3) 1.73 (7e3) 0.76 to 0.97 1.58–1.90

5.04 (9e3) 7.12 (1e2) 2.20 (4e3) 0.47 (1e3) 0.94 (3e3) 0.43 to 0.51 0.90–0.98

Table 4 Comparison of calculated properties (standard deviation shown in parenthesis) and atom charges for coesite, in MD simulations with the BOVC force fields, to experimental values. Parameter

a (Å) b (Å) c (Å) q (g/cm3) b (°) qavg (O) qavg (Si) qrange (O) qrange (Si)

Exp. [43]

7.136 12.369 7.174 2.921 120.34 – – – –

Force fields COMB10 [11]

H2 O [12] ReaxFFSiO

ReaxFFGSI SiO [13]

7.04 (1e2) 12.40 (2e2) 7.22 (9e3) 2.97 (6e3) 121.12 (1e1) 1.46 (2e4) 2.91 (4e4) 1.37 to 1.54 2.89–2.93

7.29 (2e2) 12.33 (1e2) 7.12 (7e3) 2.88 (3e3) 119.75 (8e2) 0.79 (3e3) 1.59 (5e3) 0.71 to 0.88 1.48–1.76

7.20 (1e2) 12.19 (2e2) 7.01 (1e2) 2.96 (4e3) 119.41 (1e1) 0.46 (1e3) 0.91 (2e3) 0.41 to 0.50 0.87–0.95

Table 5 Comparison of calculated properties (standard deviation shown in parenthesis) and atom charges for stishovite, in MD simulations with BOVC force fields to experimental values. Parameter

a (Å) c (Å) q (g/cm3) qavg (O) qavg (Si) qrange (O) qrange (Si)

Exp. [44]

4.180 2.666 4.283 – – – –

Force field COMB10 [11]

H2 O [12] ReaxFFSiO

ReaxFFGSI SiO [13]

4.09 (7e3) 2.97 (2e3) 4.02 (4e3) 1.47 (1e4) 2.94 (3e4) 1.43 to 1.51 2.93–2.96

4.56 (4e3) 3.04 (2e3) 3.16 (3e3) 0.53 (2e3) 1.07 (3e3) 0.48 to 0.59 1.02–1.11

5.68 (4e2) 2.54 (2e2) 2.38 (2e2) 0.41 (3e3) 0.82 (6e3) 0.23 to 0.51 0.60–0.96

the force fields can be critically evaluated based on their ability to change from a tetrahedra bonding environment (quartz, cristobalite, coesite) to an octahedra bonding environment (stishovite). Table 5 compares the calculated properties and charges of the oxygen and silicon atoms of stishovite in the MD simulations with the three BOVC force fields, with experimental values. H2 O The performed MD simulations with ReaxFFSiO and ReaxFFGSI SiO inaccurately predict the properties of the stishovite under equilibration in the NPT ensemble. The calculated properties are within H2 O 26% and 44% of experimental results for the ReaxFFSiO and

ReaxFFGSI SiO , respectively. This could be because the density and the cohesive energies of various polymorphs were determined with single point energies at different densities, which is more of an evaluation than an optimization. Thus, despite the fact that MD H2 O simulations with both ReaxFFSiO and ReaxFFGSI SiO accurately predict the properties of the other silica polymorphs, these force fields are not necessarily suitable for all phases of silica, and need to be optimized more effectively for stishovite. The MD simulation results with the COMB10 force field are accurate, to within 11% of the experimental values for the ‘‘c” lattice parameter, and 2.2% for the ‘‘a” lattice parameter (Table 5). Thus, the COMB10 force field is a more appropriate choice for modeling stishovite.

4. Summary & conclusions The accuracy and transferability of the three BOVC force fields: H2 O COMB10, ReaxFFSiO , and ReaxFFGSI SiO , are investigated for modeling the four SiO2 polymorphs, quartz, cristobalite, coesite, and stishovite. These force fields have been designed for modeling interfaces and complex chemical reactions. MD simulation results of the structural properties of crystalline silica, using these force fields, are compared to experimental values. In these simulations, the atom charges are updated each time step. The results with updating the atom charges every 10-time steps are inaccurate, since energy is not conserved. The MD simulation results with the COMB10 force field accurately predict the properties of the quartz and coesite polymorphs at ambient conditions, to within 7.3% and 1.7% of the experimental values, respectively. The simulation results for the a-cristobalite and stishovite are not as accurate; they are within 21% and 11% of the experimental values, respectively. H2 O The simulation results with ReaxFFSiO and ReaxFFGSI SiO , of the properties for quartz, cristobalite, and coesite at ambient conditions, except the BAD function, are accurate to within 5% of the H2 O experimental values. The simulation results with ReaxFFSiO of

276

B.J. Cowen, M.S. El-Genk / Computational Materials Science 111 (2016) 269–276

the BAD function for quartz are within 10% of the experimental values. H2 O The simulation results with ReaxFFSiO and ReaxFFGSI SiO inaccurately predict the properties of the high-pressure stishovite polymorph; they are within 26% and 44% of the experimental values, respectively. In addition, the results of the MD simulations H2 O with the COMB10 and the ReaxFFSiO inaccurately model the a–b transition for quartz. The results of the cell volume expansion at the a–b transition are quite different from experimental values and do not predict the experimentally measured plateau, caused by the rotation of the SiO4 tetrahedra.

Acknowledgements The financial support for this research is provided by the Institute for Space and Nuclear Power Studies at the University of New Mexico (ISNPS-UNM) and by a U.S. Nuclear Regulatory Commission (NRC) Graduate Fellowship under Grant # NRC-38-09-931 to ISNPS-UNM. The authors thank Dr. Susan Atlas of the Center for Advanced Research Computing (CARC), at the University of New Mexico for her help and suggestions, acknowledge the valuable support provided by the CARC staff and are grateful for using the supercomputer clusters at the center. References [1] J. Livage, T. Coradin, C. Roux, J. Phys.-Condens. Matter 13 (2001) R673. [2] A. Cimas, F. Tielens, M. Sulpizi, M.P. Gaigeot, D. Costa, J. Phys.-Condens. Matter 26 (2014) 244106. [3] A. Rimola, D. Costa, M. Sodupe, J.F. Lambert, P. Ugliengo, Chem. Rev. 113 (2013) 4216. [4] M. Schulz, Nature 399 (1999) 729. [5] C.J. Cochrane, P.M. Lenahan, Appl. Phys. Lett. 103 (2013) 053506. [6] B.W. van Beest, G.J. Kramer, R.A. van Santen, Phys. Rev. Lett. 64 (1990) 1955. [7] A. Pedone, G. Malavasi, M.C. Menziani, A.N. Cormack, U. Segre, J. Phys. Chem. B 110 (2006) 11780. [8] S. Tsuneyuki, M. Tsukada, H. Aoki, Y. Matsui, Phys. Rev. Lett. 61 (1988) 869. [9] A. Carre, J. Horbach, S. Ispas, W. Kob, EPL-Europhys. Lett. 82 (2008) 17001. [10] S. Munetoh, T. Motooka, K. Moriguchi, A. Shintani, Comput. Mater. Sci. 39 (2007) 334.

[11] T.R. Shan, B.D. Devine, J.M. Hawkins, A. Asthagiri, S.R. Phillpot, S.B. Sinnott, Phys. Rev. B 82 (2010) 255. [12] J.C. Fogarty, H.M. Aktulga, A.Y. Grama, A.C.T. van Duin, S.A. Pandit, J. Chem. Phys. 132 (2010) 174704. [13] A.D. Kulkarni, D.G. Truhlar, S.G. Srinivasan, A.C.T. van Duin, P. Norman, T.E. Schwartzentruber, J. Phys. Chem. C 117 (2013) 258. [14] J.G. Yu, S.B. Sinnott, S.R. Phillpot, Phys. Rev. B 75 (2007) 085311. [15] A.C.T. van Duin, S. Dasgupta, F. Lorant, W.A. Goddard, J. Phys. Chem. A 105 (2001) 9396. [16] T.R. Shan, B.D. Devine, T.W. Kemper, S.B. Sinnott, S.R. Phillpot, Phys. Rev. B 81 (2010) 125328. [17] T.R. Shan, B.D. Devine, S.R. Phillpot, S.B. Sinnott, Phys. Rev. B 83 (2011) 115327. [18] M.J. Noordhoek, T. Liang, Z.Z. Lu, T.R. Shan, S.B. Sinnott, S.R. Phillpot, J. Nucl. Mater. 441 (2013) 274. [19] W.X. Song, S.J. Zhao, J. Mater. Res. 28 (2013) 1155. [20] D. Raymand, A.C.T. van Duin, M. Baudin, K. Hermansson, Surf. Sci. 602 (2008) 1020. [21] A.C.T. van Duin, A. Strachan, S. Stewman, Q.S. Zhang, X. Xu, W.A. Goddard, J. Phys. Chem. A 107 (2003) 3803. [22] A. Yasukawa, JSME Int. J. A-Mech. M 39 (1996) 313. [23] T.T. Trinh, A.P.J. Jansen, R.A. van Santen, J. Phys. Chem. B 110 (2006) 23099. [24] B.J. Cowen, M.S. El-Genk, Comput. Mater. Sci. 107 (2015) 88. [25] S. Grazulis, D. Chateigner, R.T. Downs, A.F.T. Yokochi, M. Quiros, L. Lutterotti, E. Manakova, J. Butkus, P. Moeck, A. Le Bail, J. Appl. Crystallogr. 42 (2009) 726. [26] S. Grazulis, A. Daskevic, A. Merkys, D. Chateigner, L. Lutterotti, M. Quiros, N.R. Serebryanaya, P. Moeck, R.T. Downs, A. Le Bail, Nucleic Acids Res. 40 (2012) D420. [27] R.T. Downs, M. Hall-Wallace, Am. Mineral. 88 (2003) 247. [28] S. Plimpton, J. Comput. Phys. 117 (1995) 1. [29] H.M. Aktulga, J.C. Fogarty, S.A. Pandit, A.Y. Grama, Parallel Comput. 38 (2012) 245. [30] A.K. Rappé, W.A. Goddard, J. Phys. Chem.-Us 95 (1991) 3358. [31] E. Polak, G. Ribiere, Rev. Fr. Inform. Rech. Oper. 3 (1969) 35. [32] Y.K. Shin, T.R. Shan, T. Liang, M.J. Noordhoek, S.B. Sinnott, A.C.T. van Duin, S.R. Phillpot, MRS Bull. 37 (2012) 504. [33] I. Tajima, O. Asami, E. Sugiura, Anal. Chim. Acta 365 (1998) 147. [34] J. Fontanel, C. Andeen, D. Schuele, J. Appl. Phys. 45 (1974) 2852. [35] E.K. Chang, M. Rohlfing, S.G. Louie, Phys. Rev. Lett. 85 (2000) 2613. [36] L. Levien, C.T. Prewitt, D.J. Weidner, Am. Mineral. 65 (1980) 920. [37] M.A. Carpenter, E.K.H. Salje, A. Graeme-Barber, B. Wruck, M.T. Dove, K.S. Knight, Am. Mineral. 83 (1998) 2. [38] Y.H. Kim, M.S. Hwang, H.J. Kim, J.Y. Kim, Y. Lee, J. Appl. Phys. 90 (2001) 3367. [39] D. Herzbach, K. Binder, M.H. Müser, J. Chem. Phys. 123 (2005) 124711. [40] P. Tangney, S. Scandolo, J. Chem. Phys. 117 (2002) 8898. [41] D.R. Peacor, Z. Kristallogr. 138 (1973) 274. [42] E.C.T. Chao, E.M. Shoemaker, B.M. Madsen, Science 132 (1960) 220. [43] L. Levien, C.T. Prewitt, Am. Mineral. 66 (1981) 324. [44] L.G. Liu, W.A. Bassett, T. Takahashi, J. Geophys. Res. 79 (1974) 1160.