BKS Potential Research Articles

Mechanical behavior of alpha quartz with void defects under tension: a molecular dynamics study using different interatomic potentials

This study employs a series of molecular dynamics (MD) simulations, utilizing three commonly used interatomic potentials, i.e. van Beest, Kramer, and van Santen (BKS), Vashishta, and Tersoff to analyze the structural and mechanical characteristics within both void-free and single-void α-quartz configurations. Two distinct ensembles, NVT and NPT, were separately applied to investigate the tensile response. The validation of MD results included a comparative study of the three potentials as well as a comparison with experimental microstructural and tension studies. While BKS and Vashishta potentials accurately calculated the bond lengths, density and lattice parameters compared to the experimental values for void-free α-quartz, the results obtained with Tersoff potential exhibited relatively large deviations. The BKS potential offered an accurate description of the mechanical response of α-quartz by successfully predicting stress–strain curves. The Vashishta potential overpredicted Young’s modulus as compared to BKS. The Tersoff potential could capture the elastic deformation but was unable to predict the fracture behavior. The presence of a spherical void significantly reduced mechanical behavior of α-quartz, and the extent of this reduction was highly related to void size. When applying the BKS potential with an NVT ensemble, the ultimate tensile strengths decreased by 19% and 72% with void sizes of 2.5 and 15 Å, respectively. Equivalent stress analysis reveals that the BKS potential can effectively capture greater stress concentration around the void compared to other two potentials. Based on the comparison study, the BKS potential seems to be the most suitable one to describe α-quartz under tension in a realistic manner.

Mechanical and microstructural response of densified silica glass under uniaxial compression: Atomistic simulations**Project supported by the National Natural Science Foundation of China (Grant Nos. 51727807 and 11875318), Beijing Institute of Technology Research Fund Program for Young Scholars, and Yue Qi Young Scholar Project in CUMTB.

We investigate the mechanical and microstructural changes of the densified silica glass under uniaxial loading-unloading via atomistic simulations with a modified BKS potential. The stress–strain relationship is found to include three respective stages: elastic, plastic and hardening regions. The bulk modulus increases with the initial densification and will undergo a rapid increase after complete densification. The yield pressure varies from 5 to 12 GPa for different densified samples. In addition, the Si–O–Si bond angle reduces during elastic deformation under compression, and 5-fold Si will increase linearly in the plastic deformation. In the hardening region, the peak splitting and the new peak are both found on the Si–Si and O–O pair radial distribution functions, where the 6-fold Si is increased. Instead, the lateral displacement of the atoms always varies linearly with strain, without evident periodic characteristic. As is expected, the samples are permanently densified after release from the plastic region, and the maximum density of recovered samples is about 2.64 g/cm3, which contains 15 % 5-fold Si, and the Si–O–Si bond angle is less than the ordinary silica glass. All these findings are of great significance for understanding the deformation process of densified silica glass.

A novel method to predict the thermal conductivity of nanoporous materials from atomistic simulations

The low density (smaller than 10% of the bulk density) and the nanostructured porosity of silica aerogels provide their extremely low thermal conductivities but also impact their poor mechanical properties. Atomic scale simulation is the appropriate tool to predict the thermal and mechanical properties of such materials. For such simulations, the interatomic potential should be carefully chosen to ensure result validity but also reasonable computational times. A truncated BKS potential has been used for aerogels as it fairly reproduces the nanostructure. It allows reducing the computational time by a 3000 gain factor on the CPU time per atom per step compared to the original BKS interatomic potential while predicting correctly the mechanical properties. However, when it comes to skeletal thermal conductivity of nanoporous silica, the associated computation times are too large for a representative volume. This is due to the low thermal diffusivity of the material. Here, a new method that takes advantage of the amorphous structure of silica and the diffusive nature of phonon heat transfer at the scale of an aerogel aggregate is proposed. The time dependent temperature profile in the system obtained from Non-Equilibrium Molecular Dynamics simulations is compared to the classical solution of the thermal diffusion equation and an identification procedure is used to determine the thermal conductivity of silica aerogels.

Mechanical properties of silica glass predicted by a pair-wise potential in molecular dynamics simulations

Pair-wise BKS potential for silica glass was modified by changing the short-range cutoff to improve the prediction of mechanical properties. Classical molecular dynamic (MD) simulations were carried out by using this modified potential and better prediction of mechanical properties at ambient conditions and their changes with pressure at room temperature were shown. It was further observed that there was no major change in structural properties when different short-range cutoffs were used to relax the glass. To better understand the contribution of structure to the mechanical properties, samples quenched under various cooling rates and one structure optimized in first principles calculations were used. It was seen that the improvement of structure has a less effect than improving the interaction potential in determining the mechanical properties of BKS silica glass. The strength and weakness of the BKS potential for silica glass was discussed.

Molecular dynamics simulations of amorphous silica surface properties with truncated Coulomb interactions

Molecular Dynamics simulations have been performed to evaluate the capacity of the van Beest, Kramer and van Santen (BKS) [1] potential with truncated Coulomb interactions as proposed by A. Carré et al. [2] to reproduce amorphous silica surface properties. We compare the results obtained with the truncated BKS potential with those obtained from its full-range interaction version. Energies of (SiO2)n clusters are computed for both potentials. The energies and structural properties of small aggregates of silica (from ≈3.5nm to 7.6nm of diameter) are investigated at high temperature. Both potentials lead to the same results for those properties when considering the transition from the core to the shell of aggregates. The significant computation time saved with this cut-off on Coulomb interactions allows for large-scale simulations of silica aerogels.

On force fields for molecular dynamics simulations of crystalline silica

This paper reviews and examines interatomic potentials or force fields for molecular dynamics (MD) simulation of crystalline silica. The investigated potentials are the BKS, Pedone, Munetoh, TTAM, and CHIK. The calculated values of the lattice constants, density, radial and bond-angle distribution functions, equations of state, and phase transitions using different potentials are compared to experimental values for polymorphs of silica: quartz, cristobalite, coesite, and stishovite. Simulation results with the BKS potential accurately predict the experimental measurements to within 2%, and converge within a reasonable timeframe on an average workstation. The Pedone potential, also parameterized for other metallic oxides, computationally is slightly more expensive and is not as accurate. The simulations with both the CHIK and TTAM potentials are less accurate than with the BKS potential for modeling silica over the entire range of the phase diagram. The simulations with the Munetoh potential are by far the cheapest in terms of the modest computational requirements, but unsuitable for modeling crystalline silica. It could not produce the nature of the α–β and I–II phase transitions in quartz or the equation of state for stishovite silica, and the predicted structural properties sometimes differ from experimental values by more than 10%.

Enhanced thermal characterization of silica aerogels through molecular dynamics simulation

Porous structures of silica aerogels are generated using classical molecular dynamics, with the Tersoff potential, which has been re-parametrized for modeling silicon dioxides. This work demonstrates that this potential is superior to the widely used BKS potential in terms of characterizing the thermal conductivities of amorphous silica, by comparing the vibrational density of states with previous experimental studies. Aerogel samples of increasing densities are obtained through an expanding, heating and quenching process. Reverse non-equilibrium molecular dynamics is applied at each density to determine the thermal conductivity. A power-law fit of the results is found to accurately reflect the power-law variation found in experimental bulk aerogels. The results are also of the same order of magnitude as experimental bulk aerogels, but they are consistently higher. By analyzing the pore size distribution on different simulation length scales, we show that such a disparity is due to finite sizes of pores that can be represented, where increasing simulation length scales lead to an increase in the largest pore size that can be modeled.

Water and other tetrahedral liquids: order, anomalies and solvation

In order to understand the common features of tetrahedral liquids with water-like anomalies, the relationship between local order and anomalies has been studied using molecular dynamics simulations for three categories of such liquids: (a) atomistic rigid-body models for water (TIP4P, TIP4P/2005, mTIP3P, SPC/E), (b) ionic melts, BeF2 (TRIM model) and SiO2 (BKS potential) and (c) Stillinger–Weber liquids parametrized to model water (mW) and silicon. Rigid-body, atomistic models for water and the Stillinger–Weber liquids show a strong correlation between tetrahedral and pair correlation order and the temperature for the onset of the density anomaly is close to the melting temperature. In contrast, the ionic melts show weaker and more variable degrees of correlation between tetrahedral and pair correlation metrics, and the onset temperature for the density anomaly is more than twice the melting temperature. In the case of water, the relationship between water-like anomalies and solvation is studied by examining the hydration of spherical solutes (Na+, Cl−, Ar) in water models with different temperature regimes of anomalies (SPC/E, TIP4P and mTIP3P). For both ionic and nonpolar solutes, the local structure and energy of water molecules is essentially the same as in bulk water beyond the second-neighbour shell. The local order and binding energy of water molecules are not perturbed by the presence of a hydrophobic solute. In the case of ionic solutes, the perturbation is largely localized within the first hydration shell. The binding energies for the ions are strongly dependent on the water models and clearly indicate that the geometry of the partial charge distributions, and the associated multipole moments, play an important role. However the anomalous behaviour of the water network has been found to be unimportant for polar solvation.

Atomistic molecular dynamics simulations of shock compressed quartz

Atomistic non-equilibrium molecular dynamics simulations of shock wave compression of quartz have been performed using the so-called BKS semi-empirical potential of van Beest, Kramer, and van Santen [Phys. Rev. B 43, 5068 (1991)] to construct the Hugoniot of quartz. Our scheme mimics the real world experimental set up by using a flyer-plate impactor to initiate the shock wave and is the first shock wave simulation that uses a geometry optimised system of a polar slab in a three-dimensional system employing periodic boundary conditions. Our scheme also includes the relaxation of the surface dipole in the polar quartz slab which is an essential pre-requisite to a stable simulation. The original BKS potential is unsuited to shock wave calculations and so we propose a simple modification. With this modification, we find that our calculated Hugoniot is in good agreement with experimental shock wave data up to 25 GPa, but significantly diverges beyond this point. We conclude that our modified BKS potential is suitable for quartz under representative pressure conditions of the Earth core, but unsuitable for high-pressure shock wave simulations. We also find that the BKS potential incorrectly prefers the β-quartz phase over the α-quartz phase at zero-temperature, and that there is a β → α phase-transition at 6 GPa.

Simulation study of polymorphism and diffusion anomaly for SiO2 and GeO2 liquid

We report molecular dynamic (MD) simulations of a silica and germania liquid. We show for the first time that the liquid models constructed of the BKS and Ocffner-Elliott potential are composed of D4, D5 and D6 domains. These domains of different liquids (identical chemical composition but different density) are structurally identical. The liquid densification occurs gradually with variation in the volume of the D4, D5 and D6 domains. The void radius distributions (VRD) have been computed and they are well fitted to a Gaussian form. Both the average void radii and width of the VRD decrease with pressure. The diffusivity of Si/Ge or O depends linearly on the fraction of TOx (x = 4, 5 and 6; T is Ge or Si). The diffusion anomaly results from the very high mobility of silicon in the D5 domain compared to domains D4 and D6.

Unified interatomic potential for zircon, zirconia and silica systems

We report the development of a unified force field with high transferability and reliability to model both structures and mechanical properties of ZrSiO4, ZrO2 and SiO2 based on the success of the BKS potential for SiO2 (van Beest et al., Phys. Rev. Lett., 1990, 64, 1955). The thermal expansion, relative stability and phase transition properties are consistent with experimental data and DFT calculations. The zircon to reidite transition pressure is 6.7 GPa. Amorphization of zircon results in volume expansion of 11% and decrease in the bulk modulus of 60%. Si polymerization and Zr under-coordination were observed in amorphous ZrSiO4 and in the damage produced by a 10 keV Zr recoil in crystalline ZrSiO4. The upper limit of stored energy in amorphous ZrSiO4 is 140 kJ/mol.

Molecular dynamics simulation of ultrafast laser ablation of fused silica film

Ultrafast laser ablation of fused silica is studied using molecular dynamics simulations. Ionization and generation of free electrons, absorption of the laser energy by free electrons and energy coupling between free electrons and ions are considered. The BKS potential is applied and modified to describe molecular interactions and the effect of free electrons. Smooth particle mesh of the Ewald method (SPME) is adopted to calculate the Coulomb force. It is found that the electrostatic Coulomb force, which is caused by the ionization, plays an important role in the laser ablation process.

Interatomic Potentials of Silica Glass for Molecular Dynamics Simulation at High Temperature

Multi-body interaction and atomic charge transfer played essential roles in depicting the ionic-covalent bonds of Si—O systems. The limitations of 3-body terms in interatomic potentials for silica glass were analyzed according to molecular dynamics studies at high temperature, and it was found that pair-wise ionic potentials were better in simulating the microstructure evolution of silica glass than multi-body potentials. Moreover, the effects of atomic charge transfer on atomic self-diffusion were investigated, which showed that the activation temperature of atomic self-diffusion calculated with the Morse potential was lower than that calculated with the BKS potential due to the less atomic charge transfer in Morse potential. It was also found that the atomic self-diffusion coefficients calculated increased with the decrease in atomic charge transfer at the same temperature. Therefore, the pair-wise potentials with less atomic charge transfer should be in favor of investigating the dynamic properties of silica glass at high temperature.

Molecular dynamics simulation of ultrafast laser ablation of fused silica

Ultrafast laser ablation of fused silica is studied using molecular dynamics simulations. Ionization and generation of free electrons, absorption of the laser energy by free electrons and energy coupling between free electrons and ions are considered. The BKS potential is used to describe molecular interactions and is modified to include the effect of free electrons. Temperature, density, and material removal are computed, and the thermal and non-thermal mechanisms of ultrafast laser ablation of fused silica are discussed.

Diffusional anomaly and network dynamics in liquid silica

The present study applies the power spectral analysis technique to understand the diffusional anomaly in liquid silica, modeled using the Beest-Kramer-van Santen (BKS) potential. Molecular-dynamics simulations have been carried out to show that power spectrum of tagged particle potential energy of silica shows a regime with 1f(alpha) dependence on frequency f which is the characteristic signature of multiple time scale behaviour in networks. As demonstrated earlier in the case of water [J. Chem. Phys. 122, 104507 (2005)], the variations in the mobility associated with the diffusional anomaly are mirrored in the scaling exponent alpha associated with this multiple time scale behavior. Our results indicate that in the anomalous regime, as the local tetrahedral order decreases with temperature or pressure, the coupling of local modes to network reorganizations increases and so does the diffusivity. This symmetry-dependence of the vibrational couplings is responsible for the connection between the structural and diffusional anomalies.

Test of classical nucleation theory on deeply supercooled high-pressure simulated silica

We test classical nucleation theory (CNT) in the case of simulations of deeply supercooled, high density liquid silica, as modeled by the van Beest-Kramer-van Santen potential. We find that at density rho=4.38 gcm(3), spontaneous nucleation of crystalline stishovite occurs in conventional molecular dynamics simulations at temperature T=3000 K, and we evaluate the nucleation rate J directly at this T via "brute force" sampling of nucleation events in numerous independent runs. We then use parallel, constrained Monte Carlo simulations to evaluate DeltaG(n), the free energy to form a crystalline embryo containing n silicon atoms, at T=3000, 3100, 3200, and 3300 K. By comparing the form of DeltaG(n) to CNT, we test the ability of CNT to reproduce the observed behavior as we approach the regime where spontaneous nucleation occurs on simulation time scales. We find that the prediction of CNT for the n dependence of DeltaG(n) fits reasonably well to the data at all T studied. Deltamu, the chemical potential difference between bulk liquid and stishovite, is evaluated as a fit parameter in our analysis of the form of DeltaG(n). Compared to directly determined values of Deltamu extracted from previous work, the fitted values agree only at T=3300 K; at lower T the fitted values increasingly overestimate Deltamu as T decreases. We find that n(*), the size of the critical nucleus, is approximately ten silicon atoms at T=3300 K. At 3000 K, n(*) decreases to approximately 3, and at such small sizes methodological challenges arise in the evaluation of DeltaG(n) when using standard techniques; indeed even the thermodynamic stability of the supercooled liquid comes into question under these conditions. We therefore present a modified approach that permits an estimation of DeltaG(n) at 3000 K. Finally, we directly evaluate at T=3000 K the kinetic prefactors in the CNT expression for J, and find physically reasonable values; e.g., the diffusion length that Si atoms must travel in order to move from the liquid to the crystal embryo is approximately 0.2 nm. We are thereby able to compare the results for J at 3000 K obtained both directly and based on CNT, and find that they agree within an order of magnitude. In sum, our work quantifies how certain predictions of CNT (e.g., for Deltamu) break down in this deeply supercooled limit, while others [the n dependence of DeltaG(n)] are not as adversely affected.

Comparison of model potentials for molecular-dynamics simulations of silica

Structural, thermomechanical, and dynamic properties of pure silica SiO2 are calculated with three different model potentials, namely, the potential suggested by van Beest, Kramer, and van Santen (BKS) [Phys. Rev. Lett. 64, 1955 (1990)], the fluctuating-charge potential with a Morse stretch term for the short-range interactions proposed by Demiralp, Cagin, and Goddard (DCG)[Phys. Rev. Lett. 82, 1708 (1999)], and a polarizable force field proposed by Tangney and Scandolo (TS) [J. Chem. Phys. 117, 8898 (2002)]. The DCG potential had to be modified due to flaws in the original treatment. While BKS reproduces many thermomechanical properties of different polymorphs rather accurately, it also shows qualitatively wrong trends concerning the phononic density of states, an absence of the experimentally observed anomaly in the c/a ratio at the quartz alpha-beta transition, pathological instabilities in the beta-cristobalite phase, and a vastly overestimated transition pressure for the stishovite I --> II transition. These shortcomings are only partially remedied by the modified DCG potential but greatly improved by the TS potential. DCG and TS both reproduce a pressure-induced transition from alpha-quartz to quartz II, predicted theoretically based on the BKS potential.

Large well-relaxed models of vitreous silica, coordination numbers, and entropy

A Monte Carlo method is presented for the simulation of vitreous silica. Well-relaxed networks of vitreous silica are generated containing up to 300 000 atoms. The resulting networks, quenched under the BKS potential, display smaller bond-angle variations and lower defect concentrations, as compared to networks generated with molecular dynamics. The total correlation functions $T(r)$ of our networks are in excellent agreement with neutron scattering data, provided that thermal effects and the maximum inverse wavelength used in the experiment are included in the comparison. A procedure commonly used in experiments to obtain coordination numbers from scattering data is to fit peaks in $\mathrm{rT}(r)$ with a gaussian. We show that this procedure can easily produce incorrect results. Finally, we estimate the configurational entropy of vitreous silica.

Classical and ab-initio molecular dynamic simulation of an amorphous silica surface

We present the results of a classical molecular dynamic simulation as well as of an ab initio molecular dynamic simulation of an amorphous silica surface. In the case of the classical simulation we use the potential proposed by van Beest et al. (BKS) whereas the ab initio simulation is done with a Car-Parrinello method (CPMD). We find that the surfaces generated by BKS have a higher concentration of defects (e.g. concentration of two-membered rings) than those generated with CPMD. In addition also the distribution functions of the angles and of the distances are different for the short rings. Hence we conclude that whereas the BKS potential is able to reproduce correctly the surface on the length scale beyond approx 5 Angstroems, it is necessary to use an ab initio method to predict reliably the structure at small scales.

Computer simulations of undercooled fluids and the glass transition

Two model studies are presented that attempt to describe the static and dynamic properties of glass-forming fluids via molecular dynamics simulations: The first model is an atomistically realistic model of SiO 2, the second model provides a coarse-grained description of polymer liquids, i.e., typical `fragile' glassformers, while SiO 2 is the prototype of a `strong glassformer'. For both models, attention is given to the questions as to which range of temperatures are properties in equilibrium, and whether such simulations can help to interpret experiments and/or check theoretical predictions. While in the simulation of SiO 2 using the potential of van Beest, Kramer and van Santen (`BKS potential') only temperatures T⩾2750 K can be equilibrated, nevertheless many properties can be extracted that are in almost quantitative agreement with experiment. Sodium silicate liquids investigated at temperatures T⩾1900 K are in agreement to within a few percent with experimental data, too. At very high temperatures (T⩾2750 K) , even for SiO 2 a temperature range is found where the mode coupling theory (MCT) of the glass transition is applicable (with T c =3330 K ). The polymer melt is described by flexible polymer chains represented by a continuum bead-spring model with a Lennard–Jones interaction between the effective monomers. The typical features of undercooled macromolecular fluids are reproduced qualitatively, with the properties at small length scales (≪1 nm) described by MCT, while at larger length scales a crossover to Rouse-like dynamics occurs.