^{1}and Betsy M. Rice

^{1}

### Abstract

A new short-range pairwise numerical potential for silica is presented. The potential is derived from a single *ab initio* molecular dynamics (AIMD) simulation of molten silica using the force-matching method with the forces being represented numerically by piecewise functions (splines). The AIMD simulation is performed using the Born-Oppenheimer method with the generalized gradient approximation (BLYP) for the XC energy functional. The new effective potential includes a soft-repulsive shoulder to describe the interactions of oxygen ions at short separations. The new potential, despite being short-ranged and derived from single-phase data, exhibits a good transferability to silica crystalline polymorphs and amorphous silica. The importance of the O–O soft-repulsive shoulder interaction on glass densification under cold and shock compressions is assessed from MD simulations of silicaglass under room and shock Hugoniot conditions, respectively. Results from these simulations indicate that the appearance of oxygen complexes (primarily pairs) interacting through soft-repulsive shoulder potential occurs at 8–10 GPa, and under cold compression conditions becomes notable at 40 GPa, essentially coinciding with the transition to a Si sixfold coordination state. An analysis of changes in system structure in compressed and shocked states reveals that the O ions interacting through the soft-repulsive shoulder potential in denser states of silicaglass may create a mechanical multi-stability under elevated pressures and thus to contribute to the observed anomalous densification.

The authors wish to thank Dr. George Gazonas, Dr. Iskander Batyrev, and Dr. Scott Weingarten for helpful comments. This research was supported by the DoD High Performance Computing Modernization Program Software Application Institute for Multi-scale Reactive Modeling of Insensitive Munitions and the ARL Director under the Director's Strategic Initiative program “Multi-scale Modeling of Non-Crystalline Ceramics (Glass).” Computing support was provided by the DoD Supercomputer Resource Center located at the Army Research Laboratory.

I. INTRODUCTION

II. POTENTIAL DEVELOPMENT WITH THE MS-CG METHOD

III. FORCE-MATCHING MODELS

IV. TRANSFERABILITY OF THE FM MODEL TO SILICAPOLYMORPHS

V. GLASS UNDER PRESSURE

VI. CONCLUSIONS

### Key Topics

- Silica
- 52.0
- High pressure
- 20.0
- Mercury (element)
- 18.0
- Polymorphism
- 18.0
- Ab initio calculations
- 16.0

## Figures

Effective atom-atom forces [panel (a)] and corresponding potentials [panel (b)] in liquid SiO_{2} generated through the FM method as functions of interatomic separation: O–O (black), Si–O (red), and Si–Si (green). The dashed line corresponds to the model without the repulsive shoulder (FM-ns model). In panel (b) the dotted (FM-ws model) and dotted-dotted-dashed lines indicate the variation of the O–O repulsion in FM with block averaging along reference trajectories as discussed in the text. In panel (c) a comparison of the BKS (solid), Pedone (solid/circles), and FM (dashed) potentials is given. Due to the aphysical behavior of the BKS exp-6 at small interatomic separations, we modified the function to generate the BKS curve shown in this figure as follows: For *r* < *r* _{ infl }, where *r* _{ infl } is the position of the inflection point [*u* ^{″}(*r* _{ infl }) = 0] on the repulsive wall of the original BKS exp-6 potential, the curve is described by an exponential function whose energy and energy first derivative are the same as the original BKS exp-6 function at *r* _{ infl }. The vertical dashed line marks the location of the repulsive shoulder in the O–O FM potential.

Effective atom-atom forces [panel (a)] and corresponding potentials [panel (b)] in liquid SiO_{2} generated through the FM method as functions of interatomic separation: O–O (black), Si–O (red), and Si–Si (green). The dashed line corresponds to the model without the repulsive shoulder (FM-ns model). In panel (b) the dotted (FM-ws model) and dotted-dotted-dashed lines indicate the variation of the O–O repulsion in FM with block averaging along reference trajectories as discussed in the text. In panel (c) a comparison of the BKS (solid), Pedone (solid/circles), and FM (dashed) potentials is given. Due to the aphysical behavior of the BKS exp-6 at small interatomic separations, we modified the function to generate the BKS curve shown in this figure as follows: For *r* < *r* _{ infl }, where *r* _{ infl } is the position of the inflection point [*u* ^{″}(*r* _{ infl }) = 0] on the repulsive wall of the original BKS exp-6 potential, the curve is described by an exponential function whose energy and energy first derivative are the same as the original BKS exp-6 function at *r* _{ infl }. The vertical dashed line marks the location of the repulsive shoulder in the O–O FM potential.

Panel (a): Errors Δ*F* vs. time. Instant values of Δ*F* are shown for oxygen ions (black) and silicon ions (cyan). White lines show corresponding running time averages 〈Δ*F*〉_{ t }. Panel (b): Distributions of errors *P*(Δ*F*) for oxygen (solid/thin), silicon (dashed), and total (solid/thick).

Panel (a): Errors Δ*F* vs. time. Instant values of Δ*F* are shown for oxygen ions (black) and silicon ions (cyan). White lines show corresponding running time averages 〈Δ*F*〉_{ t }. Panel (b): Distributions of errors *P*(Δ*F*) for oxygen (solid/thin), silicon (dashed), and total (solid/thick).

King's ring size distribution [panel (a)] and all-particle RDF, *g*(*r*), (solid) with corresponding running integration number, , (dashed) [panel (b)] for the (FM-l) (red), (FM-h) (blue) structures using the FM model and for the (P-l) structure using the Pedone model (green).

King's ring size distribution [panel (a)] and all-particle RDF, *g*(*r*), (solid) with corresponding running integration number, , (dashed) [panel (b)] for the (FM-l) (red), (FM-h) (blue) structures using the FM model and for the (P-l) structure using the Pedone model (green).

Density vs. pressure at *T* = 298 K (solid lines) and along the Hugoniot locus (dashed lines) from simulations of (FM-l) structure using the FM (cyan), FM-ns (green) [panel (a)], Pedone (blue) [panel (b)], and FM-ws (magenta) [panel (c)] models. Simulations were carried out at the same pressure points (given in the text) as marked by filled circles on the cyan line. Experimental EOS obtained by cold compression (black squares) and by shock compression (red circles) is from Refs. 6 and 65, respectively. Inset to panel (c) compares the 298 K (solid) and Hugoniot (dashed) EOS from simulations of the (FM-l) (cyan) and (FM-h) (black) structures using the FM model.

Density vs. pressure at *T* = 298 K (solid lines) and along the Hugoniot locus (dashed lines) from simulations of (FM-l) structure using the FM (cyan), FM-ns (green) [panel (a)], Pedone (blue) [panel (b)], and FM-ws (magenta) [panel (c)] models. Simulations were carried out at the same pressure points (given in the text) as marked by filled circles on the cyan line. Experimental EOS obtained by cold compression (black squares) and by shock compression (red circles) is from Refs. 6 and 65, respectively. Inset to panel (c) compares the 298 K (solid) and Hugoniot (dashed) EOS from simulations of the (FM-l) (cyan) and (FM-h) (black) structures using the FM model.

Temperature along the Hugoniot locus, *T* _{Hg}(*P*), simulated using the FM (red), FM-ns (green), and Pedone (blue) models at the same pressure points as in Fig. 4.

Temperature along the Hugoniot locus, *T* _{Hg}(*P*), simulated using the FM (red), FM-ns (green), and Pedone (blue) models at the same pressure points as in Fig. 4.

Pressure dependence of *g* _{X}(*r*) at *T* = 298 K [panel (a)] and along the Hugoniot locus [panel (b)] for the (FM-l) structure using the FM model. The inset is a magnification of the region of the first peak. Arrows show the direction of changes in peak positions with pressure.

Pressure dependence of *g* _{X}(*r*) at *T* = 298 K [panel (a)] and along the Hugoniot locus [panel (b)] for the (FM-l) structure using the FM model. The inset is a magnification of the region of the first peak. Arrows show the direction of changes in peak positions with pressure.

Panel (a): *T* = 298 K (solid) and shocked (dashed) O–O coordination number within the shell *r* < 0.153 nm, , for the (FM-l) structure using the FM model. The inset is an enlargement of the low pressure region. Panel (b): *T* = 298 K (solid) and shocked (dashed) Si–O coordination number, , from FM (red, filled squares), FM-ns (blue, empty squares), and Pedone (green, empty diamonds) models. Experimental cold compression data obtained by X-ray adsorption measurements (black, filled squares) are from Refs. 7 and 9.

Panel (a): *T* = 298 K (solid) and shocked (dashed) O–O coordination number within the shell *r* < 0.153 nm, , for the (FM-l) structure using the FM model. The inset is an enlargement of the low pressure region. Panel (b): *T* = 298 K (solid) and shocked (dashed) Si–O coordination number, , from FM (red, filled squares), FM-ns (blue, empty squares), and Pedone (green, empty diamonds) models. Experimental cold compression data obtained by X-ray adsorption measurements (black, filled squares) are from Refs. 7 and 9.

Comparison of the probability density distributions of O–Si–O angles at *T* = 298 K and different pressures by the FM model for the (FM-l) structure [panel (a)] and the Pedone model for structure (P-l) [panel (b)]. The panels (c) and (d) show similar comparison for Si–O–Si angle distributions.

Comparison of the probability density distributions of O–Si–O angles at *T* = 298 K and different pressures by the FM model for the (FM-l) structure [panel (a)] and the Pedone model for structure (P-l) [panel (b)]. The panels (c) and (d) show similar comparison for Si–O–Si angle distributions.

*g* _{X}(*r*) (solid) and *g* _{OO}(*r*) (dashed) from *T* = 298 K simulations at *P* = 74 GPa for the (FM-l) structure using FM (blue) and FM-ns (red) models.

*g* _{X}(*r*) (solid) and *g* _{OO}(*r*) (dashed) from *T* = 298 K simulations at *P* = 74 GPa for the (FM-l) structure using FM (blue) and FM-ns (red) models.

Compression ρ_{0}/ρ, where ρ_{0} is ambient density vs. pressure curves for the (FM-l) structure decompressed using FM (solid) and FM-ns (dashed) models, which was previously compressed with the FM model under cold (diamonds, blue) and shock (red, triangles) conditions. Densities in samples compressed under cold (solid black, empty circles) and shock (solid red, empty squares) conditions are shown for a comparison.

Compression ρ_{0}/ρ, where ρ_{0} is ambient density vs. pressure curves for the (FM-l) structure decompressed using FM (solid) and FM-ns (dashed) models, which was previously compressed with the FM model under cold (diamonds, blue) and shock (red, triangles) conditions. Densities in samples compressed under cold (solid black, empty circles) and shock (solid red, empty squares) conditions are shown for a comparison.

*g* _{X}(*r*) for the decompressed (FM-l) sample from selected simulations using the FM model from Fig. 10. Panel (a): For samples compressed at pressures 1, 2, 3, 5, and 8.118 GPa under cold conditions. Panel (b): For samples compressed at pressure 74 GPa under cold (solid green) and shock (dashed red) conditions. The cyan lines show the original *g* _{X}(*r*) in the sample compressed at 74 GPa under cold (solid cyan) and shock (dashed cyan) conditions. In both panels the ambient structure is shown by the black solid line.

*g* _{X}(*r*) for the decompressed (FM-l) sample from selected simulations using the FM model from Fig. 10. Panel (a): For samples compressed at pressures 1, 2, 3, 5, and 8.118 GPa under cold conditions. Panel (b): For samples compressed at pressure 74 GPa under cold (solid green) and shock (dashed red) conditions. The cyan lines show the original *g* _{X}(*r*) in the sample compressed at 74 GPa under cold (solid cyan) and shock (dashed cyan) conditions. In both panels the ambient structure is shown by the black solid line.

## Tables

Coefficients of the least-squares fit for the BLYP FM SiO_{2} force field *f* _{αβ}(*r*) using the expansion in Eq. (2) with *n* _{ max } = 16. Atomic units are used. At small separations *r* < *r* _{ core }, the *f* _{αβ}(*r*) is extrapolated as . The following core radii are used: a.u., a.u., and a.u. The cutoff of 0.7038 nm must be applied to this expansion. The original numerical forces and potentials are provided in the supplementary material.^{66}

Coefficients of the least-squares fit for the BLYP FM SiO_{2} force field *f* _{αβ}(*r*) using the expansion in Eq. (2) with *n* _{ max } = 16. Atomic units are used. At small separations *r* < *r* _{ core }, the *f* _{αβ}(*r*) is extrapolated as . The following core radii are used: a.u., a.u., and a.u. The cutoff of 0.7038 nm must be applied to this expansion. The original numerical forces and potentials are provided in the supplementary material.^{66}

Structural parameters, *a, b, c, β* (length in nm and angle in deg), density ρ (kg/m^{3}), cohesive energy per SiO_{2} *E* ^{ c } (eV), elastic constants *C* _{ ij } (GPa), and bulk modulus *B* (GPa), for silica polymorphs.

Structural parameters, *a, b, c, β* (length in nm and angle in deg), density ρ (kg/m^{3}), cohesive energy per SiO_{2} *E* ^{ c } (eV), elastic constants *C* _{ ij } (GPa), and bulk modulus *B* (GPa), for silica polymorphs.

Article metrics loading...

Full text loading...

Commenting has been disabled for this content