Silicon nitride (Si3N4) is an extensively used material in the automotive, aerospace, and semiconductor industries. However, its widespread use is in contrast to the scarce availability of reliable interatomic potentials that can be employed to study various aspects of this material on an atomistic scale, particularly its amorphous phase. In this work, we developed a machine learning interatomic potential, using an efficient active learning technique, combined with the Gaussian approximation potential (GAP) method. Our strategy is based on using an inexpensive empirical potential to generate an initial dataset of atomic configurations, for which energies and forces were recalculated with density functional theory (DFT); thereafter, a GAP was trained on these data and an iterative re-training algorithm was used to improve it by learning on-the-fly. When compared to DFT, our potential yielded a mean absolute error of 8 meV/atom in energy calculations for a variety of liquid and amorphous structures and a speed-up of molecular dynamics simulations by 3–4 orders of magnitude, while achieving a first-rate agreement with experimental results. Our potential is publicly available in an open-access repository.
I. INTRODUCTION
Silicon nitride (Si3N4) is a material with a high melting point (1900 °C), which finds applications in the automotive industry, where it is used in engine fabrications;1 in the renewable energy industry, as an anti-reflective coating for solar cells;2 in the engines of the NASA space shuttle, as armorial bearing;3 and in microelectronics, as a dielectric in nanoscale metal-oxide-semiconductor field-effect transistors (MOSFETs).4 In this last application, it is particularly valued for its use in the reduction of leakage currents and the prevention of boron diffusion. Moreover, Si3N4 is more stable than other silicon-nitrogen compounds5,6 and it can crystallize in different phases, such as the trigonal α-Si3N4, the hexagonal β-Si3N4, and the cubic γ-Si3N4.
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 (<103 atoms) and short simulation times (∼101 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 Si3N4 include Tersoff,9 Billeter,10 Born–Mayer–Huggins (BMH),11 and ReaxFF.12 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,9 the interactions with hydrogen atoms,9 and interfaces with other materials, such as Si.10 However, their accuracy is insufficient for several other important applications, including modeling amorphous Si3N4, which is relevant in the semiconductor industry since the majority of thin films deposited for use in transistors are amorphous. The modeling of amorphous Si3N4 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 (∼106), 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 techniques15–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 Si3N4. 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 Si3N4.
Our proposed ML interatomic potential is validated by comparing the structural properties of our resulting amorphous Si3N4 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 Si3N4.
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.
II. METHODOLOGY
A. ML model
We implemented our ML interatomic potential within the well-established Gaussian approximation potential (GAP)25 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.
Schematic of the GAP method. This ML interatomic potential works by comparing the input atomic structure (red) to the atomic structures in its training dataset (blue) by means of a kernel function (K). This algorithm estimates the potential energy as a sum of all kernel difference estimations, weighted by the α values, learned from fitting to the training dataset.
Schematic of the GAP method. This ML interatomic potential works by comparing the input atomic structure (red) to the atomic structures in its training dataset (blue) by means of a kernel function (K). This algorithm estimates the potential energy as a sum of all kernel difference estimations, weighted by the α values, learned from fitting to the training dataset.
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.17 Our implementation of GAP and its descriptors was done using the software package QUIP;23 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; rcut 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; nmax and lmax 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.27
Parameters used to implement the GAP and its SOAP, three-body, and two-body descriptors.
Parameter . | SOAP . | Three-body . | Two-body . |
---|---|---|---|
δ (eV) | 0.4 | 1 | 4 |
rcut (Å) | 4 | 3 | 4 |
rΔ (Å) | 1 | ⋯ | ⋯ |
nmax | 8 | ⋯ | ⋯ |
lmax | 4 | ⋯ | ⋯ |
ζ | 4 | ⋯ | ⋯ |
Parameter . | SOAP . | Three-body . | Two-body . |
---|---|---|---|
δ (eV) | 0.4 | 1 | 4 |
rcut (Å) | 4 | 3 | 4 |
rΔ (Å) | 1 | ⋯ | ⋯ |
nmax | 8 | ⋯ | ⋯ |
lmax | 4 | ⋯ | ⋯ |
ζ | 4 | ⋯ | ⋯ |
B. Active-learning algorithm
Our proposed active-learning algorithm is based on an iterative re-training approach, derived from Deringer and Csányi’s original algorithm.17 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,10 in QuantumATK.28 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 α-Si3N4 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 106 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,29 using the Perdew–Burke–Ernzerhof (PBE) functional30 in conjunction with the Goedecker–Teter–Hutter (GTH) pseudopotentials,31 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 α-Si3N4 and β-Si3N4 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.
Atomic structures added to the training dataset in every iteration.
Iteration No. . | Structure type: number of configurations . |
---|---|
1 | Single atoms: 2; dimers: 267; initial MD: 565 |
2 | Non-homogeneous amorphous Si3N4: 39 |
3 | Amorphous Si3N4 with N-chains: 82 |
4 | Amorphous Si3N4 with N-clusters: 30 |
5 | Crystal α-Si3N4: 99 |
6 | Bulk Si: 87 |
7 | Amorphous Si3N4: 101 |
8 | Geometry optimized Si3N4: 18 |
9 | Amorphous Si3N4 with under-coordinated Si: 60 |
10 | Amorphous Si3N4 with over-coordinated Si: 49 |
11 | Vibration of hexagonal β-Si3N4: 82 |
12 | Amorphous Si3N4 with density 2.71 g cm−3: 141 |
Iteration No. . | Structure type: number of configurations . |
---|---|
1 | Single atoms: 2; dimers: 267; initial MD: 565 |
2 | Non-homogeneous amorphous Si3N4: 39 |
3 | Amorphous Si3N4 with N-chains: 82 |
4 | Amorphous Si3N4 with N-clusters: 30 |
5 | Crystal α-Si3N4: 99 |
6 | Bulk Si: 87 |
7 | Amorphous Si3N4: 101 |
8 | Geometry optimized Si3N4: 18 |
9 | Amorphous Si3N4 with under-coordinated Si: 60 |
10 | Amorphous Si3N4 with over-coordinated Si: 49 |
11 | Vibration of hexagonal β-Si3N4: 82 |
12 | Amorphous Si3N4 with density 2.71 g cm−3: 141 |
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 Si3N4 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.
Schematic of the iterative re-training algorithm used to develop the ML interatomic potential for Si3N4 (left) and results obtained over time for subsequent generations (right). The ML interatomic potential is used to run MD, thereby generating a trajectory of atomic configurations. The last snapshots of this trajectory are collected; their energies, forces, and/or stress tensors are recalculated with DFT and then added to the training dataset. The ML interatomic potential is thereafter re-trained and used to run a new MD, thereby repeating the process until the desired accuracy is reached. The right panel presents the improvement in the results of using the ML interatomic potential for successive re-training iterations. The earliest iterations yield unphysical heterogeneous distributions of Si and N; intermediate iterations show unphysical N-chains and N–N bonds; and the final potential overcomes all the aforementioned issues and generates realistic and defect-free structures of Si3N4.
Schematic of the iterative re-training algorithm used to develop the ML interatomic potential for Si3N4 (left) and results obtained over time for subsequent generations (right). The ML interatomic potential is used to run MD, thereby generating a trajectory of atomic configurations. The last snapshots of this trajectory are collected; their energies, forces, and/or stress tensors are recalculated with DFT and then added to the training dataset. The ML interatomic potential is thereafter re-trained and used to run a new MD, thereby repeating the process until the desired accuracy is reached. The right panel presents the improvement in the results of using the ML interatomic potential for successive re-training iterations. The earliest iterations yield unphysical heterogeneous distributions of Si and N; intermediate iterations show unphysical N-chains and N–N bonds; and the final potential overcomes all the aforementioned issues and generates realistic and defect-free structures of Si3N4.
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.
Energy, force, and stress convergence parameters for every structure type used in the training process.
Structure type . | σenergy . | σforce . | σstress . |
---|---|---|---|
Dimers | 0.01 | 0.1 | ⋯ |
Amorphous Si3N4 with N-clusters | 0.01 | 0.1 | ⋯ |
Non-homogeneous amorphous Si3N4 | 0.01 | 0.1 | ⋯ |
Crystalline Si3N4 | 0.002 | 0.02 | ⋯ |
Amorphous Si3N4 with N-chains | 0.002 | 0.02 | ⋯ |
Bulk Si | 0.002 | 0.02 | ⋯ |
Initial MD | 0.002 | 0.02 | ⋯ |
Hexagonal β-Si3N4 | 0.002 | 0.02 | 0.02 |
Amorphous Si3N4 (2.71 g cm−3) | 0.002 | 0.02 | 0.02 |
Amorphous Si3N4 with over-coordinated Si | 0.0005 | 0.005 | 0.005 |
Amorphous Si3N4 with under-coordinated Si | 0.0005 | 0.005 | 0.005 |
Final amorphous Si3N4 | 0.0005 | 0.005 | 0.005 |
Geometry optimized Si3N4 | 0.0005 | 0.005 | ⋯ |
Single atoms | 0.0001 | 0.001 | ⋯ |
Structure type . | σenergy . | σforce . | σstress . |
---|---|---|---|
Dimers | 0.01 | 0.1 | ⋯ |
Amorphous Si3N4 with N-clusters | 0.01 | 0.1 | ⋯ |
Non-homogeneous amorphous Si3N4 | 0.01 | 0.1 | ⋯ |
Crystalline Si3N4 | 0.002 | 0.02 | ⋯ |
Amorphous Si3N4 with N-chains | 0.002 | 0.02 | ⋯ |
Bulk Si | 0.002 | 0.02 | ⋯ |
Initial MD | 0.002 | 0.02 | ⋯ |
Hexagonal β-Si3N4 | 0.002 | 0.02 | 0.02 |
Amorphous Si3N4 (2.71 g cm−3) | 0.002 | 0.02 | 0.02 |
Amorphous Si3N4 with over-coordinated Si | 0.0005 | 0.005 | 0.005 |
Amorphous Si3N4 with under-coordinated Si | 0.0005 | 0.005 | 0.005 |
Final amorphous Si3N4 | 0.0005 | 0.005 | 0.005 |
Geometry optimized Si3N4 | 0.0005 | 0.005 | ⋯ |
Single atoms | 0.0001 | 0.001 | ⋯ |
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.
III. RESULTS AND DISCUSSION
A. Atomic calculations
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.
Results of testing the ML interatomic potential against DFT for forces (top) and energies (bottom), for the atomic configurations in the testing dataset, together with their respective MAE. Atomic configurations were divided into liquid (blue) and amorphous (green).
Results of testing the ML interatomic potential against DFT for forces (top) and energies (bottom), for the atomic configurations in the testing dataset, together with their respective MAE. Atomic configurations were divided into liquid (blue) and amorphous (green).
B. Structural properties
Thereafter, the GAP was used to run a melt-and-quench MD, starting from a crystalline defect-free Si3N4 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 106 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.10 The geometry optimization was done with DFT using the PBE30 functional in CP2K.29 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,9 ReaxFF,12 and Billeter.10 The amorphous Si3N4 for each of these empirical potentials was obtained by running a melt-and-quench MD in QuantumATK,28 with a Nose–Hoover NVT thermostat, starting from the exact same initial Si3N4 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 Si3N4 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/cm3, 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 Si3N4.32 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.
Top: Total RDFs for GAP (red), DFT (blue), experimental data (black), Tersoff (cyan), ReaxFF (yellow), Billeter without added terms (violet), and BHM (green). Bottom: Partial RDFs of Si–N (top), Si–Si (middle), and N–N (bottom) for GAP (color) and DFT (shadow).
Top: Total RDFs for GAP (red), DFT (blue), experimental data (black), Tersoff (cyan), ReaxFF (yellow), Billeter without added terms (violet), and BHM (green). Bottom: Partial RDFs of Si–N (top), Si–Si (middle), and N–N (bottom) for GAP (color) and DFT (shadow).
The amorphous Si3N4 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.
Si–N bond-length distributions (top) and N–Si–N bond-angle distributions (bottom) for amorphous Si3N4 structures generated with GAP (red), DFT (blue), and BMH (green).
Si–N bond-length distributions (top) and N–Si–N bond-angle distributions (bottom) for amorphous Si3N4 structures generated with GAP (red), DFT (blue), and BMH (green).
A further validation test was performed by computing the neutron scattering factor for the GAP and DFT amorphous Si3N4 structures. The results are shown in Fig. 6, where experimental data for an amorphous Si3N4 sample with a density of 2.60 g cm−3 (Ref. 32) were added for reference purposes. We note that the density of our GAP amorphous Si3N4 was 2.92 g cm−3. 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.14
Comparison of the structure factor obtained from neutron scattering for the amorphous Si3N4 generated with GAP (red) and DFT (blue), together with experimental data found in the literature.32
Comparison of the structure factor obtained from neutron scattering for the amorphous Si3N4 generated with GAP (red) and DFT (blue), together with experimental data found in the literature.32
The main application of our ML interatomic potential is to be employed in running melt-and-quench MD to model amorphous Si3N4 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.
Coordination number statistics for Si (top) and N (bottom) in the atomic structures obtained with GAP (red) and DFT (blue).
Coordination number statistics for Si (top) and N (bottom) in the atomic structures obtained with GAP (red) and DFT (blue).
C. Vibrational density of states
In order to further test our proposed ML interatomic potential, the vibrational density of states (VDOS) of the amorphous Si3N4 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.
Comparison of the calculated VDOS to DFT results. There is a good match over the entire range of relevant energy values.
Comparison of the calculated VDOS to DFT results. There is a good match over the entire range of relevant energy values.
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.
D. Computational costs
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 Si3N4 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.
Analysis of the computational cost of the proposed ML interatomic potential. Top: Running time vs the number of atoms in the system, for a simulation of 10.000 time-steps. Bottom: Performance vs different numbers of cores for a system size of 3024 atoms.
Analysis of the computational cost of the proposed ML interatomic potential. Top: Running time vs the number of atoms in the system, for a simulation of 10.000 time-steps. Bottom: Performance vs different numbers of cores for a system size of 3024 atoms.
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 Si3N4 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.
IV. CONCLUSIONS
We have presented a machine learning (ML) interatomic potential for silicon nitride (Si3N4) 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 Si3N4 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.
SUPPLEMENTARY MATERIAL
See the supplementary material for the results of further tests comparing our proposed ML interatomic potential to experimental data, DFT, and empirical potentials.
ACKNOWLEDGMENTS
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.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
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).
DATA AVAILABILITY
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.27