Efficient and Accurate Machine-Learning Interpolation of Atomic Energies in Compositions with Many Species

TL;DR

This paper introduces a simple, efficient descriptor for machine-learning potentials that scales independently of the number of chemical species, enabling accurate simulations of complex materials with many elements.

Contribution

The authors demonstrate that a constant-complexity descriptor suffices for accurate ML potentials in multi-species systems, overcoming previous quadratic scaling limitations.

Findings

Achieved around 3 meV/atom accuracy for transition-metal oxides and biomolecules with 11 species.

Proved that quadratic scaling of descriptors is unnecessary for complex compositions.

Enabled ML simulations of previously inaccessible multi-species materials.

Abstract

Machine-learning potentials (MLPs) for atomistic simulations are a promising alternative to conventional classical potentials. Current approaches rely on descriptors of the local atomic environment with dimensions that increase quadratically with the number of chemical species. In this article, we demonstrate that such a scaling can be avoided in practice. We show that a mathematically simple and computationally efficient descriptor with constant complexity is sufficient to represent transition-metal oxide compositions and biomolecules containing 11 chemical species with a precision of around 3 meV/atom. This insight removes a perceived bound on the utility of MLPs and paves the way to investigate the physics of previously inaccessible materials with more than ten chemical species.

Equations

$$E(\sigma )=\sum _i^{\text{atoms}}{\text{MLP}}_{t_i}({\tilde{\sigma }}_i^{R_{\text{c}}})$$

$${\mathrm{RDF}}_i(r)$$

$${\mathrm{ADF}}_i(\theta )$$

$${\mathrm{RDF}}_i(r)$$

$${\mathrm{ADF}}_i(\theta )$$

$$c_{\alpha }^{(2)}$$

$$c_{\alpha }^{(3)}$$

$${\mathrm{RDF}}_i(r)=\sum _{\alpha }{c}_{\alpha }^{(2)}{\varphi }_{\alpha }(r),$$

$$c_{\alpha }^{(2)}={\int }_0^{R_{\text{c}}}{\bar{\varphi }}_{\alpha }(r){\mathrm{RDF}}_i(r)\mathrm{d}r.$$

$${\mathrm{RDF}}_i(r)$$

$$c_{\alpha }^{(2)}$$

$$=\sum _{{\mathbf{R}}_j\in {\sigma }_i^{R_{\text{c}}}}{\bar{\varphi }}_{\alpha }(R_{ij})f_{\text{c}}(R_{ij})w_{t_j},$$

$$\begin{array}{rl}{T}_0(x)& =1\\ T_1(x)& =x\\ T_{n+1}(x)& =2x{T}_n(x)-T_{n-1}(x).\end{array}$$

$${\int }_{-1}^1{T}_n(x)T_m(x)\frac{\mathrm{d}x}{\sqrt{1-x^2}}=\{\begin{array}{ll}\pi & n=m=0\\ \frac{\pi }{2}& n=m\ne 0\\ 0& \text{else}\end{array}$$

$${\varphi }_{\alpha }(r)$$

$$\text{and}{\bar{\varphi }}_{\alpha }(r)$$

Full text

Atomic interaction potentials based on the interpolation of first-principles calculations with machine-learning algorithms have the potential to enable efficient linear-scaling atomistic simulations with an accuracy that is close to the reference method Lorenz et al. (2004); Behler and Parrinello (2007); Bartók et al. (2010); Rupp et al. (2012). Such machine-learning potentials (MLPs) establish a relationship between a unique descriptor and the total or atomic energy using, e.g., artificial neural networks (ANNs) Montavon et al. (2012) or Gaussian process regression (GPR, Kriging) Rasmussen and Williams (2006). However, the combined space of atomic coordinates and chemical species grows rapidly with the number of chemical species, resulting in a formal corresponding growth of the descriptor complexity and thus the complexity of the MLP. This scaling has so far limited current MLP approaches to compositions with only a few chemical species Artrith et al. (2011); Artrith and Kolpak (2014, 2015); Morawietz et al. (2016) or atomic structures Faber et al. (2016). Overcoming this limitation is a very active field of research Huo and Rupp (2017); De et al. (2016).

In this article we demonstrate that the computational complexity of MLPs does not necessarily grow with the number of chemical species, so that MLPs for materials with ten or more chemical species are in principle feasible and computationally efficient. We show that, contrary to intuition and common belief, the same model complexity that is optimal for a ternary material is also sufficient to describe a system with 11 chemical species (Fig. 1). To illustrate these concepts, we consider two different materials classes of practical relevance: cation-disordered lithium transition-metal (TM) oxides, which have recently attracted interest as high-energy-density cathode materials for Li-ion batteries Lee et al. (2014); Yabuuchi et al. (2015), and proteinogenic amino acids, i.e., the building blocks of proteins and their complexes with divalent cations Ropo et al. (2016a, b). We show that both of these high-dimensional materials systems can be accurately modeled using MLPs based on a mathematically simple and computationally efficient descriptor with constant complexity that we will introduce in the following.

In the present work, we focus on MLPs that express the total structural energy as the sum of atomic energy contributions and are in this respect similar to other many-body potentials such as embedded atom models Daw and Baskes (1984); Daw et al. (1993). However, unlike conventional potentials, the atomic energy is not confined to a rigid functional form, but is represented by a flexible non-linear machine-learning model that is trained to a descriptor of the local atomic environment. In this context, the local atomic environment $\sigma_{i}^{R_{\textup{c}}}\subset\sigma$ of an atom $i$ in a structure $\sigma$ is defined as the local structure given by the set of coordinates $\{\mathbf{R}_{1},\mathbf{R}_{2},\ldots\}$ of all atoms within a cutoff distance $R_{\textup{c}}$ from atom $i$ and the local composition, i.e., the corresponding chemical species $\{t_{1},t_{2},\ldots\}$. To be physically meaningful and transferable between equivalent structures, the descriptor needs to be invariant with respect to translation and rotation of the structure and the exchange of equivalent atoms. Several transformations for $\sigma_{i}^{R_{\textup{c}}}$ into invariant representations $\widetilde{\sigma}_{i}^{R_{\textup{c}}}$ have been proposed in the literature Behler (2011); Bartók et al. (2013); Sadeghi et al. (2013); Schütt et al. (2014); Behler (2015); Faber et al. (2015); von Lilienfeld et al. (2015), and the most commonly used methods for MLPs are the symmetry functions by Behler and Parrinello (BP) Behler and Parrinello (2007); Behler (2011) and the smooth overlap of atomic positions (SOAP) approach by Bartók, Kondor, and Csányi Bartók et al. (2013); Bartók and Csányi (2015); De et al. (2016). With an invariant descriptor $\widetilde{\sigma}_{i}^{R_{\textup{c}}}$, the total MLP energy of a structure $\sigma$ can then be expressed as

[TABLE]

.

Our approach draws inspiration from the strength of the established descriptor methods but explicitly maintains the distinction between local structure and composition by using two sets of invariant coordinates, ${}^{\{\mathbf{R}\}}\widetilde{\sigma}_{i}^{R_{\textup{c}}}$ and ${}^{\{t\}}\widetilde{\sigma}^{R_{\textup{c}}}_{i}$, that separately encode the atomic positions and species. The union of both sets, $\widetilde{\sigma}_{i}^{R_{\textup{c}}}={{}^{\{\mathbf{R}\}}\widetilde{\sigma}}_{i}^{R_{\textup{c}}}\cup{{}^{\{t\}}\widetilde{\sigma}}^{R_{\textup{c}}}_{i}$, is used as a combined descriptor for an ANN-based MLP (ANN potential). As structural descriptor ${{}^{\{\mathbf{R}\}}\widetilde{\sigma}}_{i}^{R_{\textup{c}}}$ we choose the expansion coefficients of the radial (bond length) and angular (bond angle) distribution functions in a complete basis set $\{\phi_{\alpha}\}$,

[TABLE]

and the compositional descriptor ${{}^{\{t\}}\widetilde{\sigma}}^{R_{\textup{c}}}_{i}$ is given by the expansion coefficients of the same distribution functions but with atomic contributions that are weighted differently for each chemical species. The RDF and ADF obey the invariants of the atomic energy, and basing the descriptor on an expansion in a complete basis set allows its systematic refinement by converging the number of basis functions. We implemented the descriptor into the free and open-source atomic energy network package Artrith and Urban (2016).

In general, multi-layer ANNs can reproduce any function with arbitrary precision Hornik (1991). However, the resolution of the invariant descriptor determines the maximal precision with which an ANN potential can resolve the chemical space of a given material. To determine the resolution of our combined descriptor, we trained ANN potentials to extensive reference data sets with different numbers of chemical species. We consider the resolution satisfactory if the ANN potential can reproduce the reference energies of our data sets with a precision of $\sim$3 meV/atom, which is the order of magnitude of the noise in our reference data.

Figure 1a shows the precision that can be achieved in representing Li-TM oxides with different numbers of TM species using ANN potentials based on the combined descriptor with different numbers of basis functions. The reference set for the ANN potential training comprised Hubbard-U corrected Liechtenstein et al. (1995); Anisimov et al. (1997); Jain et al. (2011) density-functional theory (DFT) energies and optimized structures of 16,047 \ceLi$M$O2 configurations in the rocksalt structure with different compositions based on nine TMs (Sc, Ti, V, Cr, Mn, Fe, Co, Ni, and Cu) and cation arrangements with up to 36 atoms. For all DFT+U calculations we employed the PBE exchange-correlation functional Perdew et al. (1996) with projector-augmented wave Blöchl (1994) pseudopotentials as implemented in VASP Kresse and Furthmüller (1996a, b). DFT energies and atomic forces were converged to 0.05 meV per atom and 50 meV/Å, respectively, gamma-centered k-point meshes with a density of 1000 divided by the number of atoms were used, and the plane-wave cutoff was 520 eV. VASP input files were generated using the pymatgen software with default parameters Ong et al. (2013). Structures with up to 5 chemical species were generated by systematic enumeration, and random atomic configurations were generated for compositions with 6–11 chemical species. Further information about the generation of these reference structures, the parameters of our DFT calculations, and the architecture of the ANNs are given in the Appendix.

As seen in Fig. 1a, the ANN potentials achieve a root mean squared error (RMSE) of $\sim$3 meV/atom relative to the DFT reference energies with a descriptor dimension of 44 (i.e., 22 basis functions). Note that, for the present work, we employed the same number of basis functions for the radial and angular expansion (i.e., 11 each), though this is not a general requirement of the methodology. Increasing the descriptor dimension beyond 52 or 60 results in a minor additional reduction of the RMSE at the cost of significantly increased computational effort. We emphasize that this RMSE is purely a quality measure of the descriptor precision and does not reflect the accuracy of the ANN potentials in simulations, which would have to be carefully validated separately.

The RMSE was evaluated after 3,000 training iterations using the LM-BFGS method Byrd et al. (1995); Zhu et al. (1997), however, with increasing number of species and increasing descriptor size the required number of training iterations to achieve convergence generally also increases. Thus, the ANN potentials for 11 chemical species and descriptor dimensions above 40 have not converged after 3,000 iterations, and the RMSEs after 5,000 iterations are shown in Fig. 1. The unconverged RMSE after 3,000 training iterations is shown in Fig. S2 in the Appendix.

Remarkably, the optimal descriptor dimension is essentially independent of the number of chemical species in the composition, and a descriptor dimension of 44 is sufficient to capture the structural and chemical features of the distinct atomic configurations in the \ceLi$M$O2 data set with up to 11 chemical species.

Figure 1b shows the equivalent analysis for the first-principles energies and structures of 45,892 conformations of the proteinogenic amino acids (5 chemical species: H, C, N, O, and S) and their complexes with the six divalent cations \ceBa^2+, \ceCa^2+, \ceCd^2+, \ceHg^2+, \cePb^2+, \ceSr^2+ (a total of 11 chemical species) by Ropo, Schneider, Baldauf, and Blum Ropo et al. (2016a) based on DFT calculations (PBE+TS-vdW Tkatchenko and Scheffler (2009)) using the FHI-aims package Blum et al. (2009). This data set was compiled specifically for the parametrization of atomic potentials and thoroughly samples the relevant conformational space Ropo et al. (2016a), an important first step towards improved force fields for proteins Piana et al. (2014). The high precision of the ANN potentials with an RMSE of $\sim$3 meV/atom for 5 and 11 chemical species indicates that our combined descriptor is not limited to crystal structures with similar atomic positions, but is also suitable to distinguish between continuous atomic arrangements.

To understand the significance of these observations, we first describe the details of the structural and compositional descriptor. We begin by expressing the atom-centered radial and angular distribution functions of Eqs. (7) and (2) in terms of discrete delta functions centered at the bond lengths between atoms $j$ and the central atom $i$, $R_{ij}=||\mathbf{R}_{j}-\mathbf{R}_{i}||$, and the bond angle $\theta_{ijk}=\angle(\mathbf{R}_{j}-\mathbf{R}_{i},\mathbf{R}_{k}-\mathbf{R}_{i})$

[TABLE]

where $f_{\textup{c}}$ is a cutoff function that smoothly goes to zero at $R_{\textup{c}}$ (in practice, we use $f_{\textup{c}}(r)=0.5[\cos(r\cdot\pi/R_{\textup{c}})+1]$). The weights $w_{t_{j}}$ and $w_{t_{k}}$ are 1 for the structural descriptor ${}^{\{\mathbf{R}\}}\widetilde{\sigma}_{i}^{R_{\textup{c}}}$ and take on species-dependent values for the compositional descriptor ${}^{\{t\}}\widetilde{\sigma}_{i}^{R_{\textup{c}}}$. Here, we followed the (Ising-model) pseudo-spin convention commonly used for lattice models Sanchez et al. (1984), i.e., $w_{l}=0,\pm{}1,\pm{}2,\ldots$ where 0 is omitted for even numbers of species. For the expansions Eqs. (7) and (2) we choose a complete orthonormal basis $\{\phi_{\alpha}\}$, i.e., $\int\phi_{\alpha}\overline{\phi}_{\beta}=1$ if $\alpha=\beta$ and [math] else. With this choice, the expansion coefficients are given by

[TABLE]

A derivation of Eqs. (10) and (6) can be found in the Appendix. The expansions are truncated at finite radial and angular orders $N_{2}$ and $N_{3}$ that determine the dimension (i.e., the complexity) and the resolution of the descriptor, i.e., ${}^{\{\mathbf{R}\}}\widetilde{\sigma}_{i}^{R_{\textup{c}}}=\{{{}^{\{\mathbf{R}\}}c^{(2)}_{1}},\ldots,{{}^{\{\mathbf{R}\}}c^{(2)}_{N_{2}}},{{}^{\{\mathbf{R}\}}c^{(3)}_{1}},\ldots,{{}^{\{\mathbf{R}\}}c^{(3)}_{N_{3}}}\}$.

For this article, we employed the Chebyshev polynomials of the first kind as basis functions (see Appendix), as they can be defined in terms of a recurrence relation that allows for highly efficient numerical evaluation of the function values and their derivatives. With this choice of basis functions, Fig. 2 shows the RDF as reconstructed based on the structural expansion coefficients $\{^{\{\mathbf{R}\}}c_{\alpha}^{(2)}\}$ for two different orders ($N_{2}=50$ and $N_{2}=150$). From comparison with Gaussian convolutions of the discrete RDF, the radial resolution of the expansion order $N_{2}=150$ is around 0.1 Å. Atomic features on smaller scales may affect the shape of the RDF but do not give rise to distinct peaks. The expansion of the ADF is completely analogous.

We note that the radial and angular BP symmetry functions Behler and Parrinello (2007); Behler (2011) can be cast into the form of Eqs. (10) and (6) but are neither orthogonal nor systematically refinable. The relationship of our structural descriptor to the coefficients of a basis set expansion is, in turn, closer in spirit to the SOAP method Bartók et al. (2010, 2013) which is based on the power spectrum of the atomic density of the local atomic environment. SOAP allows for a rigorous and systematic description of the local structure, which comes at the cost of an arithmetically (and computationally) more complex formalism. However, by limiting the descriptor to radial and angular contributions our method maintains the simple analytic nature of the BP approach that allows for a highly efficient numerical implementation and straightforward differentiation (which is required for the calculation of analytic forces and higher derivatives). Basing the radial and angular descriptors on an expansion in a complete basis set allows their systematic refinement in the spirit of the SOAP approach, though our approach is limited to two- and three-body interactions.

Also note that decomposing the local atomic environment into $n$-body contributions as done in our structural descriptor is an established and well-tested approach for lattice models such as the cluster expansion (CE) method Fontaine (1994); Ceder (1993). In CE models, the total configurational energy is expanded in a basis set consisting of site clusters $(s_{i}s_{j}\ldots)$ with increasing numbers of lattice sites $s_{\alpha}$, i.e., point clusters, pairs, trimers, …, $n$-tuples. The site clusters form a complete basis set, and the configurational averages of all equivalent clusters (the cluster correlations) are the descriptor of the CE model. Unlike MLPs, the CE energy is a linear function of the descriptor. For the case of the continuous structural energy, Thompson et al. demonstrated that a linear potential based on SOAP (which also is a complete basis of the local structure) can achieve reasonable accuracy in practice if a sufficient number of basis functions is used Thompson et al. (2014).

However, the strength of non-linear machine-learning models is that they do not require mathematically complete descriptors as long as the descriptor is able to differentiate between all relevant samples. This property is exploited, for example, in the area of image recognition and text classification Rogati and Yang (2002). In practice this means that even an incomplete descriptor of the local atomic environment may be sufficient to construct a non-linear MLP if that descriptor is able to differentiate between all relevant local atomic structures, i.e., the descriptor does not have to resolve all hypothetically possible sets of three dimensional coordinates.

This behavior is exemplified in Fig. 3, which compares the precision of an ANN potential for the \ceLi$M$O2 data set with 5 chemical species (10,175 atomic configurations) with that of a linear energy model as function of the descriptor dimension. As seen in the figure, the ANN achieves an RMSE of $\leq$3 meV/atom with descriptor dimensions of 44 (22 basis functions) and larger. Comparison with Fig. 2 shows that such a small basis set corresponds to a coarse representation of the RDF and the ADF, however, obviously this level of approximation is sufficient for the ANN potential to differentiate between all structural and compositional features in the reference set. This is not the case for the linear model whose RMSE is $>$15 meV/atom even for a descriptor dimension of 84.

In conclusion, we showed that machine-learning potentials do not require (mathematically) complete descriptors of the local atomic environment to reproduce potential energy surfaces with high precision. With this insight, we devised a combined descriptor of the local atomic structure and composition whose complexity does not scale with the number of chemical species. The method is conceptually simple and allows for highly efficient numerical implementations. The utility of the approach was demonstrated for two exemplary materials classes, lithium transition-metal oxides and amino acid complexes, each separately comprising compositions with 11 different chemical species. We showed that the potential energy landscape of both example systems can be represented with high precision by artificial neural network potentials using the combined descriptor achieving a resolution of around 3 meV/atom. Hence, machine-learning potentials are in practice not limited to compositions with small numbers of chemical species as previously argued in the literature and may be effective for the modeling of high-dimensional materials such as oxide solid solutions and peptide chains.

I Acknowledgments

This work was supported by the Office of Naval Research (ONR) under ONR award N00014-14-1-0444. This work used mainly the computational facilities of the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant no. ACI-1053575. Additional computational resources from the University of California Berkeley, HPC Cluster (SAVIO) are also gratefully acknowledged.

Appendix A The $\textup{Li}M\textup{O}_{2}$ data set

Starting point for the generation of the \ceLi$M$O2 data set of this work were the enumerated lithium transition-metal (TM) oxide configurations of reference Urban et al., 2016. The data set with 3 chemical species (Li, Ti, and O) comprised a total of 7,338 structures including the \ceTiO2 structures from reference Artrith and Urban, 2016 and additional \ceLiTiO2 configurations that were systematically enumerated based on the rocksalt structure (A sites = Li and Ti, B sites = O) up to cell sizes containing 8 cation sites using the approach by Hart et al. Hart and Forcade (2008, 2009); Hart et al. (2012). For the data set with 4 chemical species (Li, Ni, Ti, and O), 1,343 atomic configurations with compositions \ceLiNiO2 and \ceLi2NiTiO4 were additionally generated using the same enumeration methodology (giving a total of 8,681 configurations). Further, 1,494 atomic configurations with compositions \ceLiMnO2 and \ceLi2NiMnO4 were added for the set with 5 species (Li, Ti, Mn, Ni, O) to a total of 10,175 configurations. Finally, for 11 chemical species, random atomic configurations with composition \ceLi9$M$9O18 with $M$ = Sc, Ti, V, Cr, Mn, Fe, Co, Ni, Cu for all 24,310 possible compositions were generated, and 5,872 randomly selected configurations were included in the reference data set.

The complete reference data set for 11 chemical species comprises a total of 16,047 atomic configurations.

Appendix B Artificial neural network potentials

Together with the dimension of the descriptor discussed in the main text, the architecture of an artificial neural network (ANN) determines the model complexity. For the feedforward ANNs used in the present work, the architecture is given by the number of hidden layers and the number of nodes per layer employed by the ANN (see also reference Artrith and Urban (2016) for a detailed introduction). To facilitate comparison of the different structure and composition spaces on equal footing, we generally used a $N$-15-15-1 ANN architecture, i.e., an architecture with two hidden layeres containing each 15 nodes independent of the descriptor dimension $N$.

The atomic energy network (aenet) software Artrith and Urban (2016) was employed for the training of the ANN potentials with the limited-memory BFGS method Byrd et al. (1995); Zhu et al. (1997).

Appendix C Scaling behavior of existing local descriptors

The complexity of the Behler–Parrinello (BP) descriptor and the cluster-expansion basis scales at least quadratically with the number of chemical species.

As noted in the main manuscript, separate MLPs for each chemical species are constructed. For each of these MLPs, the descriptor also scales with the number of species:

For example, the angular symmetry functions for an ANN potential for a single species A describe the interactions of the central atom with two atoms of type A (A-A). For two species A and B, three interactions occur (A-A, A-B, and B-B), and for three species A, B, C, there are six (A-A, A-B, A-C, B-B, B-C, C-C). In general, for $N$ species the number of interactions is $N(N+1)/2$, i.e., the scaling is quadratic in the number of species. Further details about the symmetry function set up for multicomponent systems using the Behler-Parrinello approach along with actual parameters can also be found in reference Artrith et al., 2013.

When Ising-like pseudo spin variables are used to describe compositions, as for example in the cluster expansion (CE) method, an analogous scaling with the number of species occurs. The number of CE basis functions scales quadractically with the number of species when only pair clusters are considered. Generally, the scaling is on the order of the highest included n-body interaction, i.e., cubic for triplets, 4th order for quadruplet interactions, and so on. Mathematically, this relationship was worked out in reference Sanchez et al., 1984.

The multi-component implementation of the smooth overlap of atomic positions (SOAP) approach, is laid out in section 2.3 of reference De et al., 2016. The scaling is also quadratic, as it involves partial power spectra for each pair of species.

Appendix D Derivation of the expansion coefficients

The expansion coefficient $c^{(2)}_{\alpha}$ of the basis set expansion of the radial distribution function (RDF)

[TABLE]

where $\int\bar{\phi}_{\beta}(r)\phi_{\alpha}(r)\,\mathrm{d}r=\delta_{\beta\alpha}$ and $\{\bar{\phi}_{\alpha}\}$ is the orthogonal dual basis to $\{\phi_{\alpha}\}$, is given by

[TABLE]

Note that the RDF as defined in Eq. (4) of the main manuscript

[TABLE]

is only different from [math] for $0\leq{}r\leq{}R_{\textup{c}}$, so that the integral in Eq. (8) can be replaced by an integral over the entire space. Inserting the expression of the RDF Eq. (9) into the Eq. (8) yields

[TABLE]

which is the expression given in Eq. (6) of the main manuscript.

The derivation of the coefficients $c^{(3)}_{\alpha}$ of the angular expansion is completely analogous.

Appendix E Chebyshev polynomials of the first kind

The Chebyshev polynomials $\{T_{n}\}$ are defined by the recurrence relation

[TABLE]

The polynomials are orthogonal on the interval $[-1,1]$ with respect to a weight

[TABLE]

so that we choose the basis functions $\{\phi_{\alpha}\}$ and their duals $\{\bar{\phi}_{\alpha}\}$ on the interval $r\in{}(0,R_{\textup{c}})$ (for the radial expansion) as

[TABLE]

where $k=1/2$ for $\alpha=0$ and $k=1$ otherwise. For the angular expansion, the appropriate interval is $0\leq{}\theta<\pi$, so that $R_{\textup{c}}$ has to be replaced by $\pi$ in Eqs. (13) and (14).