LAMMPS Implementation of Rapid Artificial Neural Network Derived Interatomic Potentials

Artificial neural networks (ANNs) and other machine learning techniques are gaining prominence as a method for developing classical interatomic potentials for use in molecular statics (MS) and molecular dynamics (MD) simulation [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. Because of their ability to accurately fit large and complex data sets and because they can be rapidly optimized, ANNs have shown the potential to predict atomic energies and forces more accurately than most existing classical formalisms [6, 7]. By using density functional theory (DFT), which is relatively expensive computationally, to produce large training databases for the ANN, potentials which reproduce this data to high precision at much faster computational speeds can be created. While only a limited number of large scale dynamic simulations have been performed using these potentials, static simulations for a variety of materials have been able to reproduce ab initio forces and energies with chemical accuracy [11, 2, 12, 5].
ANNs used to model interatomic potentials behave similarly to traditional formalisms with an energy being defined for each individual atom based on the local atomic arrangement and the atomic species of it and its neighbors, mapping the local environment directly to numerical values for the energy and forces. This encoding of the local environment is referred to as the structural fingerprint. Thompson et al. developed the Spectral Neighbor Analysis Method (SNAP) [13] which uses machine-learning techniques on components of the local neighbor density. Behler and Parrinello [14, 15] defined a set of functions based on Gaussians which could be used to generate the fingerprint. This fingerprint formalism has since been used by a number of authors as the basis for an ANN. The celebrated Gaussian Approximation Potentials (GAP)[16] have been used to develop potentials with chemical accuracy for a variety of systems including tungsten [17], lithium-carbon[18], iron[19] and water[20]. Artrith and Urban used the Behler-Parrinello basis functions to create a potential for titanium dioxide using an ANN with two hidden layers with 10 neurons each which accurately described the various cystal phases of TiO₂[2]. Artrith, Urban, and Ceder have also recently proposed an efficient scheme for the generation of machine-learned potentials for composite materials with many atomic species[21]. Takahahsi et al [11] used a combination of the Gaussian formalism along with insights from the modified embedded atom method (MEAM) [22] to generate a fingerprint of several thousand values which was then used to fit a potential for titanium using linear ridge regression. Recently, a method using the embedding functions from MEAM has been used to generate potentials for titanium and zinc using a single hidden layer [23, 24]. Singraber, Behler, and Dellago [25] recently presented the n2p2 neural network package which uses the Behler symmetry functions in the construction of the fingerprint. This package has been used to develop a number of potentials including for Al-Cu [26] and Mg [27].
While the optimal structural fingerprint and network architecture for an ANN is still unclear and indeed is probably strongly dependent on the material system being examined, there is a need to be able to perform large scale dynamic calculations with these potentials. Additionally, while machine learned potentials greatly outperform ab initio calculations in terms of efficiency, they still lag behind conventional formalisms such as the embedded atom method (EAM)[28] and MEAM, which have been the standard for analysis of metals since their introduction over 25 years ago. Additionally, MEAM takes advantage of angular screening to limit the size of the neighbor list and incorporate the known shielding effect to improve performance. To this end, we have created an implementation of a rapid artificial neural network (RANN) potential formalism for use in LAMMPS [29]. This implementation accepts arbitrary Multi-layer Perceptron (MLP) architectures and structural fingerprints of arbitrary size. Sample potentials have been generated from the MEAM based structural fingerprints with angular screening, though novel functions or fingerprint types can be added independently. Section 2 of this paper reviews the formalism for ANN interatomic potentials and provides the specific structural fingerprints implemented including a discussion of the use of angular screening. Section 3 describes the framework which has been implemented in LAMMPS. Section 4 demonstrates a number of simple static and dynamic calculations using an ANN potential. In particular, conservation of energy is shown along with accurate calculation of the forces and pressure, and a comparison is done between the computational efficiency of this potential style as compared to other common styles for metal systems. It is demonstrated that, for the network architectures capable of reproducing a magnesium training database to within 1meV/atom, performance approximately one third as MEAM can be achieved. Additionally, the performance of the potential as a function of neural network architecture and neighbor list size is considered.

Most machine-learned interatomic potentials make use of MLPs to predict the energy of individual atoms in a system given their environment. Details of the training and optimization of these networks in reference to ab initio data can be found elsewhere [2, 23, 30], but the general structure of a network as an interatomic potential is as follows. The energy of a particular atom $i$, determined by its environment, is the last of N layers of the neural network. The values for any particular layer, $\mathbf{A}^{n}$, after the first is determined by the previous layer and the weight and bias matrices $\mathbf{W}^{n}$ and $\mathbf{B}^{n}$:

	$\displaystyle Z^{n}_{l_{n}}=\sum_{l_{n-1}}{W^{n}_{l_{n}l_{n-1}}A^{n-1}_{l_{n-1}}+B^{n}_{l_{n}}}$		(1)
	$\displaystyle A^{n}_{l_{n}}=g^{n}(Z^{n}_{l_{n}})$		(2)

Where $l_{n}$ is the number of neurons in layer $n$ and $g^{n}(x)$ is a nonlinear activation function. The input layer, $\mathbf{A}^{0}$ is given by the structural fingerprint of the local atomic environment. The ouput layer, $\mathbf{Z}^{N}$, will always contain a single node, representing the energy, and so $\mathbf{W}^{N}$ will always be a row vector and $\mathbf{B}^{N}$ always a single number.
A number of different functions can be used to describe this structural fingerprint. Behler and Parrinello [15] consider two types of fingerprint functions: radial symmetry functions which are a pairwise sum of Gaussians, and angular terms over all triplets of atoms. For the RANN style, however, we use the fingerprint style introduced by Dickel, Francis, and Barrett, motivated by the Modified Embedded Atom Method (MEAM) formalism [22, 31, 32], with the addition of angular screening. This style seems particularly relevant for metal potentials for the same physically motivated reasons that the MEAM formalism is effective. In this style, two different kind of input fingerprints are considered. First, simple pair interactions are considered and summed over all the neighbors of a given atom. For a given atom labeled $i$, we define a set of pair potentials interactions with the form

$$F_{n}=\sum_{j\neq i}{(\frac{r_{ij}}{r_{e}})^{n}e^{-\alpha_{n}\frac{r_{ij}}{r_{e}}}f_{c}(\frac{r_{c}-r_{ij}}{\Delta r})S_{ij}}$$

(3)

Where $j$ labels all the neighbors atoms of $i$ within a cutoff radius $r_{c}$, $n$ is an integer, different for each member of the pairwise contributions to the fingerprint, $r_{e}$ is the equilibrium nearest neighbor distance, $S_{ij}$ is an angular screening term described below, and $\alpha_{n}$ are metaparameters, which can be tuned to better optimize the potential. In principle, $\alpha_{n}$ are related to the dimensionless $\alpha$ used in the MEAM formalism which is related to the bulk modulus of the material [32]. $f_{c}(x)$ is a cutoff function which smoothly varies the weight of the interaction from 1 for atoms inside the cutoff radius to 0 once $r_{ij}>r_{c}$ where $\Delta r$ defines the width of the transition. The cutoff function used here is the same as employed in MEAM and is given by

$$r_{c}(x)=\left\{\begin{array}[]{lr}1,&x>1\\ (1-(1-x)^{4})^{2},&0\leq x\leq 1\\ 0,&x<0\end{array}\right\}$$

(4)

The second kind of fingerprint function considers three body terms, with a form similar to the partial electron densities used in MEAM. There are two parameters, $m$ and $k$ which can be varied to generate more inputs for the fingerprint:

$$G_{m,k}=\sum_{j,k}{\cos^{m}{\theta_{jik}}e^{-\beta_{k}\frac{r_{ij}+r_{ik}}{r_{e}}}f_{c}(\frac{r_{c}-r_{ij}}{\Delta r})f_{c}(\frac{r_{c}-r_{ik}}{\Delta r})S_{ij}S_{ik}}$$

(5)

where $\theta_{jik}$ is the angle between $r_{ij}$ and $r_{ik}$, m is a non-negative interger and $\beta_{k}$ is a set of metaparameters which determine the length scale of the different terms. For large numbers of neighbors, it can be more efficient to calculate this as a sum over a single list of neighbors as follows.

	$\displaystyle G_{0,k}=\left(\sum_{j}{e^{-\beta_{k}\frac{r_{ij}}{r_{e}}}f_{c}(r_{ij}})S_{ij}\right)^{2}$		(6)
	$\displaystyle G_{1,k}=\sum_{\alpha_{1}}\left(\sum_{j}\frac{r^{\alpha_{1}}_{i,j}}{r_{i,j}}e^{-\beta_{k}\frac{r_{ij}}{r_{e}}}f_{c}(\frac{r_{c}-r_{ij}}{\Delta r})S_{ij}\right)^{2}$		(7)
	$\displaystyle G_{2,k}=\sum_{\alpha_{1},\alpha_{2}}\left(\sum_{j}\frac{r^{\alpha_{1}}_{i,j}r^{\alpha_{2}}_{i,j}}{r_{i,j}^{2}}e^{-\beta_{k}\frac{r_{ij}}{r_{e}}}f_{c}(\frac{r_{c}-r_{ij}}{\Delta r})S_{ij}\right)^{2}$		(8)
	$\displaystyle G_{3,k}=\sum_{\alpha_{1},\alpha_{2},\alpha_{3}}\left(\sum_{j}\frac{r^{\alpha_{1}}_{i,j}r^{\alpha_{2}}_{i,j}r^{\alpha_{3}}_{i,j}}{r_{i,j}^{3}}e^{-\beta_{k}\frac{r_{ij}}{r_{e}}}f_{c}(\frac{r_{c}-r_{ij}}{\Delta r})S_{ij}\right)^{2}$		(9)
	$\displaystyle G_{m,k}=\sum_{\alpha_{1}}\sum_{\alpha_{2}}...\sum_{\alpha_{m}}\left(\sum_{j}\frac{r^{\alpha_{1}}_{i,j}r^{\alpha_{2}}_{i,j}...r^{\alpha_{m}}_{i,j}}{r_{i,j}^{m}}e^{-\beta_{k}\frac{r_{ij}}{r_{e}}}f_{c}(\frac{r_{c}-r_{ij}}{\Delta r}))S_{ij}\right)^{2}$		(10)

where $r^{\alpha_{m}}_{ij}$ is the $x$, $y$, or $z$ component of $r_{ij}$ for $\alpha_{m}=1,2$ or $3$, respectively. A simplification can be made by noting that many of the terms in the three body interactions are redundant. In $G_{2,k}$ for example, there will be a calculation for xy as well as for yx. By taking advantage of this redundancy, we can reduce the number of terms that need to be calculated for $G_{m,k}$ from $3^{n}$ to $\frac{(n+1)(n+2)}{2}$. The simplified expression for $G_{m,k}$ now reads:

$$G_{m,k}=\sum_{\alpha_{1}\geq\alpha_{2}\geq...\geq\alpha_{m}}\left(\frac{n!}{n_{x}!n_{y}!n_{z}!}\sum_{j}{\frac{r^{\alpha_{1}}_{i,j}r^{\alpha_{2}}_{i,j}...r^{\alpha_{m}}_{i,j}}{r_{i,j}^{m}}e^{-\beta_{k}\frac{r_{ij}}{r_{e}}}f_{c}(\frac{r_{c}-r_{ij}}{\Delta r})}\right)^{2}$$

(11)

where $n_{x}$, $n_{y}$, and $n_{z}$ are the total number of x-, y-, and z-components in the set $(\alpha_{1},\alpha_{2},...,\alpha_{m})$, respectively. Note that which of these forms (Eq. 5 or Eq.11 will be more efficiently calculated will depend on the length of the neighbor list and the magnitude of $m$. Structural fingerprints of arbitrary size can be constructed by varying the number of different $n$, $m$, and $k$ used giving more or less information to the ANN. For the potentials tested here, we have used consecutive values for n and m and totaling $n_{tot}$ and $m_{tot}$ respectively ($n\in(-1,0,1,...,n_{tot}-2)$, $m\in(0,1,...,m_{tot}-1)$). A particular ANN potential will consist of a fixed structural fingerprint, number and length of each hidden layer, weight and bias matrices for each layer and activation functions for each layer. For the potentials considered here we have used the following activation functions:

	$\displaystyle g^{N}(x)=x$		(12)
	$\displaystyle g^{n}(x)=\frac{x}{10}+\frac{9}{10}\log(e^{x}+1)\;\textrm{for}\;n<N$		(13)

The size of the weight and bias matrices will depend on the length of the fingerprint and the length of each hidden layer.
As discussed above, one of the major computational barriers in the implementation of ANN methods is the large neighbor list included in the calculation of the structural fingerprint. In order to obtain good convergence with DFT results, cutoff radii used can be larger than 12 Å[27], which for the ground state of Mg includes over 300 atoms. By comparison, recent MEAM potentials for Mg [33] have used a cutoff distance of 5.875 Å, which includes only 38 atoms in the ground state. Since the computation time scales at least linearly with the number of neighbors, minimizing the length of the neighbor list without sacrificing performance should be a key goal for optimal efficiency. The radial screening function $f_{c}$ can accomplish this, but leads to nonphysical results at large separations (see, for example [33]). Angular screening, whereby the effective interaction between atoms is reduced or eliminated by the presence of an atom located between them can effectively limit the neighbor list in a similar way without the unphysical results. Such a method has been utilized by MEAM [34], and the same screening method has been employed here. Briefly, the screening between two atoms is determined from the product of all the screening interactions by other atoms in the neighborhood:

$$S_{ij}=\prod_{k\neq i,j}S_{ikj}$$

(14)

where $S_{ikj}$ is calculated from a geometric construction considering the ellipse formed by the 3 atoms with $r_{i,j}$ one of the axes. The screening parameter, $C_{ikj}$ is then given by:

$$C_{ikj}=1+2\frac{r^{2}_{ij}r^{2}_{ik}+r^{2}_{ij}r^{2}_{jk}-r^{4}_{ij}}{r^{4}_{ij}-(r^{2}_{ik}-r^{2}_{jk})^{2}}$$

(15)

and then screening value is

$$S_{ikj}=f_{c}\left(\frac{C_{ikj}-C_{min}}{C_{max}-C_{min}}\right)$$

(16)

where $f_{c}$ is the same cutoff function used for the radial cutoff and $C_{max}$ and $C_{min}$ are metaparameters which can be tuned to determine which neighbors can be excluded from calculations.
The effect of including angular screening can be demonstrating by considering the change of an individual fingerprint as the length scale is changed continuously. For this example, we consider the value of a single fingerprint for a Mg atom in a perfect bulk hcp lattice as the lattice parameter is changed continuously. In the absence of angular screening, the value of the fingerprint will change rapidly at particular values of the lattice constant as new neighbors enter the radial screening distance. If angular screening is included with metaparameters such that, for example, only 3rd nearest neighbors are ever included regardless of lattice parameter, the change in the value is considerably smoother. Figure 1 demonstrates the difference between the two cases. Not only is the unscreened value more expensive to calculate as it requires more neighbors for most cases, it can also be seen that the predictive power of the model is hampered as the relation between two states which are physically similar – ideal 3D lattices with different lattice constants, given markedly different trends and values as the number of neighbors changes. This can be overcome through the action of the neural network and precise cancellation among different fingerprint terms, but the quality of fit and the extrapolative power of the screened functions appears more natural. Indeed unphysical deviations at relatively moderate isotropic compressions, as seen for example in [27] can be avoided entirely through angular screening as the fingerprint values change continuously in the screened case even at high compression.

For the full network, forces are calculated by taking a gradient of the energy over the entire system. This includes the individual ANN contributions from every atom in the system, so the force on atom $i$ is determined not only by the ANN output of that atom, but of all of its neighbors as well.

The ANN potential file written for LAMMPS includes all of the details of the neural network. It begins with a header which lists the types and numbers of fingerprints to be used. While we have only included the two varieties, the bond power and radial power terms, $F_{n}$ and $G_{m,k}$, in this work, the LAMMPS potential style can be expanded in a straightforward way to include other descriptors. For each of the types used, the potential file then specifies the number used and any metaparameters. For $F_{n}$ this includes the range of values for $n$ as well as the metaparameters $\alpha_{n}$, $r_{c}$, $C_{max}$, $C_{min}$ and $r_{e}$. For $G_{m,k}$, $m_{tot}$ and $k_{tot}$ are given along with the metaparameters $\beta_{k}$, $r_{e}$, $r_{c}$, $C_{min}$, and $C_{max}$. Note that while we use the same cutoff radius and angular screening parameters for both varieties, they can in principle be different as every fingerprint type has its own set of metaparameters. Next, the network architecture is specified by giving the total number of layers and the length of each layer. Figure 2 shows a sample header for the potential file. Following the header is the weight and bias matrices for each layer and finally the activation functions used for each layer. Using the potential file, the ANN subroutine then calculates the feature list for every atom and uses these to calculate the energy and force for every atom in the usual way. The subroutine and potential file also have the capability of assigning multiple ANNs to different atoms types to be used for multi-element simulations. In this way, the framework can be used for arbitrary structural fingerprints as inputs to arbitrary MLP architectures for systems with an arbitrary number of element types.
The subroutine is included as a normal pair-style in LAMMPS called ‘rann’ and can be called in the usual way:

With the pair style implemented in LAMMPS, several potentials were fit for a reference database for magnesium (Mg). The database, containing over 300,000 atomic environments was constructed using the open-source density functional theory (DFT) software, Quantum Espresso [35]. The configurations used to train the networks are shown in Table 1. These configurations included simple triaxial strain perturbations of the low energy SC, FCC, BCC, and HCP phases of Mg along with larger supercells with atoms displaced from equilibrium to mimic the effect of thermal excitation. Free surface data and vacancy data, with atoms also perturbed from equilibrium positions were also included to improve stability. Effective temperatures up to 1000K were included in the thermally excited database. Training of the network was carried out by minimizing the RMSE of the system energies compared to the DFT results using the Levenerg-Marquadt method [36, 37]. 10% of the training data from each configuration type was left out of training and used for validation. The most significant variation however, came in the size and construction of the fingerprint used as the input layer. Along with the number and type of fingerprints included, the most important factor in determining accuracy and efficiency was the number of atoms included in the neighbor list. This is a function of the cutoff radius $r_{c}$ and the angular screening parameter $C_{min}$. Table 2 summarizes the different fingerprints considered and the RMSE for both the training and validation datasets for each potential. As can be seen, even the smallest architectures had RMSE of less than 3 meV/atom and the best were less than 1 meV/atom. Table 2 also includes computational performance for a sample calculation which will be described in detail below.

We have presented an implementation of the Rapid Artificial Neural Network (RANN) interatomic potential pair style for use in the LAMMPS molecular dynamics code. The potential is suitable for static and dynamic calculations, conserving total energy and correctly calculating the pressure and forces of bulk solids. Such potential styles have been shown to be capable of reproducing training data from first principles calculations at a level which exceeds previous semi-empirical formalisms [6, 7, 10]. The developed pair style is capable of accommodating network architectures of arbitrary dimension and activation function, and the calculation of the structural fingerprint requires only singly looped summations over the neighbor list of a target atom, improving computational efficiency. Angular screening has also been introduced both to improve efficiency by limiting the number of included neighbors and to improve predictive power by introducing the physically motivated phenomenon of shielding and smoothing fingerprints as atoms move through the radial cutoff. Computational efficiency is of particular importance for this formalism as the improvements in accuracy over existing formalisms such as MEAM are greatly diminished if the runtime is closer to that of first principles calculations. As such, we demonstrate that for a network with 13 input fingerprints and a cutoff radius of 6.0 $\AA$ for the neighbor list, and angular screening restricting interactions to third nearest neighbors in bulk HCP the implementation performs at one third the speed of MEAM. Larger networks, fingerprints, or neighbor lists will hamper this performance, with the most significant cost associated with increasing the neighbor list through the length of the radial cutoff. Ultimately, a flexible, scalable formalism for ANNs, which can utilize the native parallel processing available in LAMMPS has been demonstrated with a formalism which produces high accuracy potentials which operate on speeds comparable to those of MEAM.