In operando active learning of interatomic interaction during large-scale simulations

We now turn to the question of how to construct an appropriate training set from configurations ${\{{r}_i^\ast\}}$ in which the extrapolation grade exceeds a certain threshold. One of the currently most accurate and efficient ways to compute quantum mechanical quantities, such as energies, forces and stresses, appears to be plane-wave, i.e. periodic, DFT (cf e.g. [44]). However, if we simply apply periodic boundary conditions to ${\{{r}_i^\ast\}}$, the neighborhoods of atoms at the periodic domain boundaries would differ from those in the large-scale problem. This is in particular severe for 'higher-than-zero-dimensional defects', e.g. line (dislocations), surface (grain boundaries) or volume (cracks, voids) defects. The problem is thus similar to hybrid quantum/classical mechanics (QM/MM) methods, where a small, but arbitrarily-shaped, cluster of atoms, treated with an ab initio model, is embedded in a much larger domain computed with interatomic potentials. To consider such clusters of atoms using plane-wave DFT, Woodward and Rao [5] developed two methods in the context of screw dislocations:

• (i)  
  In the first method, a vacuum region was introduced to separate the QM cluster from its periodic images. To prevent spurious free surface effects on cluster atoms, an additional buffer region replicating their true neighborhood was used in-between.
• (ii)  
  In the second method, the buffer region from method (i) was extended to the periodic domain boundaries to avoid the vacuum region.

Method (ii) allows for a smaller cell size in comparison with method (i)—at the expense of a larger buffer region. However, neither method seems to be preferable over the other in terms of accuracy and computational cost since both have been equally successfully employed in many QM/MM studies (e.g. [1, 28, 45, 46]). In principle, they therefore also apply to our algorithm, however, none of them allows one to compute the quantum mechanical energy (of the cluster) which is a global quantity corresponding to the entire periodic cell ² .

Fortunately, in the case of MLIPs, we do not stringently require exact total energies—it suffices that the MLIP is able to interpolate on the large-scale configurational space. Nevertheless, it is still necessary to prevent free surfaces and/or other types of spurious boundary effects to avoid training on artificial neighborhoods irrelevant to the large-scale problem. We therefore present an improved version of method (ii) which symmetrizes the atomic displacements at the periodic domain boundaries to ensure that the neighborhood of the atoms outside the dislocation core is close to the one they would 'feel' in the large-scale problem.

The idea is to extract the displacements $\boldsymbol{\tilde{u}}_i(\boldsymbol{r}_{0,i}) = \boldsymbol{r}_i - \boldsymbol{r}_{0,i}$, where $\boldsymbol{r}_{0,i}$ is the ideal lattice site of the ith atom, around the dislocation core and apply it to a periodic configuration, small enough to be computable by DFT but without introducing atomic neighborhoods different from those in the large-scale configuration. To that end, we proceed as illustrated in figure 2 (a)–(1) and define the rectangular region of size $l_{\rm x} \times l_{\rm y}$, containing the dislocation core atoms ${\{{r}_i^\ast\}}$ which we check on extrapolation. Next, we define three replicas of this rectangular region and translate them according to figure 2 (a)–(2) by $l_{\rm x}$ with respect to the x-axis and/or $l_{\rm y}$ with respect to the y-axis.

Assuming that the dislocation line is parallel to the z-axis, we then define the mirrored displacement field

$$\begin{equation} \boldsymbol{u}(\boldsymbol{x}) = \begin{cases} \boldsymbol{\tilde{u}}(x,y,z) &\boldsymbol{u}\in\{{r}_{0,i}^\ast \}, \\ \boldsymbol{\tilde{u}}(-x-l_{\textrm{x}},y,z) &\boldsymbol{u}\in\{{r}_{0,i}^1 \}, \\ \boldsymbol{\tilde{u}}(x,-y-l_{\textrm{y}},z) &\boldsymbol{u}\in\{{r}_{0,i}^2 \}, \\ \boldsymbol{\tilde{u}}(-x-l_{\textrm{x}},-y-l_{\textrm{y}},z) &\boldsymbol{u}\in\{{r}_{0,i}^3 \}, \end{cases} \end{equation} \tag{ 5 }$$

where $\{{r}_{0,i}^1 \}$–$\{{r}_{0,i}^3 \}$ are the ideal lattice sites in the translated regions. Imposing this displacement gives the configurations $\{{r}_i^1 \}$–$\{{r}_i^3 \}$ shown in figure 2 (a)-(2). The union of ${\{{r}_i^\ast\}}$ and $\{{r}_i^1 \}$–$\{{r}_i^3 \}$ is clearly periodic and from the figure it can be seen that this procedure creates minimal spurious effects in the sense that the atomic neighborhoods are not 'too far' from the large-scale problem, indicated by the centrosymmetry parameter which is less than 0.04 in the vicinity of the periodic domain boundaries—which is more than two orders of magnitude smaller than for the core atoms. This configuration appears similar to the well-known quadrupolar cell (cf e.g. [47, 48]), yet it is emphasized that our configuration fully includes the elastic far-field contribution due to external boundary conditions intrinsic to the large-scale problem.

We may further reduce this cell to the dipole-like configuration framed by the continuous line in figure 2 (a)-(3) along similar lines as for the quadrupolar cells if the displacement is point symmetric with respect to the dislocation positions. To see why this choice is indeed justified we refer to our analysis in Appendix B.

It should be noted that this construction applies independently of the position of the dislocation core as shown in figure 2 (b). That is, if the dislocation starts moving, e.g. due to some far-field applied stress, we detect the dislocation core position and proceed in the same way as above by shifting the rectangular region accordingly.

Even though the proposed symmetrization technique ensures that no artificial neighborhoods occur in the training configuration, we would like to point out that there is potential freedom in modifying the positions of the boundary atoms: if, e.g. the symmetry of the displacement field breaks due to additional long-range effects, one may introduce an intermediate step which minimizes the centrosymmetry parameter, or, alternatively, the extrapolation grade, of atoms residing in some boundary layer.

Important remark. We find it important to emphasize that, despite the fact that the dislocation lines are extremely (some may say 'unphysically') close to each other in the training set, the fitted model is still expected to (and does in fact, as we will see in next section) give accurate predictions. The reason for such training sets being representative of configurations with normal, 'physical' distance between dislocations rests on the main assumption behind the success of machine-learning potentials: that the perturbations to electronic degrees of freedom decrease fast and therefore the electronic interaction can be well-described by the relative positions of atoms in a small neighborhood (like 5 Å). Hence, if the atomic neighborhoods in the training set are representative of those in a large-scale simulation, then the potential is expected to give accurate predictions. The difference in the long-range elastic energy between training and actual configurations will be resolved by the atomistic deformation which is explicitly represented in the problem of interest.

In particular, total energies and forces of the small periodic training configurations are not used to compute any physical quantities (defect formation energies, energy barriers etc). On the contrary, they are computed in the large-scale simulation in which we are free to choose any boundary condition (e.g. flexible boundary conditions [5, 49–51]). Hence, there is no need for the application of correction methods (e.g. [52, 53]) employed when operating solely on small, fully periodic DFT cells due to the interaction of defects with their periodic images through elastic fields . This is another advantage of our method.

For the dislocation core relaxation we set the bounds $\gamma_{\rm {min}}$ and $\gamma_{\rm \textrm{max}}$, within which we consider configurations as potential candidates for the training set, to 1.1 and 10, respectively. This choice is close to the optimum (cf [34]) and should therefore yield the best possible accuracy that one could achieve with an MTP in practice. Using this parameter set, the active learning algorithm converged after 63 iterations with a training set size of 40 configurations. The training errors reported in table 1 are in the range of the numerical approximation of DFT codes themselves. Therefore, the MTP essentially reproduces the training data without introducing further errors.

The differential displacement map for the relaxed core structure computed with the MTP is shown in figure 3. It can immediately be seen that our algorithm predicts the compact core structure from previous established DFT studies (cf e.g. [48, 57]). To quantify this result, we show the energy difference between the initial linear elastic and the final relaxed configuration in table 2. This is of course not feasible to compute with a reference DFT calculation so that we compare our result to the 135-atom periodic configuration. The energy difference is ≈ 2× larger for the periodic problem, which is expected since there are two dislocations in this system. In addition, we have computed the energy difference for the screw-cylinder configuration using the EAM4 potential as a reference model—where we can compute the exact solution. Here, the error of the MTP-EAM4 is ≈ 2% which confirms the excellent agreement between both models and provides further evidence that our algorithm reliably extracts the core energy; note that an error of 1–2 % also occurs for the periodic problem (compare the EAM4 and MTP-EAM4 values in table 2).

Zoom In Zoom Out Reset image size

It is interesting to note that it was 'easier' for the MTP to reproduce DFT than EAM4—the difference between MTP and DFT is less than 1 meV per Burgers vector, while the difference between MTP and EAM4 is of the order of 10 meV per Burgers vector. Speculatively, we attribute it to the fact that the MTP has a flexible functional form which makes it easy to reproduce the smooth underlying DFT energy (at least, when the configurational space is not too large), whereas fitting an EAM potential introduces an oscillatory behavior commensurate with atomic shell radii in the pair potential function as shown in figure 4. Such a nonsmooth behavior appears to be more difficult to fit with a smooth MTP radial basis than fitting the DFT energy.

Zoom In Zoom Out Reset image size

Using the relaxed configuration from the previous section, we now compute the Peierls barrier for screw dislocation motion under zero stress by calculating the minimum energy path using the nudged elastic band [62, NEB] method. We create the final image by translating the displacement field of the relaxed core structure by $\sqrt{2/3}\,\boldsymbol{a_0}$ along the $\langle 112 \rangle$ direction, i.e. to the next easy-core position on the {110} plane. To create the initial guess for the intermediate images, we use the same displacements according to their reaction coordinate.

For this problem we have chosen the extrapolation bounds $\gamma_{\rm {min}} = 2$ and $\gamma_{\rm \textrm{max}} = 10$. Using an NEB with 5 images, the active learning algorithm converged after 38 iterations with a training set size of 22 configurations. Therefore, our algorithm only required a full DFT calculation on the 135-atom symmetrized core configuration (cf figure 2 (4)) roughly every 9th iteration. The training errors in table 3 are slightly worse than for the core relaxation, which is expected since the training set is now more diverse due to the intermediate core configurations. However, a maximum error in the total energy of 5.6 meV (per configuration) is still a very small quantity since the Peierls barrier is more than an order of magnitude higher (cf table 4).

The corresponding energy curve is shown in figure 5 (a). Again, since a reference DFT calculation is unfeasible, we have compared our result to the one obtained from a periodic calculation using solely DFT. According to the convergence study by Grigorev et al [63], the difference between a periodic configuration of 135 atoms and an effectively infinite configuration should be in the range of only a few percent. From table 4, it follows that the difference is ≈ 4% and thus in agreement with their analysis. For testing purposes, we have also carried out calculations with the EAM4 potential as our reference model, as in the previous section. Here, the energy difference between the periodic and the cylinder configuration is ≈ 6%. Moreover, the Peierls energies for the EAM4 and the MTP-EAM4 coincide almost perfectly which further validates our training methodology. In addition, the differential displacement map for the relaxed central image is shown in figure 5 (b) which demonstrates that the MTP-DFT also correctly predicts the same intermediate configuration observed in our and recent DFT studies.

Zoom In Zoom Out Reset image size

Finally, we test our algorithm for the screw-cylinder problem subject to a far-field applied stress in order to predict the Peierls stress, that is, the stress under which the dislocation starts moving. To this end, we apply a deformation corresponding to the shear stress $\boldsymbol{\tau}_{\rm yz}$ to all atoms and subsequently relax the system. During relaxation, we detect the position of the dislocation core and update the extrapolation region accordingly (cf figure 2 (b)). For this problem, we have used $\gamma_{\rm \textrm{min}} = 2$ and $\gamma_{\rm {max}} = 10$ as for the barrier calculations.

We have applied an initial stress of ≈ 1.74 GPa and increased it in 0.1 GPa-steps upon convergence. The dislocation then moved to the next easy-core position when the stress was 2.34 GPa. Hence, we define the obtained Peierls stress as 2.29± 0.05 GPa. This value is in the same range as previous DFT studies as shown in table 5. Nevertheless, we would like to remind the reader that these previous studies use indirect methods; e.g. the method in [64] is based on fitting a Kocks-type function using the Peierls barriers from several periodic calculations subject to different applied stresses as fitting parameters. They are therefore computationally significantly more expensive, require prior knowledge of the slip plane and are also prone to errors due to the analytic fit.

In order to judge on the accuracy of our method, we compute the Peierls stress with the EAM3 potential as the reference model which allows us to perform an exact large-scale calculation ⁴ . Following the same procedure as described above, except that we now increment the applied stress by 0.01 GPa, the predicted Peierls stress for the MTP-EAM3 agrees well with our large-scale reference calculation showing an error of only 6%–8% (cf table 5). Moreover, this corresponds to the relative fitting errors for the MTP-EAM3 shown in table D10. Given that the training errors for the MTP-EAM3 and the MTP-DFT (table 6) are in the same range, it is likely that the DFT Peierls stress can also be predicted with the MTP-DFT up to an error of 6–8%. Should we need a higher accuracy, we may simply increase the number of basis functions of the MTP but note as well that errors up to 10% are absolutely tolerable when the Peierls stress is used to calibrate a higher-scale model, such as discrete dislocation dynamics.

To validate the correct functioning and demonstrate the capabilities of the active learning algorithm, we show the evolution of the dislocation position (with respect to the origin) and the training set size as a function of the iteration in figure 6 (a). Having crossed the first Peierls barrier after ≈ 550 iterations, the training set contained 20 configurations, increased to 22 while crossing the second Peierls barrier, and then continued gliding without triggering new DFT calculations. This behavior is expected since the dislocation moves by alternating from ① easy-core ② intermediate-core ③ easy-core $\rightarrow$ ... etc as shown in figure 6 (b) and does therefore not encounter any significantly new, i.e. extrapolative, configurations.

Zoom In Zoom Out Reset image size

We have developed an active learning algorithm for large-scale atomistic simulations using MTPs, a class of MLIPs. Our algorithm allows for a novel way of performing large-scale simulations with DFT accuracy fully automatically by letting active learning find the right training data during the simulations—without the need for any prior (passive) training. We have applied the algorithm to model screw dislocation motion in bcc tungsten and shown that the MTP is able to reproduce known mechanically relevant properties from the literature, such as core structure, transition state barriers and Peierls stress, alongside with a significant reduction in the necessary amount of DFT calculations compared to standard methods.

More specifically, we have shown that active learning is able to reliably detect a minimum set of sufficiently distinct training configurations. This is in particular demonstrated in section 4.3, where the dislocation moves through the effectively infinite crystal: in the beginning, active learning requests new training data during the change of the core structure but stops when the dislocation repetitively encounters the same configurations after crossing the first Peierls valley. Moreover, using our symmetrization strategy to convert extrapolative regions into periodic configurations, the MTP gave accurate predictions for forces as well as energy differences.

Our algorithm can in principle be applied by a user in the same way as a standard large-scale atomistic simulation. Therefore, no additional expert knowledge is required in order to extract physical quantities, such as the Peierls stress, from indirect DFT methods (cf section 4.3). This opens new prospects for using MLIPs in general as a robust tool for DFT-accurate simulations of extended defects.

Our example in section 4.3 further shows that we can include far-field applied stresses which influence the core structure of the dislocation and its the elastic field. We therefore expect that our algorithm is directly applicable to related problems involving interactions between dislocations and different types of defects, e.g. vacancies or other dislocations.

While we have assumed that nonlinearities are confined to the small-sized extrapolation region(s), we have demonstrated that active learning significantly reduces the amount of actual DFT calculations. If the size of the defect core increases to several hundreds of atoms, e.g. in the case of dissociated dislocations in face-centered-cubic lattices, the saved computing capacities may thus be invested in a few larger calculations in order to compute energies; combined with many small cluster calculations in order to compute forces.

Further, we think that our methodology will also be useful for inherently three-dimensional problems, e.g. curved dislocations. For this class of problems, we currently envision two possible approaches. The first possibility is to train a MLIP on several selected straight dislocations with different character angles. The idea is then to switch off active learning and use the MLIP as a standard interatomic potential. The second possibility is to proceed as we do in the present work, i.e. we carve out a collection of atoms in the vicinity of the dislocation core and construct a periodic configuration. This procedure is illustrated in figure 7 for a three-dimensional dislocation in a bcc lattice moving by means of the kink-pair nucleation mechanism (see, e.g. [66]). Both directions will be explored in future work.

Zoom In Zoom Out Reset image size

It is finally worthwhile to emphasize that our methodology does not exclude multi-component materials involving alloying elements and chemical impurities, e.g. hydrogen or helium, as most MLIPs already include the necessary descriptors (e.g. [10, 37]). The discovery of new materials mainly comes through the search for novel alloys and we think that investigating corresponding mechanisms using MLIPs in combination with active learning will be one of the most auspicious potential applications since the configurational space will likely be too vast to be explored with only passively selected training sets.

The basis functions from equation (2) for MTPs are given by the following form [7]

$$\begin{equation} B_\alpha(\{{r}_{ij}\}) = \prod_{l = 1}^k M_{\mu_l(\alpha),\nu_l(\alpha)}(\{{r}_{ij}\}), \end{equation} \tag{ A1 }$$

where

$$\begin{equation} M_{\mu_l,\nu_l}(\{{r}_{ij}\}) = \sum_{\boldsymbol{r}_{ij} \in \{{r}_{ij}\}} f_{\mu_l}(|\boldsymbol{r}_{ij}|) \, \underbrace{\boldsymbol{r}_{ij} \otimes \cdots \otimes \boldsymbol{r}_{ij}}_{\nu_l \; \mathrm{times}}, \end{equation} \tag{ A2 }$$

are the moment tensor descriptors. Here, the $f_{\mu_l}$'s are radial basis functions which vanish for $|\boldsymbol{r}_{ij}| > r_{\rm cut}$, where $r_{\rm cut}$ is a predefined cut-off radius. We characterize the degree of an MTP by its level. That is, an MTP of level d implies that the $B_\alpha$'s comprise all basis functions which satisfy $(2\mu_1 + \nu_1) + (2\mu_2 + \nu_2) + ... \le d$.

For all numerical examples, conducted in this work, we have used an MTP of level d = 16 with a cut-off radius $r_{\rm cut} = 5$ Å.

In the case of non-linear MTPs, as we use them in the current work, the radial basis functions $f_{\mu_l}$ depend on an additional parameter vector $\underline{c}_{\mu_l}$ as follows (cf [37])

$$\begin{equation} f_{\mu_l}(|\boldsymbol{r}_{ij}|;\underline{c}_{\mu_l}) = \sum_n c^n_{\mu_l} Q^n(|\boldsymbol{r}_{ij}|), \qquad \mathrm{with} \; Q^n(|\boldsymbol{r}_{ij}|) = T_n(|\boldsymbol{r}_{ij}|)(r_{\rm cut} - |\boldsymbol{r}_{ij}|)^2, \end{equation} \tag{ A3 }$$

where T_n is the nth order Chebysev polynomial and the function $(r_{\rm cut} - |\boldsymbol{r}_{ij}|)^2$ ensures a smooth transition to 0 as $|\boldsymbol{r}_{ij}| \rightarrow r_{\rm cut}$.

Hence, the basis functions have a non-linear dependence on $\underline{c}_{\mu_l}$. Therefore, we need to modify the active learning algorithm to take this non-linearity into account. For this purpose, under the assumption that the vector $\underline{\theta}$ now comprises the basis coefficients and the $\underline{c}_{\mu_l}$ parameters, we linearize equation (3) with respect to some initial $\bar{\underline{\theta}}$ to first order such that

$$\begin{equation} \boldsymbol{\epsilon}(\{{r}_i\};\underline{\theta}) - \boldsymbol{\epsilon}^{\rm qm}(\{{r}_i\}) \approx \boldsymbol{\epsilon}(\{{r}_i\};\bar{\underline{\theta}}) - \boldsymbol{\epsilon}^{\rm qm}(\{{r}_i\}) + \frac{\partial \boldsymbol{\epsilon}(\{{r}_i\};\underline{\theta})}{\partial \underline{\theta}}\bigg\vert_{\underline{\theta} = \bar{\underline{\theta}}} (\underline{\theta} - \bar{\underline{\theta}}) = L(\underline{\theta}) = 0. \end{equation} \tag{ A4 }$$

The system matrix associated with the new loss function L is then the n × m Jacobian matrix

$$\begin{equation} \underline{\underline{J}} = \begin{pmatrix} \frac{\partial \boldsymbol{\epsilon}(\{{r}_{ij}\}_1;\underline{\theta})}{\partial \theta_1}\Big\vert_{\underline{\theta} = \bar{\underline{\theta}}} & \cdots & \frac{\partial \boldsymbol{\epsilon}(\{{r}_{ij}\}_1;\underline{\theta})}{\partial \theta_m}\Big\vert_{\underline{\theta} = \bar{\underline{\theta}}} \\ \vdots & \ddots & \vdots \\ \frac{\partial \boldsymbol{\epsilon}(\{{r}_{ij}\}_n;\underline{\theta})}{\partial \theta_1}\Big\vert_{\underline{\theta} = \bar{\underline{\theta}}} & \cdots & \frac{\partial \boldsymbol{\epsilon}(\{{r}_{ij}\}_n;\underline{\theta})}{\partial \theta_m}\Big\vert_{\underline{\theta} = \bar{\underline{\theta}}} \end{pmatrix}, \end{equation} \tag{ A5 }$$

where n is the total number of atomic neighborhoods in the training set which is usually significantly larger than the size of the parameter vector m.

We proceed as in section 2 by finding an m × m submatrix $\underline{\underline{A}}$ of $\underline{\underline{J}}$ using the maxvol algorithm. The extrapolation grade γ of a neighborhood ${\{{r}_i^\ast\}}$ is now defined—analogously to equation (4)—as the maximum absolute element of the vector

$$\begin{equation} \begin{aligned} \underline{c} & = \underline{b}^{T} \underline{\underline{A}}^{-1}\\ & = \begin{pmatrix} \frac{\partial \boldsymbol{\epsilon}(\{{r}_{ij}^\ast\};\underline{\theta})}{\partial \theta_1}\Big\vert_{\underline{\theta} = \bar{\underline{\theta}}} & \cdots & \frac{\partial \boldsymbol{\epsilon}(\{{r}_{ij}^\ast\};\underline{\theta})}{\partial \theta_m}\Big\vert_{\underline{\theta} = \bar{\underline{\theta}}} \end{pmatrix} \begin{pmatrix} \frac{\partial \boldsymbol{\epsilon}(\{{r}_{ij}\}_1;\underline{\theta})}{\partial \theta_1}\Big\vert_{\underline{\theta} = \bar{\underline{\theta}}} & \cdots & \frac{\partial \boldsymbol{\epsilon}(\{{r}_{ij}\}_1;\underline{\theta})}{\partial \theta_m}\Big\vert_{\underline{\theta} = \bar{\underline{\theta}}} \\ \vdots & \ddots & \vdots \\ \frac{\partial \boldsymbol{\epsilon}(\{{r}_{ij}\}_m;\underline{\theta})}{\partial \theta_1}\Big\vert_{\underline{\theta} = \bar{\underline{\theta}}} & \cdots & \frac{\partial \boldsymbol{\epsilon}(\{{r}_{ij}\}_m;\underline{\theta})}{\partial \theta_m}\Big\vert_{\underline{\theta} = \bar{\underline{\theta}}} \end{pmatrix}^{-1}. \end{aligned} \end{equation} \tag{ A6 }$$

Since we perform structural relaxation and compute energy barriers, we need to fit the MTP with respect to energies and forces. Therefore, we compute the hyperparameters $\underline{\theta}$ by minimizing the loss functional

$$\begin{equation} L(\underline{\theta}) = C_{\rm e} \sum_{\{{r}_{i}^{\rm tr} \} \in \mathscr{T}} \left( \boldsymbol{\Pi}(\{{r}_{i}^{\rm tr} \}\underline{\theta}) - \Pi^{\rm qm}(\{{r}_{i}^{\rm tr} \}) \right) + C_{\rm f} \sum_{\{{r}_{i}^{\rm tr}\} \in \mathscr{T}} \sum_{\boldsymbol{r}_{i} \in \{{r}_{i}^{\rm tr} \}} \left( \boldsymbol{f}_{\rm r}(\{{r}_{i}^{\rm tr} \}\underline{\theta}) - \boldsymbol{f}^{\rm qm}_{\rm r}(\{{r}_{i}^{\rm tr} \}) \right) + C_{\rm r} (\underline{\theta}^{\mathsf{T}} \cdot \underline{\theta}), \end{equation} \tag{ A7 }$$

where $\boldsymbol{f}_{\rm r}$, $\boldsymbol{f}^{\rm qm}_{\boldsymbol{r}}$ are the MTP, respectively, the quantum mechanical forces on an atom $\boldsymbol{r}_i$ corresponding to a specific configuration $\{{r}_i^{\rm tr} \}$ in the training set $\mathscr{T}$. The C-constants are regularization parameters which we set to

$$\begin{align} C_{\rm e} = 1, \qquad C_{\rm f} = 0.01, \qquad C_{\rm r} = 10^{-6}, \end{align} \tag{ A8 }$$

throughout this work.