Machine learning interatomic potential for silicon-nitride (Si₃N₄) by active learning

Silicon nitride (Si₃N₄) is a material with a high melting point (1900 °C), which finds applications in the automotive industry, where it is used in engine fabrications;¹ in the renewable energy industry, as an anti-reflective coating for solar cells;² in the engines of the NASA space shuttle, as armorial bearing;³ and in microelectronics, as a dielectric in nanoscale metal-oxide-semiconductor field-effect transistors (MOSFETs).⁴ In this last application, it is particularly valued for its use in the reduction of leakage currents and the prevention of boron diffusion. Moreover, Si₃N₄ is more stable than other silicon-nitrogen compounds^5,6 and it can crystallize in different phases, such as the trigonal α-Si₃N₄, the hexagonal β-Si₃N₄, and the cubic γ-Si₃N₄.

In order to study the properties of this material on atomistic scales, molecular dynamics (MD) simulations are frequently employed, in which the behavior of the atomic system is simulated by considering its interatomic interactions over time. Even though MD is a widely used technique in academia and industry, researchers employing it face a long-standing dilemma: accuracy vs efficiency.^7,8 On the one hand, ab initio methods, such as density functional theory (DFT), offer high accuracy; however, their prohibitively high computational costs limit their application to only relatively small systems (<10³ atoms) and short simulation times (∼10¹ ps). On the other hand, empirical potentials, which are hand-crafted parametrizations of the potential energy surface (PES), require much less computational resources; however, their accuracy is not sufficient for many applications. Examples of previously designed empirical potentials for Si₃N₄ include Tersoff,⁹ Billeter,¹⁰ Born–Mayer–Huggins (BMH),¹¹ and ReaxFF.¹² These empirical potentials are orders of magnitude faster than DFT and accurate enough for the specific applications for which they were developed, such as studying diffusion,⁹ the interactions with hydrogen atoms,⁹ and interfaces with other materials, such as Si.¹⁰ However, their accuracy is insufficient for several other important applications, including modeling amorphous Si₃N₄, which is relevant in the semiconductor industry since the majority of thin films deposited for use in transistors are amorphous. The modeling of amorphous Si₃N₄ is typically done using the melt-and-quench technique, which is the direct motivation for the development of this universally applicable machine learning (ML)-based potential. Melt-and-quench is a technique used to model amorphous materials by running extensive MD simulations. Hereby, the atomic system is heated above its melting point, equilibrated in the liquid phase, and subsequently quenched to room temperature in an amorphous solid.

In order to access realistic amorphous structures, ab initio accuracy combined with comparatively low quenching rates of 5 K/ps or below is required.^13,14 This translates into higher numbers of time-steps (∼10⁶), which is out of range for DFT due to its high computational costs. On the other hand, empirical potentials are a very useful tool to run these extensive MD simulations due to their low computational costs; however, their accuracy is not always near that of ab initio methods.

In this context, ML interatomic potentials have gained a lot of interest in the recent past due to their attractive cost-to-accuracy ratio, as they are very often able to combine ab initio level accuracy with computational times similar to those of empirical potentials. These interatomic potentials work by training a ML model, such as kernel-based techniques^15–17 or neural networks,^18–20 to act as a computationally more efficient surrogate model of the PES. Moreover, their popularization has been further aided by readily available software packages.^21–24 

In this work, we present a general-purpose ML interatomic potential for Si₃N₄. Our approach combines an established ML model with a simple active-learning algorithm in a bootstrap process. This allowed us to semi-automatically develop an accurate ML interatomic potential, starting from an inexpensive empirical potential, poorly optimized for amorphous Si₃N₄.

Our proposed ML interatomic potential is validated by comparing the structural properties of our resulting amorphous Si₃N₄ models to those of DFT calculations, empirical potentials, and experimental data found in the literature. To the best of the authors’ knowledge, this is the first time a ML interatomic potential has been developed for Si₃N₄.

This paper is divided into the following parts: Sec. II presents the ML model, the dataset creation, and the proposed active-learning algorithm; Sec. III presents and discusses the results of validating our proposed ML interatomic potential for MD simulations against empirical potentials and experimental data and for atomic calculations against DFT, where the accuracy and computational costs are considered; finally, Sec. IV presents a summary of this work.

We implemented our ML interatomic potential within the well-established Gaussian approximation potential (GAP)²⁵ method, represented in Fig. 1. This is a kernel-based method, successfully used in previous developments of ML interatomic potentials for different materials.^15–17 Similar to other ML models employed in interatomic potentials, GAP computes the total potential energy of a given atomic system as the sum of a local energy contribution from each atom. GAP assumes that similar local atomic environments make similar contributions to the total potential energy. Therefore, it computes the local energy contribution of a given atom by comparing its local atomic environment to the local atomic environments of the atomic configurations available in its training dataset. This comparison is made by a kernel function, which gives name to this kind of method.

Similar to other ML interatomic potentials, GAP requires a descriptor to act as an intermediary between the atomic structures and the ML model itself. A descriptor is a mathematical function that maps the three-dimensional, or Euclidean, coordinates of an atomic system to a ML-compatible representation, in the form of a matrix or vector. This new representation must be invariant with respect to translations and rotations of the whole system, as well as permutations of identical atoms. Descriptors can be either global (describing the entire atomic system) or local (describing the local atomic environment of a given atom in the system), with the latter being more commonly used in the development of ML interatomic potentials for bulk systems.

In this work, we paired the GAP with a set of three descriptors: smooth overlap of atomic positions (SOAP),^26, two-body, and three-body.¹⁷ Our implementation of GAP and its descriptors was done using the software package QUIP;²³ its most relevant parameters are summarized in Table I. Broadly speaking, δ scales the kernel per descriptor, thereby determining the importance of each descriptor in the GAP computations; r_cut is the descriptor cutoff, which determines the distance at which neighbor atoms are considered to describe the local atomic environments; r_Δ is the cutoff transition width, which determines the distance needed for the descriptor cutoff to be smoothly taken to zero; n_max and l_max are the number of angular and radial basis functions for the SOAP descriptor, respectively; and ζ is the power the kernel is raised to. For the complete list of parameters and further details on the implementation, the reader is referred to our open-access repository, where the ML interatomic potential and the training dataset are publicly available.²⁷ 

Our proposed active-learning algorithm is based on an iterative re-training approach, derived from Deringer and Csányi’s original algorithm.¹⁷ Our algorithm allows us to start with a simple and inexpensive empirical potential for an initial sample of the PES of interest and eventually build up a ML interatomic potential with near ab initio accuracy. This approach aims at saving computational time since it avoids running an initial MD simulation with DFT. Moreover, due to its iterative nature, this active-learning algorithm automatically runs and improves the ML interatomic potential until the desired accuracy is reached.

For the generation of the first training dataset, we run an MD simulation with the inexpensive empirical Billeter potential,¹⁰ in QuantumATK.²⁸ During this initial run, an N–N repulsive term was included to suppress the formation of N–N bonds, with the aim of reducing the number of under-coordinated atoms in the resulting amorphous structures.

This initial MD started from a defect-free α-Si₃N₄ crystal containing 224 atoms at 300 K, which was thereafter heated above its melting point to 5000 K and subsequently quenched back to 300 K in 10⁶ steps, using a step size of 0.5 fs, resulting in a quench rate of 4.7 K/fs. The resulting trajectory was subsequently sub-sampled, leading to 565 atomic configurations. Furthermore, single-atom configurations for Si and N (needed for the GAP to learn the correct free energies) were added to the training dataset, together with 267 dimers (Si–Si, Si–N, and N–N).

The energies and forces for these atomic configurations were re-calculated with DFT in the software package CP2K,²⁹ using the Perdew–Burke–Ernzerhof (PBE) functional³⁰ in conjunction with the Goedecker–Teter–Hutter (GTH) pseudopotentials,³¹ as these were found to be suitable for the system of interest.

Once the initial training dataset was ready, a GAP was trained on it. This GAP was used to run a series of sequential melt-and-quench MD simulations. After each MD, the last snapshots from the trajectory were saved and their energies, forces, and/or stress tensors were re-calculated using DFT and added to the training dataset. The GAP was re-trained on the extended training dataset before running the next MD, and the process was repeated until the desired accuracy was reached. We applied this technique to α-Si₃N₄ and β-Si₃N₄ with different densities, as well as to bulk Si, in order to provide a thorough sampling of the PES of interest. Table II shows the structure types added in each iteration.

The selection of the atomic configurations from which to start running the MD is manually done by the user, based on the application of interest. Thereafter, the sampling of the PES is automatically done by the algorithm, thereby constituting a semi-automatic active learning method.

The results presented in the right panel of Fig. 2 show how GAPs from earlier iterations of the training process were prone to creating non-homogeneous atomic structures when used to model amorphous Si₃N₄ and GAPs from intermediate iterations still yielded chains and clusters of Si and N. These unphysical features were not present in the quenched structures after the final iteration of the re-training cycle.

Table III shows the energy, force, and stress σ parameters, which determine the convergence criteria in the GAP training process used for each structure type. We note that the DFT stress tensors were only calculated for the last five iterations.

The training process of the GAP consists in an optimization problem: finding the set of α values that produces the lowest errors for the atomic configurations in the training dataset. If not specified, the convergence parameters will be the same for all atomic configurations, in which case the GAP would overfit to the atomic configurations with the highest energies and/or forces. Therefore, the data are divided into four blocks, each with distinct σ parameters, as a measure to avoid overfitting to configurations with higher energies and/or forces.

The resulting GAP was validated against energies and forces from DFT for the atomic structures in a testing dataset; the results are shown in Fig. 3. An equivalent plot for the empirical potentials is shown in the supplementary material. The testing dataset contains 1000 liquid and amorphous atomic configurations, randomly selected from the initial MD, none of which was used in the training process. The deviations were estimated by computing the mean absolute error (MAE) between both calculation methods. The energy and force calculations yielded a MAE of 8 meV/atom and 0.26 eV/Å, respectively. Moreover, the clear linear correlation between GAP and DFT indicates no systematic error.

Thereafter, the GAP was used to run a melt-and-quench MD, starting from a crystalline defect-free Si₃N₄ system composed of 224 atoms. The initial velocities were drawn from a Maxwell–Boltzmann distribution at 300 K. Then, the system was heated above the melting temperature, at 5000 K, thereby becoming a liquid. An equilibration phase continued at 5000 K for 10.000 steps, followed by a quenching phase over 10⁶ steps, which drove the system back to a temperature of 300 K. A time-step of 0.5 fs was used, resulting in a quenching rate of 4.7 K/fs.

A DFT reference atomic structure was obtained by running a geometry optimization of the final snapshot of the melt-and-quench MD performed with the Billeter empirical potential.¹⁰ The geometry optimization was done with DFT using the PBE³⁰ functional in CP2K.²⁹ This technique to build a reference DFT structure was preferred over running a DFT MD due to the high computational costs of the latter. As will be shown in this section, this reference structure correlates well with experimental data.

Both the GAP and DFT atomic structures were first compared by computing their total and partial radial density functions (RDFs); the results are presented in Fig. 4. The top panel compares the total RDF (i.e., all atomic species were considered) for GAP, DFT, and the empirical potentials Tersoff,⁹ ReaxFF,¹² and Billeter.¹⁰ The amorphous Si₃N₄ for each of these empirical potentials was obtained by running a melt-and-quench MD in QuantumATK,²⁸ with a Nose–Hoover NVT thermostat, starting from the exact same initial Si₃N₄ atomic configuration used for the GAP MD. The same recipe was used as in the GAP MD, with identical time-step and quench-rate. Moreover, an amorphous Si₃N₄ structure generated with the BMH empirical potential, obtained from Ref. 11, was added to the comparison. This BMH structure contains 112 atoms, has a density of 2.9 g/cm³, and was generated using a quench rate of 0.1 K/fs. In addition, we validated our results by comparing them with experimental data on amorphous Si₃N₄.³² A comparison between the RDF of the BMH structure and a GAP structure produced by a melt-and-quench MD under the same conditions can be found in the supplementary material.

The amorphous Si₃N₄ structures created with the ReaxFF and Billeter empirical potentials yield high peaks for short interatomic distances, which do not correlate with experimental data. The Tersoff and BMH empirical potentials yield good agreement with the DFT reference, as well as the experimental data, over the entire range of relevant interatomic distances. The GAP correlation to these references is slightly better than that of the empirical potentials, particularly if the height of the first and secondary peaks is considered.

Figure 5 compares the distributions of Si–N bond-lengths and N–Si–N bond-angles of these atomic structures. The distributions in the structure generated with the BMH empirical potential are added for reference purposes. The distributions in the structures generated with Tersoff, Billeter, and ReaxFF can be found in the supplementary material. They are not included here, as the high errors found in them make the comparison unnecessary. The bond-length distribution is virtually identical for GAP and DFT; both share minimum and maximum values at 1.6 and 2.4 A, respectively, as well as a single peak at 1.79 A for DFT and 1.78 A for GAP. On the other hand, for the bond-angle distributions, the minimum and maximum values for DFT/GAP are 72.0°/75.2° and 176.3°/176.8°, respectively. The mean angle for GAP is 109.0°, while for DFT, it is 109.2°. Overall, the GAP yields a better agreement with the DFT reference than the BMH empirical potential, which slightly overestimates the angles and underestimates the bond-lengths.

A further validation test was performed by computing the neutron scattering factor for the GAP and DFT amorphous Si₃N₄ structures. The results are shown in Fig. 6, where experimental data for an amorphous Si₃N₄ sample with a density of 2.60 g cm⁻³ (Ref. 32) were added for reference purposes. We note that the density of our GAP amorphous Si₃N₄ was 2.92 g cm⁻³. The similarity between the GAP results and the experimental data is remarkable for low values of q. For higher values, the discrepancy is more noticeable; however, the GAP and DFT results remain in good agreement, thereby hinting at limitations at the theory level and not in the training of the ML interatomic potential. If more accurate results were needed for this, or any other specific application, the same framework could be used with a different theory level.¹⁴ 

The main application of our ML interatomic potential is to be employed in running melt-and-quench MD to model amorphous Si₃N₄ with vanishing defect densities, starting from different initial configurations, which vary in crystal type and density. Therefore, to analyze the defects in the structure generated with our potential, its coordination number was computed and compared to that of the DFT reference structure. Figure 7 shows the coordination numbers of Si atoms (top) and N atoms (bottom). The qualitative similarity between both coordination number distributions is clear.

In order to further test our proposed ML interatomic potential, the vibrational density of states (VDOS) of the amorphous Si₃N₄ structure generated using the melt-and-quench technique was computed.

The finite difference calculations were carried out by displacing a single atom by 0.01 Å along each coordinate axis. The resulting dynamical matrix was then diagonalized in mass-weighted coordinates to obtain the system’s normal modes and the corresponding phonon energies.

The results of these calculations are presented in Fig. 8. As can be seen, there is excellent agreement between the VDOS spectra obtained by DFT calculations and our ML interatomic potential over the whole relevant energy range.

The good match is partially due to including forces in the training process, as this provides the ML interatomic potential with the explicit value of the first derivative of the PES with respect to the atomic coordinates and a much more accurate representation of the PES curvature, reducing the numerical errors when computing the VDOS.

As mentioned earlier, the main advantage of ML interatomic potentials over ab initio methods when running MD is their more attractive accuracy-to-cost ratio. Maintaining a level of accuracy comparable to that of DFT, while being orders of magnitude computationally cheaper, allows for the use of ML interatomic potentials in MD simulations with system sizes and simulations times previously out of reach.

This section will analyze the computational costs of the proposed ML interatomic potential, for which we ran two different tests. In the first test, we analyzed how the time needed to run an MD simulation scales with the number of atoms in the system. Therefore, we considered an Si₃N₄ structure and we ran an MD simulation with 112, 378, 896, 1750, and 3024 atoms. The MD simulation was performed in Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) using 20 cores, and it consisted of 10.000 steps at 1000 K, with a time-step of 1 fs. The results can be seen in the top panel of Fig. 9, where a nearly linear scaling of the simulation time and the number of atoms can be observed.

The second test aimed at analyzing how parallelization can be used to reduce the computational times needed to run MD simulations with our proposed ML interatomic potential. Therefore, we used LAMMPS to re-run an MD simulation, which consisted in an Si₃N₄ system with 3024 atoms, evolving over 10.000 times at 1000 K, with a time-step of 1 fs. The simulation was repeated with 2, 4, 8, 16, and 32 cores, and the results are presented in the bottom panel of Fig. 9, where a nearly linear correlation between performance and the number of cores can be seen.

By running an MD simulation with DFT for the same system using 4 cores, we estimated its performance to be about 0.0011 (K steps/h). Our GAP yields a performance of 1.28 (K steps/h) under the same circumstances, thereby making it about 3–4 orders of magnitude faster than DFT for this particular scenario. However, the computational cost of DFT scales cubically with the number of atoms in the system, as shown in the supplementary material. Therefore, the performance gap between DFT and our ML interatomic potential grows with system size.

We have presented a machine learning (ML) interatomic potential for silicon nitride (Si₃N₄) using the Gaussian approximation potential (GAP) method. The training was done using an active-learning algorithm, consisting of iterative re-trainings. Following our approach, a GAP was used to run molecular dynamics (MD). The resulting atomic trajectory was sub-sampled, and the energies, forces, and/or stress tensors of the final atomic configurations were re-calculated using density functional theory (DFT). These data were added to the training dataset, the GAP was re-trained, and the process was repeated until the desired accuracy was reached. By applying this method to a variety of different atomic systems, we were able to thoroughly sample the potential energy surface (PES) of interest. Compared to energies from DFT calculations, our ML potential yields a mean absolute error (MAE) of ∼8 eV/atom. Modeling amorphous Si₃N₄ with our GAP yields remarkable agreement with DFT and experimental data. Moreover, our GAP was approximately 3–4 orders of magnitude faster than DFT, thereby making it an attractive enabler for much larger and more realistic structures, allowing for MD simulations over much larger atomic system sizes and time scales.

See the supplementary material for the results of further tests comparing our proposed ML interatomic potential to experimental data, DFT, and empirical potentials.

This project received funding from the European Union’s Horizon 2020 research and innovation program under grant agreement No. 871813, within the framework of the project Modeling Unconventional Nanoscaled Device FABrication (MUNDFAB). The financial support provided by the Austrian Federal Ministry for Digital and Economic Affairs; the National Foundation for Research, Technology and Development; and the Christian Doppler Research Association is gratefully acknowledged. The authors acknowledge the support from the Vienna Scientific Cluster (VSC). G.S. acknowledges the support of the Exascale Computing Project (17-SC-20-SC), a collaborative effort of the U.S. Department of Energy Office of Science and the National Nuclear Security Administration.

The authors have no conflicts to disclose.

Diego Milardovich: Conceptualization (lead); Data curation (equal); Investigation (equal); Methodology (equal); Software (lead); Validation (equal); Visualization (lead); Writing – original draft (lead). Christoph Wilhelmer: Data curation (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Writing – review & editing (equal). Dominic Waldhoer: Data curation (equal); Investigation (equal); Methodology (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal). Lukas Cvitkovich: Resources (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). Ganesh Sivaraman: Formal analysis (equal); Resources (equal); Validation (equal); Writing – review & editing (equal). Tibor Grasser: Funding acquisition (lead); Resources (equal); Supervision (equal); Writing – review & editing (equal).

The GAP potential is stored in a GitHub repository, in an XML format, together with the training dataset of atomic structures used in this paper.²⁷