Atomic energy mapping of neural network potential

TL;DR

This paper investigates how neural network potentials infer atomic energies from total energies, revealing vulnerabilities and proposing methods to improve atomic energy mapping through careful training set selection and invariant point monitoring.

Contribution

It introduces a framework for analyzing atomic energy mapping in neural network potentials and demonstrates how to enhance accuracy by strategic training practices.

Findings

Neural network potentials can have incorrect atomic energy mappings despite accurate total energy training.

Invariant points in feature space are useful for assessing atomic energy mapping.

Careful training set selection and monitoring improve atomic energy inference.

Abstract

We show that the intelligence of the machine-learning potential arises from its ability to infer the reference atomic-energy function from a given set of total energies. By utilizing invariant points in the feature space at which the atomic energy has a fixed reference value, we examine the atomic energy mapping of neural network potentials. Through a series of examples on Si, we demonstrate that the neural network potential is vulnerable to 'ad hoc' mapping in which the total energy appears to be trained accurately while the atomic energy mapping is incorrect in spite of its capability. We show that the energy mapping can be improved by choosing the training set carefully and monitoring the atomic energy at the invariant points during the training procedure.

Equations

$$E_{\text{tot}}^{\text{DFT}}=\sum _i{E}_{\text{at}}({\mathcal{R}}_i;R_{\text{c}}),$$

$$\begin{array}{rl}& E_{\text{tot}}^{\text{DFT}}=E_{\text{kin}}+E_{\text{XC}}+E_{\text{Coul}}\\ & =-\frac{1}{2}\int {\nabla }_{\mathbf{r}}^2\rho (\mathbf{r},{\mathbf{r}}^{\prime }){|}_{\mathbf{r}={\mathbf{r}}^{\prime }}d{\mathbf{r}}^{\prime }+\int \rho (\mathbf{r}){\epsilon }_{\text{XC}}(\rho (\mathbf{r}),\nabla \rho (\mathbf{r}))d\mathbf{r}\\ & +\frac{1}{2}\int \frac{\rho (\mathbf{r})\rho ({\mathbf{r}}^{\prime })}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d\mathbf{r}d{\mathbf{r}}^{\prime }-\sum _i\int \frac{q_i\rho (\mathbf{r})}{|\mathbf{r}-{\mathbf{r}}_i|}d\mathbf{r}+\sum _{i>j}\frac{q_i{q}_j}{|{\mathbf{r}}_i-{\mathbf{r}}_j|},\end{array}$$

$$\begin{gathered} E_{\text{Coul},i}=\frac{1}{2}\sum _{j\ne i}\int \frac{{\rho }_{\text{tot},i}(\mathbf{r}){\rho }_{\text{tot},j}({\mathbf{r}}^{\prime })}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d\mathbf{r}d{\mathbf{r}}^{\prime } \\ +\frac{1}{2}\int \frac{{\rho }_i(\mathbf{r}){\rho }_i({\mathbf{r}}^{\prime })}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d\mathbf{r}d{\mathbf{r}}^{\prime }-\int \frac{q_i{\rho }_i(\mathbf{r})}{|\mathbf{r}-{\mathbf{r}}_i|}d\mathbf{r}. \end{gathered}$$

$${\rho }_{\text{at},i}(\mathbf{r},{\mathbf{r}}^{\prime })={\rho }_{ii}(\mathbf{r},{\mathbf{r}}^{\prime })+\frac{1}{2}\sum _{j\ne i}^{|{\mathbf{r}}_j-{\mathbf{r}}_i|<R_{\text{c}}^3}{\rho }_{ij}(\mathbf{r},{\mathbf{r}}^{\prime }).$$

$$E_{\text{kin},i}=-\frac{1}{2}\int {\nabla }_{\mathbf{r}}^2{\rho }_{\text{at},i}(\mathbf{r},{\mathbf{r}}^{\prime }){|}_{\mathbf{r}={\mathbf{r}}^{\prime }}d{\mathbf{r}}^{\prime }.$$

$$E_{\text{at},i}=E_{\text{kin},i}+E_{\text{XC},i}+E_{\text{Coul},i},$$

$$E_{\text{at},i}⟶E_{\text{at}}^{\text{DFT}}(\mathcal{R};R_{\text{c}}),$$

$$E_{\text{tot}}^{\text{DFT}}=\sum _i{E}_{\text{at}}^{\text{DFT}}({\mathcal{R}}_i;R_{\text{c}})\simeq \sum _i{E}_{\text{at}}^{\text{NN}}({\mathbf{G}}_i;r_{\text{c}}).$$

Full text

Recently, machine-learning (ML) approaches to developing interatomic potentials are attracting considerable attention because it is poised to overcome the major shortcoming inherent to the classical potential and first-principles method, i.e., difficulty in potential development and huge computational cost, respectively. Favored ML models are the neural network Behler and Parrinello (2007); Zhang et al. (2018) and Gaussian process Bartók et al. (2010). In particular, the high-dimensional neural network potential (NNP) suggested by Behler and ParrinelloBehler and Parrinello (2007) is attracting wide interests with applications demonstrated over various materials encompassing metals,Artrith and Kolpak (2014); Boes and Kitchin (2017); Eshet et al. (2012) insulators,Artrith and Urban (2016); Li et al. (2017) semiconductors,Sosso et al. (2013); Behler et al. (2008) and molecular clusters.Kolb et al. (2016)

While the methodological advances are under rapid progress Gastegger et al. (2018); Li and Ando ; Imbalzano et al. (2018); Hajinazar et al. (2017); Onat et al. (2018); Pham et al. (2016); Jeong et al. (2018), the conceptual foundation of NNP is still elusive, partly due to the black-box nature of the neural network. Furthermore, NNP infers atomic energies while it is trained over total energies that are sums of atomic energies. This obscures the nature of training procedure and makes it difficult to assess learning quality. Responding to this, in this Letter, we try to address basic questions on NNP such as ‘what is the intelligence of NNP?’ and ‘how the learning quality is determined?’. At variance with general views, we show that the core of training procedure in NNP is to infer the reference atomic energy grounded on the density functional theory (DFT), from the given relationship between the structure and total energy. With examples on Si, we demonstrate that NNP is prone to ad hoc mapping in which the total energy is trained accurately but the atomic energy mapping is incorrect. We also show that the reference atomic energies can serve as a tool to assess the learning quality of NNP.

Most ML potentials are based on the representability of the DFT total energy ($E_{\text{tot}}^{\text{DFT}}$) as a sum of the atomic energy ($E_{\text{at}}$) that depends on the local environment within a certain cutoff radius ($R_{\text{c}}$):

[TABLE]

where $i$ is the atom index and $\mathcal{R}_{i}$ is the collection of relative position vectors of atoms lying within $R_{\text{c}}$ from the $i$th atom. (For simplicity, we assume a unary system that is large enough that various cutoff spheres in the following discussions do not self-overlap under periodic boundary conditions and wave functions are effectively real-valued.) As is well known, the total energy can be expressed by integration of the local energy density, although it is not unique.[][; AppendixHandreferencestherein.]24-MartinBook Then, by partitioning the space into non-overlapping atomic volumes, one can assign energies to each atom whose sum equals to the total energy.Yu et al. (2011); Popelier (2015) In consideration of locality or ‘nearsightedness’ of the electronic structure,Kohn (1996) which empowers the $\mathcal{O}(N)$ approach Goedecker (1999), it would be formally viable to define $E_{\text{at}}(\mathcal{R}_{i};R_{\text{c}})$ in Eq. (1) within DFT, which depends only on local features and so is transferable. In the below, we elaborate on this explicitly, with a particular attention to the transferable range.

Within the semilocal density approximation, $E_{\text{tot}}^{\text{DFT}}$ can be expressed in terms of the one-electron density matrix $\rho(\mathbf{r},\mathbf{r^{\prime}})$ and the electron density $\rho(\mathbf{r})=\rho(\mathbf{r},\mathbf{r})$:

[TABLE]

where the atomic unit is used, $\varepsilon_{\text{XC}}$ is the exchange-correlation energy density, $q_{i}$ and $\mathbf{r}_{i}$ are the ionic charge and position of the $i$th atom, respectively. Under the assumption that $\mathcal{O}(N)$ methods, in particular the divide-and-conquer (DAC) approachYang (1991); Yang and Lee (1995), work well for given systems, we will explicitly show that i) each energy term can be split without any loss into atomic contributions that are defined locally around each atomic site, and ii) the atomic energy depends only on nearby atoms such that it is transferable to other systems as long as local environments are maintained.

We start with partitioning the space into atomic cells without gaps or overlapping (for instance, Voronoi cells). Let $V_{i}$ be the cell enclosing the $i$th atom. We define $\rho_{i}(\mathbf{r})$ as $\rho_{i}(\mathbf{r})=\rho(\mathbf{r})[\mathbf{r}\in V_{i}]\,$ where […] is the Iverson bracket whose value is 1 (0) when the logical proposition in the bracket is true (false). Obviously, $\rho(\mathbf{r})=\sum_{i}\rho_{i}(\mathbf{r})$. It is easily seen that $E_{\text{XC}}$ is the sum of the atomic exchange-correlation energy ($E_{\text{XC},i}$) that is obtained by substituting $\rho_{i}(\mathbf{r})$ for $\rho(\mathbf{r})$ in the integrand of $E_{\text{XC}}$. As is assumed in the DAC method,Yang (1991) the charge density at a certain point is influenced by only nearby atoms if the local chemical potential of electrons is fixed. This means that $\rho_{i}(\mathbf{r})$, and hence $E_{\text{XC},i}$ is affected by atomic arrangements within a certain cutoff ($R_{\text{c}}^{1}$) from $\mathbf{r}_{i}$.

Next, we define the total charge density in $V_{i}$: $\rho_{\text{tot},i}(\mathbf{r})=q_{i}\delta(\mathbf{r}-\mathbf{r}_{i})-\rho_{i}(\mathbf{r})$. It is straightforward to show that $E_{\text{Coul}}$ can be expressed as a summation of the atomic Coulomb energy, $E_{\text{Coul},i}$, defined as follows:

[TABLE]

The first term on the right-hand side of Eq. (3) is long-ranged, which is incompatible with the finite cutoff. However, electrostatic interactions are effectively screened or cancelled in many condensed phases, so it would be a reasonable approximation to ignore them beyond a certain cutoff ($R_{\text{c}}^{2}$). Thus, we omit the Coulomb interaction between $\rho_{\text{tot},i}$ and $\rho_{\text{tot},j}$ if $\lvert\mathbf{r}_{j}-\mathbf{r}_{i}\rvert>R_{\text{c}}^{2}$. Since $\rho_{i}(\mathbf{r})$ and $\rho_{\text{tot},i}(\mathbf{r})$ are influenced by atoms within $R_{\text{c}}^{1}$ (see above), $E_{\text{Coul},i}$ depends on atoms inside $R_{\text{c}}^{1}+R_{\text{c}}^{2}$ (neglecting the volume of $V_{i}$). To note, some implementations of NNP explicitly describe the long-range Coulomb potential, separately from short-ranged atomic energies.Artrith et al. (2011); Morawietz et al. (2012)

As the last step, we discuss the locality of $E_{\text{kin}}$. Since $\rho(\mathbf{r},\mathbf{r}^{\prime})$ decays exponentially with $\lvert\mathbf{r}-\mathbf{r}^{\prime}\rvert$ in insulators and metals at finite temperatures Goedecker (1999), one can neglect $\rho(\mathbf{r},\mathbf{r}^{\prime})$ when $\lvert\mathbf{r}-\mathbf{r}^{\prime}\rvert$ is bigger than a cutoff ($R_{\text{c}}^{3}$), which is utilized in the density-matrix-based DAC method Yang and Lee (1995). Therefore, for a given position $\mathbf{r}$, $\rho(\mathbf{r},\mathbf{r}^{\prime})$ is determined by the atomic configurations within a cutoff distance ($R_{\text{c}}^{4}$) from $\mathbf{r}$, which should be larger than $R_{\text{c}}^{3}$. With the projected density matrix $\rho_{ij}(\mathbf{r},\mathbf{r}^{\prime})=\rho(\mathbf{r},\mathbf{r}^{\prime})[\mathbf{r}\in V_{i}][\mathbf{r}^{\prime}\in V_{j}]$, we define the atomic density matrix $\rho_{\text{at},i}(\mathbf{r},\mathbf{r}^{\prime})$ as follows:

[TABLE]

It can be shown that $\rho(\mathbf{r},\mathbf{r}^{\prime})=\sum_{i}\rho_{\text{at},i}(\mathbf{r},\mathbf{r}^{\prime})$ and $\rho_{\text{at},i}(\mathbf{r},\mathbf{r}^{\prime})$ depends only on the atomic arrangements within $R_{\text{c}}^{4}$ from the $i$th atom (neglecting the volume of $V_{i}$). The atomic kinetic energy is then given in the following:

[TABLE]

Since the kinetic-energy operator is linear, the sum of the atomic kinetic energy is equivalent to the total kinetic energy.

Combining the above analyses, the atomic energy of the $i$th atom formally derives from the DFT calculations:

[TABLE]

and $E_{\text{tot}}^{\text{DFT}}=\sum_{i}E_{\text{at},i}$. By evaluating $E_{\text{at},i}$ in various structures, one can obtain in principle the atomic energy as a continuous function of the local environment:

[TABLE]

where $R_{\text{c}}=\max(R_{\text{c}}^{1}+R_{\text{c}}^{2},R_{\text{c}}^{4})$. Note that the atomic energy is not unique because it depends on the way to define atomic cells.

The existence of $E_{\text{at}}^{\text{DFT}}$ implies that the objective of the present machine learning is to identify $E_{\text{at}}^{\text{DFT}}$ when only total energies are informed. This is at variance with the conventional view that NNP is merely an interpolation of given total energies. Behler ; Khorshidi and Peterson (2016) Mathematically, the neural network has the capability to infer the underlying function when only sums of function values are provided. (See examples in Supplemental Material 111See Supplemental Material below for the capability of neural network learning from the sum, atomic energy mapping of Ni nanocluster, the training set and computational details, the estimation of prediction uncertainty, the analysis on the distance to training set, and the metric to quantify connectivity in $\mathbf{G}$ space. on a piecewise cubic spline and the embedded atom potential.) To reduce the huge dimension of $\mathcal{R}$ and obtain $E_{\text{at}}$ in a computationally feasible way, two approximations are adopted. First, the cutoff radius is reduced from $R_{\text{c}}$, which should be fairly large for high accuracy, to $r_{\text{c}}$ that is usually chosen to be 6-7 Å. This is a reasonable range because the chemical influence rapidly diminishes beyond this boundary. Second, the local environment is described by feature vectors whose dimension is significantly lower than for $\mathcal{R}$. The popular choices are smooth-overlap-of-atomic-positions (SOAP) Bartók et al. (2013) or symmetry function vectors ($\mathbf{G}$) Behler (2011). These feature vectors also automatically incorporates rotational and translational invariance inherent to the atomic energy. Here, we employ the symmetry function. Thus,

[TABLE]

The accuracy of NNP therefore hinges on how close $E_{\text{at}}^{\text{NN}}$ obtained through machine learning is to the reference $E_{\text{at}}^{\text{DFT}}$ over the configurational space spanned by the given training set. However, since $E_{\text{at}}^{\text{NN}}$ is fitted to the total energies, rather than directly to $E_{\text{at}}^{\text{DFT}}$, the ML procedure does not necessarily guarantee sufficient accuracies in $E_{\text{at}}^{\text{NN}}$. That is to say, $E_{\text{at}}^{\text{NN}}$ can reproduce total energies in the training set precisely but deviate significantly from $E_{\text{at}}^{\text{DFT}}$. Indeed, we will demonstrate that NNP is vulnerable to such ‘ad hoc’ energy mapping, which leads to incorrect total energies in related configurations and undermines the transferability of NNP. To note, the ad hoc mapping should be distinguished from the gauge-dependent degree of freedom in the energy density.Martin (2004)

Even though the existence of $E_{\text{at}}^{\text{DFT}}$ was shown formally in the above, the actual calculation of $E_{\text{at}}^{\text{DFT}}$ would be highly costive. (We note a recent effort to directly train NNP over atomic energies. Huang et al. (2019)) Furthermore, there exist an infinite number of valid $E_{\text{at}}^{\text{DFT}}$, making it hard to grade the energy mapping of $E_{\text{at}}^{\text{NN}}$. However, there are invariant points in the $\mathbf{G}$ space at which $E_{\text{at}}^{\text{DFT}}$ is uniquely defined without any degree of freedom. For instance, in the crystalline Si, all the atoms are equivalent, and so the total energy per atom is simply equal to $E_{\text{at}}^{\text{DFT}}$ for the corresponding $\mathbf{G}$. Transforming lattice vectors of the unit cell also results in similar conditions. In the below, utilizing these special $\mathbf{G}$ points, we will analyze the atomic energy mapping for three examples on Si that are progressively more complicated.

As the first example, we train NNPs for crystalline Si by adopting an in-house code named SIMPLE-NNLee et al. . The training (validation) set consists of 350 (150) MD snapshots of the 64-atom cubic supercell under the NVT condition of 1000 K and the equilibrium volume at 0 K. After training, the root-mean-squared error (RMSE) in the total energy and atomic force is 0.9 (0.9) meV/atom and 0.10 (0.11) eV/Å for the training (validation) set, respectively. (See the Supplemental Material Note (1) for further details in DFT calculations and NNP training.) Atomic vibrations during MD give rise to local expansion or compression. As a result, atomic configurations around certain Si atoms resemble those in the crystalline phase under hydrostatic pressures, which forms the equation of state (EOS) and corresponds to invariant $\mathbf{G}$ points explained above. This is confirmed by principal-component analysis (PCA) and measuring the shortest distances from the invariant $\mathbf{G}$ points to the training set. (See Supplemental Material Note (1).) This implies that atomic energies at the $\mathbf{G}$ points along EOS are learnable although they do not belong to the training set. Therefore, if atomic energies are properly mapped, NNP should be able to predict correctly the energy-volume relation at 0 K.

Figure 1 compares EOS inferred by the as-trained NNP (blue line) with DFT results (black dots). The light shade means prediction uncertainty evaluated by ensembles of NNP Peterson et al. (2017) (see also Supplemental Material Note (1)). The squared bracket indicates the range of the volume whose $\mathbf{G}$ is in close proximity to the training set. Interestingly, NNP predicts correctly the energy at equilibrium but energies at other volumes significantly deviate from the DFT curve with errors far bigger than RMSE in the total energy. That is to say, NNP predicts the total energy correctly but the atomic energy is markedly wrong, which corresponds to the ad hoc energy mapping. From the continuity in $E_{\text{at}}$, the incorrect energy mapping should affect other training points neighboring invariant $\mathbf{G}$ points, implying that the ad hoc mapping extends over a significant portion of the training set.

The ad hoc mapping in the above example happens because the training set consists of structures with a fixed volume. This condition constrains the local expansion and contraction to occur concurrently within the same structure. Consequently, any additional atomic energy that varies linearly with the volume does not affect the total energy, and so the slope of EOS at the equilibrium volume becomes an arbitrary number. The ad hoc mapping in this case can be resolved by considering structures with different volumes or including virial stress in the loss function. For instance, the red line in Fig. 1 shows EOS predicted with NNPs that are trained with MD snapshots from NPT ensembles at 1000 K and zero pressure. During MD, the supercell expands or shrinks, avoiding the exact cancellation among the local volume changes. As a result, it is seen that NNPs can predict the slope and curvature of EOS reasonably.

The second example concerns a surface model of Si. The training set consists of MD trajectories of Si(100)-(2$\times$2) symmetric slab in the NVT condition at a certain temperature between 100 and 1000 K [see Fig. 2(a)]. To assess the learning quality, we compare atomic energies for the geometry relaxed at 0 K with DFT. Unlike crystalline Si in the previous example, the reference $E_{\text{at}}^{\text{DFT}}$ is not available directly. Nevertheless, Si atoms inside the slab (blue atoms) have neighborhood similar to that in the crystal (see a dashed circle). Therefore, $E_{\text{at}}$ in this region should be close to the crystalline $E_{\text{at}}^{\text{DFT}}$ at the equilibrium volume [$E_{\text{at}}^{\text{DFT}}(\text{bulk})$]. Since $E_{\text{tot}}^{\text{DFT}}$ is available for the whole structure, the average $E_{\text{at}}^{\text{DFT}}$ for the surface region (red atoms) can be obtained as $[E_{\text{tot}}^{\text{DFT}}-N_{\text{b}}\cdot E_{\text{at}}^{\text{DFT}}(\text{bulk})]/N_{\text{s}}$, where $N_{\text{b}}$ and $N_{\text{s}}$ are the number of atoms in the bulk and surface regions, respectively. By taking the difference in averaged values of $E_{\text{at}}^{\text{NN}}$ and ${E}_{\text{at}}^{\text{DFT}}$ in each region, one can quantify average mapping errors, $\Delta\bar{E}_{\text{at}}(\text{bulk})$ and $\Delta\bar{E}_{\text{at}}(\text{surface})$, respectively.

Figure 2(b) presents $\Delta\bar{E}_{\text{at}}(\text{bulk})$ and $\Delta\bar{E}_{\text{at}}(\text{surface})$ for NNPs trained over MD trajectories at different temperatures. At a low temperature of 100 K, the mapping error is $-108$ and 76 meV/atom for bulk and surface regions, respectively, which is far bigger than RMSE (0.3 meV/atom). This is another example of ad hoc energy mapping; NNP correctly predicts the total energy because errors in the atomic energy mapping cancel with each other. In Fig. 2(b), it is intriguing that the mapping error gradually decreases as the temperature in the training set increases, and at the high temperature of 1000 K, the magnitude of mapping errors becomes comparable to RMSE in the total energy (3 meV/atom).

To understand the temperature-dependent mapping error, we examine in Fig. 2(c) the distribution of training points in the $\mathbf{G}$ space using PCA on the training sets at 100 and 1000 K. It is seen that at 100 K, the training points corresponding to the bulk and surface region are well separated. In contrast, energetic vibrations at 1000 K result in much broader distribution of training points such that bulk and surface regions are slightly connected. (Other combinations of principal axes show similar behaviors.) As schematically drawn in Fig. 2(d), if clusters of training points are separate as in 100 K, the machine learning is prone to ad hoc mapping because any cancelling offsets give almost the same total energy and atomic forces. On the other hand, at higher temperatures with every region connected to some degrees, $E_{\text{at}}$ at intermediate configurations helps adjust the energy offset between basins. In Supplemental Material Note (1), we define a metric that can quantify the connectivity in the $\mathbf{G}$ space and show that the mapping error is sufficiently small when the connectivity is higher than a cutoff.

Albeit simple, the above cases substantiate the ad hoc mapping that originates from limitations in the training set. In practice, a single training set usually encompasses diverse structures such as bulk, surfaces, and defects, and chances are that the ad hoc mapping can be avoided in principle. Nevertheless, the error-cancelling energy offsets as in Fig. 2(d) are omnipresent, which can go unnoticed if the training procedure is monitored by RMSE only. To show this, we generate a training set from MD simulations of a 239-atom Si nanocluster with Wulff-constructed {100}, {110}, and {111} facets at 1000-1700 K. [See Fig. 3(a) for the structure relaxed at 0 K.] The analysis on the connectivity (see above) confirms that training points are well connected.

In Fig. 3(b), we plot RMSE for the total energy and force with respect to the training epoch. It also shows $\Delta\bar{E}_{\text{at}}(\text{bulk})$ and $\Delta\bar{E}_{\text{at}}(\text{surface})$ for the (100)-(2$\times$2) surface model in Fig. 2(a). The analysis similar to Fig. S3 shows that the $\mathbf{G}$ points in the (100)-(2$\times$2) slab model are in the vicinity of training points, and hence they are learnable. Therefore, NNP is expected to predict surface and bulk energies in reasonable agreement with DFT results. In Fig. 3(b), it is seen that RMSE remains almost constant after about 100 epochs while $\Delta\bar{E}_{\text{at}}(\text{bulk})$ and $\Delta\bar{E}_{\text{at}}(\text{surface})$ converge at much slower rates. This indicates a risk in concluding the training convergence in terms of RMSE, and supports $E_{\text{at}}$ at invariant $\mathbf{G}$ points as alternative convergence parameters. Obviously, if the crystalline structures are included in the training set, $\Delta\bar{E}_{\text{at}}(\text{bulk})$ would converge as fast as RMSE, but this does not guarantee the proper energy mapping at other training points. Therefore, we suggest to collect invariant $\mathbf{G}$ points as a separate test set for monitoring the atomic energy mapping, rather than including them in the training set, at least in the initial stage of training.

After a sufficient number of epochs, the surface energies for (100)-(2$\times$2), (110)-(2$\times$1), and (111)-(2$\times$1) slab models that are fully relaxed by NNP agree with DFT results within 8%. (The corresponding errors by NNP trained up to 200 epochs are within 20%.) It is intriguing that just one type of structure (nanocluster) can train NNP over such a wide range of configurations when the energy mapping is correct. This implies that NNPs with proper mapping are more transferable than those with ad hoc mapping, which may contribute to improving the stability of MD simulations.Jeong et al. (2018) It will be also useful in developing general-purpose NNPs.Bartók et al. (2018) Finally, we find that monitoring the energy mapping is helpful in selecting training parameters such as the regularization parameter of the neural network.

In conclusion, we showed that the aim of training NNP is to learn the atomic energy function defined at the DFT level from total energies, and the transferability of NNP lies in the accuracy of atomic energy mapping. The invariant $\mathbf{G}$ points with the unique ${E}_{\text{at}}^{\text{DFT}}$ provided ways to examine the atomic energy mapping. Several examples confirmed that NNP is vulnerable to ad hoc mapping due to limitations in the training set and/or certain choices of computational parameters. The energy mapping can be improved by choosing the training set carefully and monitoring the atomic energy at the invariant points during the training procedure. By clarifying what NNP actually learns, the present work will contribute to constructing accurate and transferable machine-learning potentials.