4G-HDNNP

The overall structure of the 4G-HDNNP is shown schematically in Fig. 2 for an arbitrary binary system. Like in 3G-HDNNPs the total energy consists of a short-range part, which, as we will see below, requires in addition non-local information, and an electrostatic long-range part, which is not truncated,

Fig. 2

Schematic structure of a 4G-HDNNP for a binary system.

For a binary system containing N_a atoms of element a and N_b atoms of element b the total energy consists of a short-range energy E_short, which is a sum of atomic energies E_i, and a long-range electrostatic energy E_elec computed from atomic charges Q_i. The atomic charges are determined by a charge equilibration method using environment-dependent atomic electronegativies χ_i expressed by atomic neural networks (red). These charges are then used to calculate the electrostatic energy and in addition serve as non-local input for the short-range atomic neural networks (blue) yielding the E_i. The geometric atomic environments are described by atom-centered symmetry function vectors G_i, which depend on the Cartesian coordinates R_i of the atoms and serve as inputs for the atomic neural networks.

Like in the CENT approach the atomic electronegativities are local properties defined as a function of the atomic environments using atomic neural networks. As in 2G- and 3G-HDNNPs there is one type of atomic neural network with a fixed architecture per element in the system making all atoms of the same type chemically equivalent, while the specific values of the electronegativities depend on the positions of all neighboring atoms inside a cutoff sphere of radius R_c. The positions of the neighboring atoms inside this sphere are specified by a vector G_i of atom-centered symmetry functions⁴¹, which ensures the translational, rotational and permutational invariance of the electronegativities.

To predict the atomic charges, which are represented by Gaussian charge densities of width σ_i taken from the covalent radii of the respective elements, a charge equilibration scheme⁴² is used. In this scheme, the charge is distributed among the atoms in an optimal way to minimize the energy expression

Highly optimized algorithms are available for systems of linear equations, which can be efficiently solved for small and medium-sized systems containing up to about ten thousand atoms in a few seconds on modern hardware. For larger systems the cubic scaling of the standard algorithms can pose a bottleneck. In that case one could resort to using iterative solvers for which the most expensive part of each iteration is a matrix vector multiplication involving the matrix A. This corresponds to the evaluation of the electrostatic potential at each atoms position for which numerous low-complexity algorithms, such as fast multipole methods, are known. In this way it is possible to reduce the effort from cubic to nearly linear scaling providing access to very large systems.

Overall, this process is like in the CENT³², but the main difference is in the training process. While in CENT only the error with respect to the DFT energies is minimized, the atomic charges obtained during the charge equilibration process serve merely as intermediate quantities, which do not have a strict physical meaning. In the 4G-HDNNP proposed in this work, the charges are trained directly to reproduce reference charges from DFT, which therefore are qualitatively meaningful although one should be aware that atomic partial charges are not physical observables and different partitioning schemes can yield different numerical values⁴³.

Once the atomic electronegativities have been learned, a functional relation between the atomic structure and the atomic partial charges is available. The intermediate global charge equilibration step ensures that these charges depend on the atomic positions, chemical composition and total charge of the entire system, and thus in contrast to 3G-HDNNPs non-local charge transfer is naturally included.

In a second step, the local atomic energy contributions yielding the short-range energy according to

The short-range atomic neural networks are then trained to express the remaining part of the total energy E_ref according to

In summary, in contrast to the CENT method, the short-range interactions are not described through the charges resulting from the charge equilibration process but are described by separate short-range neural networks, which enables a more accurate description of the total energy.

Overview of test systems

In the following subsections we demonstrate the limitations of ML potentials based on local properties only and show how they can be overcome by the 4G-HDNNP. For this purpose we use a set of non-periodic and periodic systems, which cover a wide range of typical situations in chemistry and materials science. The non-periodic systems consist of a covalent organic molecule, a small metal cluster and a cluster of an ionic material covering very different types of atomic interactions. These examples demonstrate the simultaneous applicability of a single 4G-HDNNP to systems of different total charges and the correct description of long-range charge transfer and the associated electrostatic energy. As a periodic system we have chosen a small gold cluster adsorbed on a MgO(001) slab, which is a prototypical example for heterogeneous catalysis. We show that in contrast to established ML potentials, the 4G-HDNNP is able to reproduce the change in adsorption geometry of the cluster if dopant atoms are introduced in the slab far away from the cluster. In all cases, the 4G-HDNNP PES is very close to the results obtained from DFT.

While in theses examples we do not explicitly investigate the transferability of the potentials to different systems, we expect that the 4G-HDNNP in general provides an improved transferability compared to 2G and 3G ML potentials due to the underlying physical description of the global charge distribution and the resulting electrostatic energy. This expectation is supported by the fact that even traditional charge equilibration schemes with constant electronegativities are known to work well across different systems⁴⁴. Furthermore, for the related CENT approach a broad transferability has already been demonstrated for different atomic environments³³.

A benchmark for organic molecules

The first model system we study is a linear organic molecule consisting of a chain of ten sp-hybridized carbon atoms terminated by two hydrogen atoms as shown in Fig. 3a. Molecules of this type have been studied before in electronic structure calculations^45–47. For this molecule we will now demonstrate the applicability of 4G-HDNNPs to systems with long-range charge transfer induced by protonation, which changes the total charge and the local structure in a part of the system. Since the majority of existing machine learning potentials rely on local structural information only without explicit information about the global charge distribution and total charge, they are not simultaneously applicable to both neutral and charged systems.

Fig. 3

Charge redistribution in organic molecules.

a DFT-optimized structures of C₁₀H₂ (left) and C₁₀H$$ {}_3^{+}$$ (right) with atom IDs. Carbon and hydrogen atoms are colored in gray and white, respectively. The dashed circle shows the cutoff radius of the left carbon atom defining its chemical environment. b shows the atomic partial charges obtained from DFT. The unscaled and scaled 3G-HDNNP charges are displayed in c, while the 4G-HDNNP charges are shown in d.

This is different for 4G-HDNNPs, which naturally include the correct long-range electrostatic energy for any global charge present in the training set. Because of the protonation of the terminal carbon atom, its hybridization state changes to sp² and the electronic structure of the resulting C₁₀H$$ {}_3^{+}$$ cation is modified even at very large distances along the whole molecule, which is reflected in the differences of the DFT charges of the molecules in Fig. 3b, which have been structurally optimized by DFT. The geometries of both molecules are given in the supplementary tables.

Using a data set containing both molecules, we have constructed 2G-, 3G-, and 4G-HDNNPs using a cutoff radius R_c = 4.23 Å as illustrated by the circle in Fig. 3a for the example of the left carbon atom. In Fig. 3c we show the atomic partial charges obtained with the 3G-HDNNP in two forms: first as unscaled charges directly obtained from the atomic neural network fits without any constraint for the correct total charge of the system, and second rescaled to ensure total charges of zero or one, respectively. It can be seen that the scaling process does not significantly improve the 3G-HDNNP charges.

The atoms in the left half of the molecule are far from the added proton such that their atomic environments differ only slightly due to the DFT geometry optimization. In addition, in the training set a lot of basically identical environments but different atomic charges are present for these atoms, which results in high fitting errors due to the contradictory information. As a consequence the neural networks assign averaged charges to these atoms, which differ qualitatively from the DFT reference charges of both systems. For instance, the 3G-HDNNP partial charges on atom 2, i.e., the left carbon atom, are almost identical in both molecules although they are very different in DFT. Note that the predicted charges of atoms 1-6 in C₁₀H₂ and C₁₀H$$ {}_3^{+}$$ would be even exactly identical if the latter molecule would not have been relaxed after protonation. The charges obtained with the 4G-HDNNP shown in Fig. 3d, on the other hand, match the DFT charges very accurately for both molecules, as they can be distinguished in this method.

The inaccurate charges obtained with the 3G-HDNNP lead to a poor quality of the potential energy surface, and the same is observed for the short-range only 2G-HDNNP. In Table 1 we compare the errors of the total energies as well as the mean errors of the atomic charges and forces of all HDNNP generations for the DFT-optimized structures. It can be seen that the errors of all quantities obtained for the 4G-HDNNP are much lower than for the 2G- and 3G-HDNNPs. Further, we note that in several cases the energies obtained by the 3G-HDNNP are even worse than for the 2G-HDNNP, as the unphysical charge distribution to some extent prevents the accurate representation of the energy.

Table 1

Energy and charge error obtained for the organic molecules. Energy error (meV/atom) and mean errors of the atomic charges (10⁻³ e) and forces (eV/Å) of C₁₀H₂ and C₁₀H$$ {}_3^{+}$$ with respect to DFT obtained with the different HDNNP generations for the DFT-optimized structures. For the 3G-HDNNP the results for scaled and unscaled charges are given.

		Energy	Charges	Forces
	2G-HDNNP	0.684	—	0.095
C10H2	3G-HDNNP (unscaled)	1.255	19.72	0.430
	3G-HDNNP (scaled)	2.193	10.76	0.138
	4G-HDNNP	0.463	4.820	0.032
	2G-HDNNP	0.922	—	0.127
C10H$$ {}_3^{+}$$	3G-HDNNP (unscaled)	0.046	17.82	0.658
	3G-HDNNP (scaled)	1.425	17.72	0.259
	4G-HDNNP	0.176	5.048	0.042

To investigate the forces in more detail, in Fig. 4 we plot the individual atomic forces in both molecules using the 2G-HDNNP and the 4G-HDNNP for the DFT-optimized structures. For all atoms in both molecules the 4G-HDNNP yields very low-force errors, with an average error of only 0.037 eV/Å underlining the quality of this PES. However, for the 2G-HDNNP the forces acting on the left half of C₁₀H$$ {}_3^{+}$$ and on all atoms in C₁₀H₂ the force errors are significantly larger. The reason is again the 2G-HDNNP cannot distinguish both molecules for these atoms, and the force errors are only low close to the extra proton in C₁₀H$$ {}_3^{+}$$, which can be recognized as a distinct local structural feature in the atomic environments of the right half of this molecule.

Fig. 4

Force errors of the HDNNPs for the organic molecules.

2G- and 4G-HDNNP forces for the atoms in the DFT-optimized structures of C₁₀H₂ and C₁₀H$$ {}_3^{+}$$ (indicated in a and b, respectively).

Interestingly, the relatively high errors of the 2G-HDNNP forces are not matched by high energy errors, which instead are surprisingly low and smaller than 1 meV/atom for both molecules. This suggests that the total energy predicted by 2G-HDNNPs may benefit from error compensation in the atomic energies in that the atomic energies in the right half of C₁₀H$$ {}_3^{+}$$ are adjusted to compensate the deficiencies of the atomic energies in the left half of the molecule.

Metal clusters: Ag₃

In this example, we investigate a small metal cluster, Ag₃, in two different charge states. The potential energy surface of small clusters is strongly influenced by the ionization state of the cluster and the ground state can differ as a function of the total charge of the cluster^48–51. Owing to the small system size there are no long-range effects, and the full system is included in each atomic environment. Therefore, in principle 2G-HDNNPs should be perfectly suited to describe the PES of Ag₃, but this is only true as long as the total charge of the system does not change, since for a combination of data with different total charges, like Ag$$ {}_3^{+}$$ and Ag$$ {}_3^{-}$$, in the training set the unique relation between atomic positions and the energy is lost. The minimum-energy structures of both cluster ions obtained from DFT are shown in Fig. 5a along with the atomic partial charges. After training a 2G-HDNNP and a 4G-HDNNP to data containing both types of clusters, we have reoptimized the geometries by the respective HDNNP generation. As expected, the minima obtained with the 2G-HDNNP (Fig. 5b) are identical for both charge states, but do not agree with any of the DFT structures. The 4G-HDNNP on the other hand, which in addition to the structural information also takes the total charge and the resulting partial charges into account, is able to predict the minima and also the atomic partial charges of both systems with very high accuracy (Fig. 5c). In this case, the inability of the 2G-HDNNP to distinguish between clusters is also apparent from the energy errors with respect to DFT. While the energy errors for Ag$$ {}_3^{-}$$ and Ag$$ {}_3^{+}$$ obtained from the 4G-HDNNP are only about 1.166 meV/atom and 0.320 meV/atom, respectively, the errors of the 2G-HDNNP are 0.605 and 2.017 eV/atom and thus several orders of magnitude larger. The 3G-HDNNP using scaled charges performs even worse and errors of 0.713 and 5.721 eV/atom are obtained. This is due to the non-physical electrostatic contribution calculated from the incorrectly predicted charges.

Fig. 5

Optimized geometry and atomic charges of Ag clusters.

Structures and atomic partial charges of Ag$$ {}_3^{+}$$ and Ag$$ {}_3^{-}$$ optimized with DFT in a, the 2G-HDNNP in b and the 4G-HDNNP in c. The numbers give the root mean squared displacement (RMSD) in Å compared to the respective DFT minima. The partial charges in b are shown for illustration purposes only and have been obtained from a scaled 3G-HDNNP.

NaCl cluster ions

As the last non-periodic example we select a system with mainly ionic bonding, which is a positively charged Na₉Cl$$ {}_8^{+}$$ cluster, and we analyze the changes of the PES, if a neutral sodium atom is removed. The initial structure of the cluster ion has been obtained from a DFT geometry optimization and is shown in Fig. 6. The sodium atoms are shown in purple, blue, and brown, while the chlorine atoms are displayed in gray. We then construct a second system by removing the brown sodium atom from the cluster while keeping the positions of the remaining atoms fixed. Since the overall positive charge of the cluster is maintained, the charge is redistributed throughout the new Na₈Cl$$ {}_8^{+}$$ cluster ion.

Fig. 6

Optimized structure of the Na₉Cl$$ {}_8^{+}$$ cluster.

Sodium atoms are shown in purple, blue and brown, chlorine atoms in gray. The arrow indicates the direction along which the blue sodium atom is moved for the energy and force plots in Fig. 7a and 7b. The position of this atom is defined by the Na–Na distance indicated as dashed line.

To investigate the consequences of this change in the electronic structure on the PES, we compute and compare the energies and forces when moving the blue sodium atom along a one-dimensional path indicated by the arrow in Fig. 6 for both cluster ions. The distance to the closest neighboring sodium atom highlighted as dashed line is used to define the structure.

Figure 7 shows the energies for both systems obtained with DFT, as well as the 2G-, 3G- and 4G-HDNNPs. All energies are given as relative energies to the minimum DFT energy of the respective cluster ion and refer to the full systems. First, we note that the positions of the DFT minima differ by more than 0.1 Å, i.e., depending on the presence of the very distant brown atom the blue atom adopts different equilibrium positions. The 2G-HDNNP, however, is unable to distinguish these minima and instead the same local minimum Na–Na distance is found for both systems, which is approximately the average value of the two DFT minima. We note that the 2G-HDNNP energy curves of the two systems are not identical but there is an energy offset, as some of the atomic environments in the right part of the systems differ yielding different atomic energies. Since these environments do not change when moving the blue atom this offset is constant. For the 3G-HDNNP the same qualitative behavior is observed, and two very similar but not identical minima are found for both systems. Still, in case of the 3G-HDNNP the energy offset between both systems is not merely a constant anymore, as the long-range electrostatic interactions between the blue and the brown atom in Na₉Cl$$ {}_8^{+}$$ are position-dependent. We note that in spite of these qualitative differences with respect to DFT, the 2G- and 3G-HDNNP curves show only a deviation of about 1 meV per atom from the DFT curves. This is very small and in the typical order of magnitude of state-of-the-art ML potentials, and in the present case this apparently high accuracy hides the qualitatively wrong minima. Finally, the 4G-HDNNP energies for both systems are very accurate and the energy curves match the corresponding DFT curves very closely. Both distinct local minima are correctly identified and at the right positions.

Fig. 7

Relative energies and forces of the NaCl clusters.

a Relative energies of all potentials with respect to the DFT minima of the Na₈Cl$$ {}_8^{+}$$ and the Na₉Cl$$ {}_8^{+}$$ clusters as a function of the Na–Na distance and b forces acting on the blue sodium atom for the the path shown in Fig. 6. For the 3G-HDNNP unscaled charges have been used in this plot.

Next, we turn to the forces shown in Fig. 7b. The results are fully consistent with our discussion of the energy curves. The DFT forces acting on the displaced atom are different for both cluster ions and well reproduced by the 4G-HDNNP. The 2G-HDNNP forces of both systems are exactly identical due to the constant offset between both energy curves (Fig. 7a), while the 3G-HDNNP forces of both systems are slightly different due to the additionally included long-range electrostatics.

Au₂ cluster on MgO(001)

As example for a periodic system we choose a diatomic gold cluster supported on the MgO(001) surface. Similar systems have attracted attention because of their catalytic properties for reactions like carbon monoxide oxidation, epoxidation of propylene, water-gas-shift reactions, and the hydrogenation of unsaturated hydrocarbons⁵². Theoretical^53,54 as well as experimental studies⁵⁵ have shown that the geometry of these clusters can be modified by the introduction of dopant atoms into the oxide substrate. This ability to control the cluster morphology is of great interest, as it can enhance the catalytic activity of the system⁵⁴. 2G-HDNNPs have been used before to study the properties of supported metal clusters^56–58, but systems as complex as doped substrates to date have remained inaccessible, since long-range charge transfer between the dopant and the gold atoms is crucial to achieve a physically correct description of these systems.

For Au₂ at MgO(001) there are two main adsorption geometries, an upright “non-wetting” orientation of the dimer attached to a surface oxygen and parallel to the surface in a “wetting” configuration, in which the two Au atoms reside on two Mg atoms. DFT optimizations of the positions of the gold atoms with fixed substrate for the doped and undoped surfaces reveal that the presence of the dopant atoms changes the relative stability of both structures. On the pure MgO support (Fig. 8a) the minimum-energy structure is “non-wetting”, while a flat “wetting” geometry is more stable if the MgO is doped by three aluminum atoms (Fig. 8b) corresponding to 2.86% of the slab. The Al dopant atoms were introduced into the 5th layer, resulting in a distance of >10 Å from the gold atoms. Despite this large separation, we found that by doping the charge on the Au₂ cluster is reduced (becomes more negative) by about 0.2 e compared to the same geometry for the undoped surface. This change in the electronic structure does not only lead to a switching in the energetic order of the geometries but also to a change of the bond-length between the gold atoms and the substrate.

Fig. 8

Geometry of Au₂ clusters on undoped and doped MgO(001) surface.

Au₂ cluster in the non-wetting geometry on the undoped a and the wetting geometry on Al-doped b MgO(001) surface represented by a periodic (3 × 3) supercell. Au atoms are shown in yellow, O in red, Mg in green and Al in blue. The configuration of the gold cluster has been optimized by DFT for a fixed substrate. The structure visualization for periodic systems was carried out using VESTA⁶⁷.

The energy difference (E_wetting − E_non-wetting) between the wetting and non-wetting configurations calculated with different methods on a doped substrate are −2.7 meV for DFT, 375 meV for the 2G-HDNNP and −41 meV for the 4G-HDNNP. On an undoped substrate we obtained 929 meV for DFT, 375 meV for the 2G-HDNNP and 975 meV for the 4G-HDNNP. These numbers were obtained after the positions of the gold atoms were optimized. In case of the 2G-HDNNP, both optimizations yield the same structure. For the 2G-HDNNP the energy differences for the doped and undoped systems are exactly the same as the dopant atoms are outside the local chemical environments of the gold atoms. Thus, the 2G-HDNNP cannot take the change of the PES by doping into account. The DFT and 4G-HDNNP results agree in that there is a slight preference for the wetting configuration for the doped surface, while in the undoped case the non-wetting configuration is clearly more stable.

An analysis of the PES for the case of the non-wetting geometry for the doped and undoped slabs is given in Fig. 9, which shows the energies relative to the minimum DFT energies of the respective systems as a function of the distance between the bottom Au atom and its neighboring oxygen atom for DFT, the 2G-HDNNP and the 4G-HDNNP. The energy curves of the 4G-HDNNP and DFT are very similar and can resolve the different equilibrium bond lengths for the doped (4G-HDNNP: 2.342 Å; DFT: 2.332 Å) and undoped (4G-HDNNP: 2.177 Å; DFT: 2.190 Å) substrates. The 2G-HDNNP yields the same adsorption geometry with a bond-length of 2.256 Å in both cases, while the energies substantially differ from the DFT values with the main effect of the dopant being a constant energy shift between both substrates, similar to what we have observed in the presence or absence of the additional sodium atom in the NaCl cluster.

Fig. 9

Energies and forces for the gold cluster.

a Relative energy and b sum of forces acting on the Au₂ cluster for the cluster adsorbed at the MgO(001) substrate for the non-wetting geometry for the Al-doped and undoped cases. The local minima of the energy curves are marked with a dot. The Au–O bond-length refers to the distance between the Au closest to the surface and its neighboring oxygen atom.

Neural network potentials

The HDNNPs reported in this work have been constructed using the program RuNNer^59–61. Atom-centered symmetry functions⁴¹ have been used for the description of the atomic environments within a spatial cutoff radius set to 8–10 Bohr depending on the system. For a given system, the same parameters of the symmetry functions and the same atomic neural network architectures have been used for the different generations of HDNNPs being compared, and the parameters and cutoff radii for all systems can be found in supplementary tables. The functional forms of the symmetry functions are given in ref. ⁴¹. In all examples, the atomic neural networks consist of an input layer with the number of symmetry functions ranging from 12 to 54 depending on the specific element and system, two hidden layers with 15 neurons each, and an output layer with one neuron providing either the atomic short-range energy or electronegativity. Forces have been obtained as analytic energy derivative. The activation functions in the hidden layers and the output layer were the hyperbolic tangent and the linear function, respectively.

In all cases 90% of the available reference data was used for training the HDNNPs while the remaining 10% of the data points were used as an independent test set to confirm the reliability of PESs and detect possible over-fitting. Energies and forces were used for training the short-range atomic neural networks.

Moreover, a screening of the short-range Coulomb electrostatic interaction was applied in order to facilitate the fitting of the short-range energies and forces obtained from Eq. (9)²³. The inner and outer cutoff radius for screening of the electrostatic interaction have been set to 1.69–2.54 Å and the cutoff of the symmetry functions, respectively. The widths of the Gaussian charge densities in Eq. (4) have been set to the covalent radii of the elements. All the details of the training process and the validation strategies for HDNNPs in general can be found in recent reviews^60,61.

The HDNNP-based geometry optimizations were performed using simple gradient descent algorithms and the numerical threshold of the forces was set to 10⁻⁴ Ha/Bohr ≈ 0.005 eV/Å, which is the same convergence used in the DFT calculations used for validating the HDNNP results.

DFT calculations

The DFT reference data has been generated using the all-electron code FHI-aims⁶² employing the Perdew–Burke–Ernzerhof⁶³ (PBE) exchange-correlation functional with light setting. The total energy, sum of eigenvalues, and charge density for all systems except Au₂-MgO were converged to 10⁻⁵ eV, 10⁻² eV, and 10⁻⁴ e, respectively. For the Au₂-MgO systems stricter settings have been applied by multiplying each criterion by a factor 0.1 in combination with a 3 × 3 × 1 k-point grid. Spin polarized calculations have been carried out for the Au₂-MgO, NaCl and Ag₃ systems. Reference atomic charges were calculated using Hirshfeld population analysis⁴⁰. In principle any other charge partitioning scheme could be used in the same way.

The data set of the C₁₀H₂/C₁₀H$$ {}_3^{+}$$ molecules and the Ag₃ clusters have been constructed by performing Born-Oppenheimer molecular dynamics⁶⁴ simulations for each system at 300 K with 5000 steps at a time step of 0.5 fs. A Nosé-Hoover thermostat⁶⁵ was applied to run simulations in the canonical (NVT) ensemble, and the effective mass was set to 1700 cm⁻¹. In addition, the trajectory path during the geometry relaxations up to a numerical convergence of 0.001 eV/Å of the forces was also added to the data set to have sufficient sampling close to equilibrium structures. The geometry optimization of the Ag$$ {}_3^{-}$$ system has been terminated when reaching forces below 0.0015 eV/Å.

In case of the NaCl cluster and the Au₂ cluster at the MgO surface the reference data set consists of two structurally different types of systems, and half of the data set was dedicated to each of the two cases. We performed a random sampling along the trajectories depicted in Figs. 7 and 9 and added further Gaussian distributed displacements to ensure sufficient sampling of the PES in the vicinity of the structures of interest. For the NaCl cluster we used Gaussian displacements with a standard deviation of 0.05 Å. As in the Au₂-MgO system we only investigated the change in geometry of the Au₂ cluster, while the MgO substrate remained fixed during all geometry relaxations, we used a smaller magnitude of the Gaussian displacements for the substrate than for the cluster. A standard deviation of 0.02 Å was used for the substrate and 0.1 Å was used for the gold cluster. Half of the data set consists of structures with an undoped substrate, while the other half includes a doped substrate. Half of the samples of each substrate configuration were generated with the Au₂ cluster in its wetting configuration, and the other half with the cluster in its non-wetting configuration. The total number of reference data points for the NaCl cluster and Au₂-MgO slab is 5000, while the the Ag₃ clusters and the organic molecule it is 10,019 and 11,013, respectively.