Improved description of atomic environments using low-cost polynomial functions with compact support

First applications of NNPs in computational chemistry date back as far as 1995 [15]. A considerable advance in the prediction of molecular forces and energies was made in 2007, when Behler and Parrinello [17] proposed to use HDNNP in combination with a fictitious decomposition of the potential energy $V$ from $N$ individual atoms:

$$\begin{equation} V = \sum_i^N V_i(\mathbf{G}_i). \end{equation} \tag{ 1 }$$

The atomic contribution $V_i$ depends on the chemical environment of atom $i$, which is described by a tuple of $N_j$ input values $\mathbf{G}_i = \{G_1, \dots, G_j \}$. The forces on the $k$th atom follow from the gradient:

$$\begin{equation} -\nabla_k V = - \sum_i^N \sum_j^{N_j} \frac{\partial V_i}{\partial G_{i,j}} \nabla_k G_{i,j}. \end{equation} \tag{ 2 }$$

The $\mathbf{G}_i$ are used as the input layer of a feed-forward neural network (figure 1) that yields the scalar $V_i$. The network itself is constituted by multiple layers of nodes, which are connected to each other by virtue of weights $a_{nm}$. Starting from a set of input values $(\mathbf{G}_i)$, the data is propagated to node $n$ of hidden layer $k$ as follows:

$$\begin{equation} y_n^k = f_a \left( b_n^k + \sum_m^{M_l} a_{nm}^{lk} y_m^l \right), \end{equation} \tag{ 3 }$$

where $l = k-1$, $y_n^k$ represents the value of the node and $f_a$ is a non-linear activation function. $b_n^k$ denotes the bias associated to $y_n^k$ and $M_l$ is the number of neurons in layer $l$. One independent neural network per element is used. Together with the biases $b_n^k$, the weights $a_{nm}$ are the fitting parameters of the HDNNP.

Zoom In Zoom Out Reset image size

The atomic environment descriptors $\mathbf{G}_i$ are commonly called SFs. In order to ensure permutational symmetry, the SFs are identical for all atoms of the same element. However, in order to account for varying chemical environments (bond lengths, van-der-Waals radii, relevant bond angles), they can differ between different elements. An appropriate choice of SFs is crucial in order to reliably discriminate between unique points on the PES.

In the case of Behler–Parrinello HDNNPs, the input layer consists of scalars obtained from a predetermined set of atom-centred SFs. The choice of SFs should discriminate between all relevant structural features, such that every distinct point on the PES is described by a unique tuple of SF values. In practice, this is achieved by combining a set of spherically symmetric and angle-dependent functions that have to be parametrised appropriately. Some examples are given in figure 2. For the radial part, Behler and Parrinello suggested:

$$\begin{equation} f\kern1.8pt_\mathrm{B}^\mathrm{rad} (r,r_s) = e^{-\eta(r-r_s)^2}, \end{equation} \tag{ 4 }$$

where $r$ is an interatomic distance and $r_s$ its reference point. An angular dependency can be introduced by a function of the form:

$$\begin{equation} f\kern1.8pt_\mathrm{B}^\mathrm{ang} (\vartheta,\lambda,\zeta) = \left(1 + \lambda \cos(\vartheta) \right)^\zeta, \end{equation} \tag{ 5 }$$

where $\vartheta$ is the angle formed by any three atoms, bound by $[0,\pi]$. This function appropriately has $\partial f\kern1.8pt_\mathrm{B}^\mathrm{ang} / \partial \vartheta \lvert_{\vartheta = 0} = \partial f\kern1.8pt_\mathrm{B}^\mathrm{ang} / \partial \vartheta \lvert_{\vartheta = \pi} = 0$. For $\lambda = 1$, maximum and minimum of the function lie at 0 and $\pi$, respectively; for $\lambda = -1$ the inverse holds. Choosing $\zeta \gt 1$ allows to contract the function, resulting in a more rapid decay. Values of $\zeta \lt 1$ are not possible, since this would result in non-zero derivatives at the boundaries.

Zoom In Zoom Out Reset image size

While the form of $f{\,\,}^\mathrm{ang}_\mathrm{B} (\vartheta,\lambda,\zeta)$ guarantees zero derivatives at the boundaries for appropriate choices of $\lambda$ and $\zeta$, the exponential in $f{\,\,}^\mathrm{rad}_\mathrm{B} (r,r_c)$ is finite everywhere. Therefore, a cutoff function with cutoff radius $r_c$ is introduced, which ensures that the radial term and its derivatives are zero for all $r \gt r_c$. Common choices of cutoff function $f_c(r,r_c)$ include:

$$\begin{align} f_c(r,r_c) = \begin{cases} \frac{1}{2} \left[ \cos\left(\frac{\pi r}{r_c}\right) + 1 \right] \quad &\mathrm{for} \quad r \leqslant r_c \\ 0 \quad &\mathrm{for} \quad r \gt r_c \end{cases}, \end{align} \tag{ 6 }$$

or

$$\begin{align} f_c(r,r_c) = \begin{cases} \tanh^3 \left(1 - \frac{r}{r_c}\right) \quad &\mathrm{for} \quad r \leqslant r_c \\ 0 \quad &\mathrm{for} \quad r \gt r_c \end{cases}. \end{align} \tag{ 7 }$$

With this choice of functions, the radial SF $G_i$ centred on an atom $i$ is given by a sum over all its neighbours j that are within $r_c$ of $i$:

$$\begin{equation} G_i^\mathrm{rad}(r_c,r_s) = \sum_{j \ne i} f\kern1.8pt^\mathrm{rad}_{\mathrm{B}} (r_{ij},r_s) f_c(r_{ij},r_c). \end{equation} \tag{ 8 }$$

Note that such a product may have a width and a maximum different from $f\kern1.8pt^\mathrm{rad}_\mathrm{B}$ alone, see also figure 2. The angular SF is constructed from a product of angular and radial terms, with the sum encompassing all neighbours $j$ and $k \ne j$. Depending on the radial term, two types of angular SFs can be distinguished, namely narrow $n$,

$$\begin{align} G_i^\mathrm{ang.n} (r_c,r_s,\lambda,\zeta) & = 2^{(1-\zeta)} \sum_{j \ne i, k \gt j} f\kern1.8pt^\mathrm{rad}_{\mathrm{B}} (r_{ij}, r_s) f\kern1.8pt^\mathrm{rad}_{\mathrm{B}} (r_{ik}, r_s) f\kern1.8pt^\mathrm{rad}_{\mathrm{B}} (r_{jk}, r_s) \end{align} \tag{ 9 }$$

$$\begin{align} &\qquad\quad\quad\qquad\quad\times f\kern1.8pt^\mathrm{ang}_{\mathrm{B}} (\vartheta_{ijk},\lambda,\zeta) f_c(r_{ij},r_c) f_c(r_{ik},r_c) f_c(r_{jk},r_c), \end{align} \tag{ 10 }$$

which also constrains the interatomic distance between atoms $j$ and $k$, and the wide form w:

$$\begin{equation} G_i^\mathrm{ang.w} (r_c,r_s,\lambda,\zeta) = 2^{(1-\zeta)} \sum_{j \ne i, k \gt j} f\kern1.8pt^\mathrm{rad}_{\mathrm{B}} (r_{ij}, r_s) f\kern1.8pt^\mathrm{rad}_{\mathrm{B}} (r_{ik}, r_s) f\kern1.8pt^\mathrm{ang}_{\mathrm{B}} (\vartheta_{ijk},\lambda,\zeta) f_c(r_{ij},r_c) f_c(r_{ik},r_c), \end{equation} \tag{ 11 }$$

the evaluation of which is more computationally expedient due to the lack of a radial term for $jk$. These SFs have successfully been applied in a wide variety of contexts. Note that the distribution of SF values at $\lambda = -1$ is not symmetric with respect to $\lambda = +1$ due the presence of a radial function [21]. Figure 3 displays the dependency of BPSFs on the position of atomic neighbour(s).

Zoom In Zoom Out Reset image size

Recently, Singraber et al have introduced a simple expression for $f_c$ based on polynomials [27]:

$$\begin{align} f_c^\mathrm{poly} (r,r_c) = \begin{cases} f\kern1.8pt^\mathrm{poly2} \left( \frac{r}{r_c} \right) \quad &\mathrm{for} \quad 0 \leqslant r \leqslant r_c \\ 0 \quad &\mathrm{for} \quad r \gt r_c, \end{cases} \end{align} \tag{ 12 }$$

with the polynomial:

$$\begin{equation} f\kern1.8pt^\mathrm{poly2} (x) = x^3 ( x (15 - 6x) - 10 ) + 1, \end{equation} \tag{ 13 }$$

with $f\kern1.8pt^\mathrm{poly2}(0) = 1$ and $f\kern1.8pt^\mathrm{poly2} (1) = 0$. At $x=0$ and $x=1$, equation (13) has continuous derivatives up to second order, which ensures a continuous, smooth transition to outside of the compact set. Similar polynomials with continuous derivatives up to even higher order can be easily constructed [27], but their use does not result in any practical benefit [18]. The main advantage of polynomial cutoff functions as introduced in [27] lies in their expedient evaluation compared to expressions based on hyperbolic tangent or cosine, reducing the overhead due to the cutoff function. However, the computation of costly exponentials is still needed for all radial components. Hence, a grouping strategy was proposed in [27] which avoids multiple evaluations of the same exponential term, resulting in considerable speedups.

On the one hand, the summation over the arguments of $G^\mathrm{ang}$ still involves a substantial number of exponential functions that have to be explicitly evaluated even when a grouping strategy is used. On the other hand, the present choice of $f\kern1.5pt^\mathrm{ang}$ implies that, in order to specifically sample angles which are not centred close to either 0 or $\pi$, a combination of several values of $\lambda$ and $\zeta$ is necessary, since introducing an angular shift $\vartheta = 0$ would inevitably violate the boundary condition $\partial f\kern1.5pt^\mathrm{ang} / \partial \vartheta = 0$ at the bounds of the angular domain. In the following, we will attempt to overcome this issue by constructing a simple set of SFs based on polynomials with compact support which can be used to indiscriminately describe both angular and radial dependencies.

In order to do so, we generalise equation (12) to a generic form, expressed in terms of the boundaries $x_\mathrm{min},x_\mathrm{max}$ of the underlying domain. In particular, for our polynomial SFs $f\kern1.7pt_\mathrm{PSF}$, we also allow for arguments $x \lt 0$:

$$\begin{align} f\kern1.6pt_\mathrm{PSF} (x,x_0,\Delta x) = \begin{cases} f\kern1.8pt^\mathrm{poly2} \left( \frac{x-x_0}{\Delta x} \right) \quad &\mathrm{for} \quad x_0 \leqslant x \leqslant x_0 + \Delta x \\ f\kern1.8pt^\mathrm{poly2} \left( \frac{x_0-x}{\Delta x} \right) \quad &\mathrm{for} \quad x_0 - \Delta x \leqslant x \lt x_0 \\ 0 \quad &\mathrm{elsewhere,} \end{cases} \end{align} \tag{ 14 }$$

where $x_0 = \frac{1}{2}(x_\mathrm{max}+x_\mathrm{min})$ and $\Delta x = \frac{1}{2}(x_\mathrm{max}-x_\mathrm{min})$. The resulting function is symmetric with respect to $x_0$ and has continuous derivatives — up to an order determined by the construction of $f\kern1.8pt^\mathrm{poly2}$ — on all $\mathbb{R}^3$. It is possible to further influence the decay of the radial function by appropriately modifying its argument such that the SF decays more rapidly around its maximum while also approaching 0 more slowly, thus breaking the point symmetry around $f\kern1.7pt_\mathrm{PSF} = 0.5$. We therefore define the asymmetric PAS as:

$$\begin{align} f\kern1.2pt_\mathrm{PAS} (x,x_0,\Delta x) = \begin{cases} f\kern1.8pt^\mathrm{poly2} \left( \frac{2(x-x_0)}{\Delta x} - \left( \frac{x-x_0}{\Delta x} \right)^2 \right) \quad &\mathrm{for} \quad x_0 \leqslant x \leqslant x_0 + \Delta x \\ f\kern1.8pt^\mathrm{poly2} \left( \frac{2(x_0-x)}{\Delta x} - \left( \frac{x_0-x}{\Delta x} \right)^2 \right) \quad &\mathrm{for} \quad x_0 - \Delta x \leqslant x \lt x_0 \\ 0 \quad &\mathrm{elsewhere,} \end{cases} \end{align} \tag{ 15 }$$

which offers the same benefits as $f\kern1.7pt_\mathrm{PSF}$ for arguments ∈ [0, 1]. By construction, $f\kern1.2pt_\mathrm{PAS}$ resembles the shape of the product of a local, narrow Gaussian function with a cutoff function and can therefore be used to adapt existing BPSF-setups to polynomial form. The functional forms of $f\kern1.2pt_\mathrm{PAS}$ and $f\kern1.7pt_\mathrm{PSF}$ are depicted in figure 4.

Zoom In Zoom Out Reset image size

The corresponding radial SFs can now be constructed as:

$$\begin{equation} G_i^\mathrm{rad/P} = \sum_{j \ne i} f_\mathrm{P} (r_{ij},r_s,\Delta r), \end{equation} \tag{ 16 }$$

with $f_\mathrm{P}$ either $f\kern1.7pt_\mathrm{PSF}$ or $f\kern1.2pt_\mathrm{PAS}$. Their narrow angular counterparts read:

$$\begin{equation} G_i^{\mathrm{ang.n/P}} = \sum_{j \ne i, k\gt j} f_\mathrm{P} (r_{ij},r_s,\Delta r) f_\mathrm{P} (r_{ik},r_s,\Delta r) f_\mathrm{P} (r_{jk},r_s,\Delta r) f\kern1.7pt_\mathrm{PSF} (\vartheta_{ijk},\vartheta_0,\Delta \vartheta), \end{equation} \tag{ 17 }$$

while the wide angular functions become:

$$\begin{equation} G_i^\mathrm{ang.w/P} = \sum_{j \ne i, k\gt j} f_\mathrm{P} (r_{ij},r_s,\Delta r) f_\mathrm{P} (r_{ik},r_s,\Delta r) f\kern1.7pt_\mathrm{PSF} (\vartheta_{ijk},\vartheta_0,\Delta \vartheta), \end{equation} \tag{ 18 }$$

where the radial component $f_\mathrm{P}$ can again be described by both $f\kern1.8pt_\mathrm{PAF}$ or $f\kern1.7pt_\mathrm{PSF}$. Both types of angular SF can be centred anywhere within $[0,\pi]$, provided that $\partial f\kern1.7pt_\mathrm{PSF} / \partial \delta \lvert_{\vartheta = 0} = 0$ and $\partial f\kern1.7pt_\mathrm{PSF}/\partial \delta \lvert_{\vartheta = \pi} = 0$. This enables chemically intuitive choices of angular maxima, i.e. 60^∘ for a double bond or 109, 4^∘ for a carbon-carbon single bond. Figure 5 shows the dependency of PAS-type SFs on the position of atomic neighbour(s).

Zoom In Zoom Out Reset image size

Equations (14) and (15) do not involve any exponentials, and the most expensive arithmetic operation, an inverse cosine, has a modest computational footprint.

All HDNNPs were trained using the n2p2 package [26, 27, 31]. HDNNPs for the prediction of enthalpies of formation in the QM9 [30] database were obtained using the protocol of [21] with five-fold cross validation, but without dataset normalisation. For the remaining setups, 10 training runs per system are carried out with an independent seed for the random number generator each. If the value $G_i$ of a given SF never exceeds 10⁻³ over the whole test set, the function is discarded before training. The choice of $\mathbf{G}$ is detailed in the supporting information (available online at stacks.iop.org/MLST/2/035026/mmedia). Results discussed in this text refer to the best set of weights out of 10 independent runs (5 for the QM9 set), which are defined as the weights that yield the lowest force root mean-square deviation (RMSD) with respect to a test set constituted by 10% of training-set structures that were randomly removed before training. If several sets of weights yield force RMSDs that are identical to a few %, the set that displays the lowest relative error in energy predictions is used.

The systems investigated here are depicted in figure 6. The training sets for water and copper sulphide correspond to those published in [5, 26], respectively. Forces and energies for an isolated DMABN molecule, liquid ethyl benzene and anisole (16 molecules in a periodic box) were obtained from high-temperature Car-Parrinello MD runs performed with the CPMD code [32] using the SCAN exchange-correlation functional [33]. The amino group of DMABN was rotated using a slowly growing harmonic restraint, giving rise to a set of non-equilibrium conformers. Details on the computational setup can be found in the supporting information.

Zoom In Zoom Out Reset image size

Test sets consists of 10% randomly chosen configurations that are removed from the training set before training and which are then used to validate the predictive power of the NNP.

All MD simulations using HDNNPs were carried out using the n2p2-interface [27, 31] to the LAMMPS code [34, 35]. System setups are detailed in the supporting information.

We first assess the accuracy that can be attained with polynomial SFs by comparing the efficiency and accuracy of both training and productive simulations using identically set up HDNNPs. Two setups of polynomial SFs will be used: one, PSF, consisting entirely of symmetric functions of the form 14, and one mixed setup, PAS, including both symmetric and asymmetric radial components. We will compare against a reference setup using BPSFs, which has successfully been used in [26] to describe the PES of water and several ice phases. PSF and PAS were constructed by plotting the BPSFs and using their peaks and the full width at half maximum as rough indicators for the polynomial input. For sake of simplicity, the only symmetric radial terms (PSF-type) in the PAS set were those with maximum width; all others were chosen to be of the asymmetric type.

First, we compare the learning curves of both methods to assess the efficiency of the training procedure. Here, we define the learning curve as the RMSD between trained energy and forces and their reference as a function of the number of training steps (epochs). Curves for Behler–Parrinello type and polynomial SF setups are plotted in figure 7 and are virtually indistinguishable, suggesting that both types of symmetry functions can be trained with equal efficiency. RMSDs at the last training step are comparable between different methods, with the largest remaining force deviation being observed for the PSF setup. Conversely, the lowest force RMSD is observed for the PAS setup. BPSFs yield a training error intermediate between PAS and PSF data.

Zoom In Zoom Out Reset image size

The same conclusion holds when analysing the distribution of both force and energy errors within the test set: Force RMSDs are 43.7 meV Å⁻¹ for BPSFs, 48.5 meV Å⁻¹ for PSF and 43.2 meV Å⁻¹ for PAS, respectively. Energy RMSDs range from 1.17 meV/atom for BPSFs to 0.95 meV/atom for PAS. In practical applications, these values could be considered identical. In figure 7, we histogram both the relative force error for HDNNPs trained with Behler–Parrinello type and polynomial SFs. The resulting distributions for the test set compare well between PSF, PAS and BPSF: All relative errors show a rapid decay, indicating an excellent quality of the fit. PAS and BPSF show very similar relative error distributions, whereas the distribution of PSF skews very lightly towards the right; however, we will show that this is of no practical relevance. Figure 8 shows predicted forces against their reference values; ideally, all values should lie on a diagonal. Use of PAS and PSF improves the accuracy of the most extreme force values at either end of the graph, which lie closer to the diagonal than for BPSF. This indicates that, although force RMSDs differ only little between the different SF setups, qualitative improvements to conventional BPSFs are possible by use of a PSF or PAS set.

Zoom In Zoom Out Reset image size

Figure 9 provides a practical comparison of the predictive power of the different HDNNPs. Both the radial distribution function (left panel) as well as the mean square displacement as a function of time, from which the diffusion coefficient can be derived, (right panel) are practically indistinguishable between BPSF, PAS and PSF. These results further support that polynomial SFs are a viable alternative to the commonly used Behler–Parrinello type functions. The excellent agreement between radial distribution functions and mean-squared displacements obtained from MD suggest that the small differences in force RMSD over the test set do not entail significant changes in the quality of the predictions made by the HDNNP. Notably, the mean-squared displacements (MSD) for PSF and BPSF are virtually identical.

Zoom In Zoom Out Reset image size

All SF setups used for water were based on the optimised, hand-selected set of [5]. Such a selection by hand may not always be possible. In the following, we shall further validate the accuracy of PSF and PAS by comparing their performance to the accuracy of a semi-automatically generated set of BPSFs for copper sulphide. It has been show in [26] that a HDNNP with BPSFs can be successfully used to observe structural phase transitions in this compound.

Figure 10 compares absolute and relative force errors obtained using BPSFs of [26] as well as our PSF and PAS. Force RMSDs are 51.9 meV Å⁻¹ for BPSF, 56.7 meV Å⁻¹ for PSF and 52.3 meV Å⁻¹ for PAS, which corresponds roughly to a 10% difference between PSF and BPSF. Again, the performance of BPSF and PAS is on par, whereas use of PSF leads to slightly larger errors, which is also reflected in the relative error distribution skewing lightly to the right. Overall, however, this effect is expected to remain negligible in practical applications. Deviations from the diagonal in the graph comparing prediction and reference forces (figure 10, right panel) do not exhibit visible differences, further indicating that the general performance of all SF types investigated here remains comparable and that the differences in force RMSDs observed here are not expected to influence the quality of the predictions in a significant manner. Energy RMSDs remain similar for all methods, spanning a range from 0.99 meV/atom for PAS to 1.06 meV/atom for BPSF.

Zoom In Zoom Out Reset image size

In the following, we shall compare the performance of BPSFs generated according to [24] with our PSFs using an isolated (gas-phase) DMABN molecule and liquid ethyl benzene. We note that optimal SF setups could possibly be found for all SF types investigated here; however, since this is a time-intensive procedure, the main purpose of this comparison is to compare the out-of-the-box performance of some generic sets of SFs.

To this end, we introduce a simple generation scheme for radial SF parameters given a minimal and maximum width, $\Delta_\mathrm{min}$ and $\Delta_\mathrm{max}$, and a maximum cutoff radius $r_c$. For $N_0$ SFs which are symmetric w.r.t. to the origin, we choose the width $\Delta^i_0$ parameters of the $i$th out of $N_0$ PSF or PAS as:

$$\begin{align} \Delta^0_0 & = r_c \end{align} \tag{ 19 }$$

$$\begin{align} \Delta^i_0 & = \Delta_\mathrm{max} - \frac{(i-1)(\Delta_\mathrm{max} - \Delta_\mathrm{min})}{N_0-2} \quad ; \quad i \geqslant 1 \end{align} \tag{ 20 }$$

For a given number $N_s$ of SFs which are centred at $r_s^i \ne 0$, we obtain shifts and widths $\Delta^i_s$ as:

$$\begin{align} \Delta^i_s & = \frac{2r_c}{N_s+1} \end{align} \tag{ 21 }$$

$$\begin{align} r_s^i & = \frac{(i-1)\Delta^i_s}{2}. \end{align} \tag{ 22 }$$

This partitioning is reminiscent to the one introduced by Gastegger et al in [21]. Angles were chosen to reflect common structural features encountered in organic chemistry: linearity for triple bonds, 60^∘ and 109.1^∘ angles to describe double- and single carbon–carbon bonds, respectively. This resulted in the following angular functions: A wider set centred on $\vartheta_0^{1} = 0^\circ$, $\vartheta_1^{1} = 90^\circ$ and $\vartheta_3^{1} = 180^\circ$ with a width of $\Delta \vartheta^{1,2,3} = 90^\circ$ each; and an additional set centred on $\vartheta_0^{4} = 60^\circ$, $\vartheta_0^{5} = 120^\circ$ and $\vartheta_0^6 = 109.1^\circ$, corresponding to a more specific sampling of environments typical of double- and single bonds, with a width of $\Delta \vartheta^{4,5,6} = 60^\circ$ each. For sake of comparison, the cutoff function distance for BPSFs and the maximum extent of PAS and PSF were set to 15 a.u.

Structures in the DMABN-set cover regions of the PES thermally accessible at 1320 K, including a forced rotation around the amino group along a non-equilibrium path. This results in many forces with large magnitude which are well outside of the equilibrium regime in which simple harmonic potentials can also yield good approximations. The DMABN set therefore serves as an assessment for the transferability of our polynomial SFs to structures which are considerably off-equilibrium.

Force RMSDs on the test set amount to 63.2 meV Å⁻¹ for BPSF, 56.6 meV Å⁻¹ for PSF and 54.0 meV Å⁻¹ for PAS. Note that this amounts to a difference exceeding 15% between worst (BPSF) and best (PAS) set. However, the distribution of relative errors remain comparable between all SFs, as evidenced by the force histogram in the right panel of figure 11. Energy RMSDs range from 0.64 meV/atom for PSF to 0.72 meV/atom for BPSF, with the overall error remaining low for all types of functions. The force predictions obtained here imply that for a generic test set containing a significant amount of off-equilibrium structures, generic PAS and PSF sets are able to outperform BPSFs.

Zoom In Zoom Out Reset image size

We will now show that the same considerations also hold for an organic liquid in which van der Waals interactions dominate. For a test set of liquid ethyl benzene, using the same SF setups as for DMABN, we find force prediction RMSDs of 89.5 meV Å⁻¹ for BPSFs, 68.9 meV Å⁻¹ for PSF and 61.7 meV Å⁻¹ for PAS, respectively. This corresponds to a difference in accuracy approaching 50% between BPSF and PAS. Whereas these considerable differences are not evident from the diagonal force-to-force plot in the left panel of figure 12, they become visible in the force histogram in the right panel: Histograms for PSF and PAS are very similar, with the only significant difference being the peak around 0%, but BPSFs skew considerably to the right, and much more so than what was observed for PSF in the case of water or Cu₂S. Energy RMSDs are about 0.05 meV/atom for all SFs studied here, which is of excellent accuracy. For liquid anisole, all SF setups result in about equally accurate predictions, as evidenced by figure 13: observed force prediction RMSDs amount to 100.8 meV Å⁻¹ for BPSFs and a comparable 106.5 meV Å⁻¹ for PSF, with the lowest deviation again being observed at 98.7 meV Å⁻¹ for PAS. All SFs yield energy RMSDs that can be considered identical, within a small range from 0.08 meV/atom for BPSF to 0.11 meV/atom for PAS and PSF. Overall, the trends obtained for force predictions in organic molecules indicate that our generic PAS and PSF setups on average exhibit superior transferability between systems, while at the same time retaining lower force errors on the test sets compared to generic BPSFs. Importantly, whereas BPSFs performed considerably worse than polynomial functions for liquid ethyl benzene, neither PAS nor PSF showed a considerably inferior performance with respect to BPSF in any of the other systems investigated here. Energy predictions are much more robust with respect to the choice of system, they remain of comparable accuracy for all SFs investigated here.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to pinpoint possible improvements when solely energies are to be treated, we have constructed a set of wACSF based on PAS and have used this set to predict enthalpies of formation in the QM9 database. In order to render the training more challenging, and contrary to [21], enthalpies in the dataset were neither normalised nor rescaled. PAS were set up according to the procedure outlined in section 4.3, employing a maximum width of 8 Å and using a random choice of 10 000 structures for training (and testing) as outlined in [21], the remainder being used as a validation set. Mean absolute errors (MAE) over predicted enthalpies of formation for the remaining ≈124 000 molecules of the database were averaged over 5 independent partitionings into train, test and validation sets as in [21].

As shown in table 1, while performance on the training set is comparable between the best BPSF set reported in [21] and PAS, the quality of the predictions on test and validation sets within a given selection of 10 000 test- and training structures considerably improves when PAS are used. For the random partitionings into training and test data studied here, errors obtained from HDNNPs trained with PAS are about 5-fold smaller than their BPSF counterparts. It should be noted that once BPSF are used to train a HDNNP on an appropriately normalised and rescaled training set, their performance approaches that of PAS on a non-normalised test set [21]. This indicates that the robustness of the predictions substantially increases when PAS are used, as test and validation error are substantially lower even when the dataset is not subjected to any normalisation.

The choice of SFs has an influence on both, the computational footprint of training as well as the cost of force predictions during production runs. For the former, if SFs can be calculated for all structures in the training set and then stored, it is mainly the number of SFs that influences performance: more SFs imply more weights to be optimised and as such result in higher execution times. Memory requirements will grow as well as the number of SFs increases. Even if a (generic) set of two SFs might comprise the same number of input functions, pruning of functions that never exceed a certain threshold (typically 10⁻³) can result in different set sizes during training. In productive MD runs, SFs have to be recalculated for every new structure and hence, performance will not only be determined by their sheer number, but also by the number and type of floating point operations involved in their computation. In the following, we will analyse both performance-related aspects.

Table 2 lists the cost of running 100 steps of MD for a periodic box containing 360 water molecules using the HDNNPs investigated here both for a straightforward calculation of all SFs as well as using an optimised algorithm tailored for distributed memory parallelisation on core processing units (CPUs) as reported in [27]. The former serves as a measure for the cost of the floating-point operations associated to every SF type; this corresponds to the limit of all N out of N SFs being independent and differently parametrised. In comparison with the former, the latter indicates the benefit of underlying algorithmic optimisations; this however requires that a subset of the N SFs contain identical parameters and cutoff functions in order to benefit from a grouping and caching strategy. The effective speedups obtained by using PSF or PAS for the water setup reach almost a factor of 5 when no grouping or caching strategy is applied. This indicates that the number of expensive floating point operations is significantly reduced when purely polynomial SFs are used. Vice versa, SF grouping and caching has a much lower influence of execution times of PSF and PAS, where those optimisations reduce the computational cost by about 15%. As would be expected by the larger number of terms involved in the evaluation of PAS, their computation is slightly more expensive than the evaluation of PSF. However, the difference is small enough to be negligible in practice.

In contrast, BPSFs—for which grouping and caching strategies were developed in the first place—benefit from speedups of up to a factor of 3 when grouping and caching strategies are introduced. Note that, by construction, PAS and PSF are not intended to require use of optimisation strategies, which is reflected in comparably modest speedups upon introduction of grouping and caching. Still, even in those cases, PSF and PAS are almost a factor of 2 faster than BPSF. In addition, it should be noted that PAS yield the lowest RMSD on water with the smallest total number of SFs; they therefore do not only lower the computational cost, but also the memory footprint with respect to BPSFs.

Table 3 compares the numbers of SFs effectively used during training after pruning all SFs that never exceed a magnitude of 10⁻³. This provides a measure for computational cost and memory requirements. Apart from ethyl benzene, in which significantly more SFs of BPSF-type were removed during pruning, numbers are comparable between BPSF, PSF and PAS. It should be noted that all systems except water were trained with automatically generated sets of SF parameters, which also accounts for the different extent of functions with maximum amplitude of $\lt\kern-1.6pt10^{-3}$ in these setups. In this context, it is interesting to note that in spite of the structural similarities between anisole and ethyl benzene, for BPSFs, there is a significant difference in the number of pruned SFs. For PAS and PSF, no such difference can be observed; this can also account for the larger errors in force RMSDs that are observed when ethyl benzene is trained with a set of BPSFs. Within the hand-picked set of SFs for water, differences between PSF and PAS are notable, with the latter allowing for a reduction of about 15% with respect to the former. BPSFs lie in between. This suggests that, for hand-selected sets, PAS allow for an even smaller number of SFs to be used.

Overall, the performance of PSF and BPSFs was comparable for all systems studied here. Force RMSDs for water and Cu₂S were slightly higher for PSF; however, this was not reflected in the dynamic and structural properties of water obtained by 2 ns of HDNNP-MD, where results using PSF and BPSFs showed excellent agreement. On the other hand, for DMABN and liquid ethyl benzene, PSF considerably outperformed BPSFs. For all systems, PAS were able to outperform either PSF or BPSF. In particular, when used in the wACSF scheme [21], PAS greatly improve performance over BPSF for data sets that are not normalised, indicating superior robustness. We therefore strongly recommend PAS for future HDNNP setups, since they allow to speed up productive MD runs by a factor of 1.8 and further improve accuracy of the HDNNPs. In particular, due to their floating point cost being lower by about a factor of almost 5, we particularly recommend use of PAS for implementations where caching and grouping strategies are not feasible. This can, for instance, be the case when offloading computation of SFs to graphical processing units, where the resolution of if-statements and loop breaks associated to optimisation schemes can considerably slow down computation with respect to purely arithmetic operations. Not least, thanks to their small computational footprint, PAS can also be used in programs that have not undergone extensive optimisation, and their simple structure facilitates implementation and calculation of derivatives up to an arbitrary order defined by the underlying polynomial.