Machine learning (ML) potentials and force fields [1][2][3] have revolutionized atomistic modeling by facilitating larger and longer simulations crucial for modeling dynamic and kinetic properties. 4,5 General-purpose ML potentials (e.g., ANI-2x 6 , OrbNet Denali, 7 AIQM1 8 ) model chemical (local) interactions and can be useful for subsets of chemical space. These approaches assist molecular screening but require enormous data sets of hundreds of thousands of structures. Alternatively, ML potentials can be tailored to specific systems to improve desired simulation reliability. This requires that models be retrained for each system, so training must involve minimal human involvement and computations to be practical.
Size transferability to hundreds of molecules is paramount for useful ML potentials. Most ML potentials rely on local descrip-tors [9][10][11][12] or graph neural networks (GNNs) 13,14 that partition total properties into atomic contributions. Local descriptors have been successful in numerous applications, but they inherently neglect or limit complicated non-local interactions by enforcing atomic radial cutoffs. For example, a recent study showed that a deep neural network potential's predictions of liquid water properties are sensitive to training data relevant to the thermodynamic state point. 15 Global descriptors (such as the Coulomb matrix and pairwise atomic distances) impose no such constraints and capture interactions at all scales. 16,17 Still, they are usually restricted to the same number of atoms.
Gradient-domain ML (GDML) uses a global descriptor and has demonstrated remarkable success in many chemical applications with monomers or dimers. [18][19][20][21][22] Moreover, GDML only needs energies and forces of approximately 1000 structures to accurately learn the potential energy surfaces of molecules 19 and periodic materials. 17 The global descriptor limits GDML to the same system it was trained on, whether a single molecule or a chemical reaction. Size-transferable GDML for molecular ensembles would provide rapidly trained force fields for high-quality molecular simulations involving solvents.
Many-body expansions (MBEs) should enable size-transferable GDML because systems with non-covalent clusters are naturally described in terms of n-body interactions. 23,24 Data-driven, many-body potentials (e.g., MB-pol) have already been widely successful in modeling aqueous systems. [25][26][27] This expansion is formally exact if all N-body interactions are accounted for with sufficient accuracy and precision. However, the expansion is typically truncated to the third order due to combinatorics. One can avoid truncating the expansion and include all contributions by using a classical many-body polarization model (e.g., a Tholetype model as in MB-pol 25 ). We expect training on fundamental n-body interactions found in clusters would extend GDML force fields to be useful for bulk liquid simulations. Alternative approaches exist; for example, Gaussian Approximation Potential (GAP) 3,11 was extended to liquid methane by decomposing energies into fundamental interactions (e.g., repulsion, dispersion, and electrostatic contributions) and different scales. 28 This is another rigorous approach that requires considerable effort with large numbers of quantum chemical calculations.
MBEs share characteristics with local descriptors but provide several key advantages. First, n-body interactions are more efficiently treated on a molecular basis instead of an atomic basis. Second, errors associated with MBE truncation can be corrected using a variety of schemes. For example, using long-range physical models to capture induction and dispersion effects. 29,30 Alternatively, one could use ONIOM-style 31 approaches such as molecule-in-molecules (MIM) 32 and molecular tailoring approach (MTA) 33 where low-cost calculations on the whole structure (i.e., supersystem) are used to capture all long-range interactions. Third, these n-body contributions can be observed in relatively small clusters. Local descriptors require data on large clusters to achieve similar levels of size transferability. 34 This opens the door for many-body GDML (mbGDML) force fields trained on high levels of theory that scale poorly with system size, such as CCSD(T). In addition, mbGDML naturally incorporates intramolecular/monomer flexibility, which is extremely challenging for analytical potentials.
Thus, mbGDML should provide size-transferable force fields trained on highly accurate quantum chemical methods. To evaluate this, we developed an automated framework in Python to facilitate training and application of mbGDML force fields (available at github.com/keithgroup/mbGDML). GDML force fields for water (H 2 O), acetonitrile (MeCN), and methanol (MeOH) were trained on 1000 structures for 1-, 2-, and 3-body interactions. GAP and SchNet 35 were also evaluated in this many-body framework. The size transferability of mbGDML was further assessed against a highly promising graph neural network, Neural Equivariant Interatomic Potentials (NequIP). 13 Reference structures from the literature were used to benchmark energy and forces predictions. The following sections demonstrate mbGDML energy and force accuracies within 0.38 kcal/mol per monomer and 0.06 kcal/(mol Å) per atom for structures containing up to 20 monomers (120 atoms). The MBE framework itself contributes 14% to 83% of these errors depending on the system. Error cancellation dramatically improves relative energy predictions of mbGDML to less than 3 kcal/mol and achieves fair to excellent agreement with solvent radial distribution functions (RDFs).
The MBE represents the total system energy, E, composed of N noncovalently connected (i.e., non-intersecting) fragments as the sum of n-body interaction energies: 36
Here, N is the number of monomers; i, j, k are monomer indices; E
(1) i is the energy of monomer i; and DE (n) i, j, ... represents the nbody interaction energy contribution of the fragment containing monomers i, j, ... with lower order (< n) contributions removed. For example, the 2-body contribution of the fragment containing monomers i and j is
and the 3-body contribution with monomers i, j, and k are
Equation 1 is exact when all n-body contributions up to N are accounted for with exact accuracy and precision. This equation also holds for properties expressed as a derivative of energy (i.e., gradients).
The xTB program 37 v6.4.0 was used to run MD simulations of the three solvents at 500 K. Small clusters containing up to three molecules were sampled from these simulations to generate data sets for training. Higher temperatures provided configurations relevant at lower temperatures with the added benefit of sampling high-energy regions. 2 GFN2-xTB, 38 a semiempirical quantum mechanics method, was used as a compromise between the cost of quantum chemical methods and potentially not having classical force field parameters for species of interest. Furthermore, simulation accuracy is not a significant concern because only reasonable geometries are desired at this stage.
Equation 1 represents the MBE framework where individual GDML force fields are trained on intramolecular (i.e., 1-body) and intermolecular (i.e., 2-and 3-body) energies and forces. Energies and forces were calculated with ORCA v4.2.0 39,40 using secondorder Møller-Plesset perturbation (MP2) theory, 41 the def2-TZVP basis set, 42 and the frozen core approximation. This level of theory was chosen for its efficiency and accuracy for noncovalent interactions, but future applications of mbGDML are recommended to use the highest levels of quantum chemical theory available for training data. The Resolution of Identity (RI) approximation was only used for calculations containing 16 or more monomers. Additional calculation details and discussion can be found in the Electronic Supplementary Information (ESI).
We evaluated mbGDML, mbGAP, and mbSchNet on tetramers (4mers), pentamers (5mers), and hexamers (6mers) from the literature. [43][44][45] These test structures have minimal higher-order (> 3-body) contributions that increase with the number of monomers. Furthermore, many-body ML (mbML) potentials considered here implement a distance-based cutoff for 2-and 3-body contributions (see the ESI for more details). Small clusters allow us to determine whether errors are from the underlying MBE framework or ML predictions.
ML potentials discussed here are trained on small data sets of only 1000 structures to showcase GDML data efficiency. Training sets were determined through an iterative training procedure to reduce the maximum model error. 46 GAP and SchNet models were trained on the same training sets as GDML for a fair comparison. In theory, training sets could be tailored for GAP and SchNet to reduce errors; however, a cursory attempt did not substantially improve results. We reiterate that GAP and SchNet normally require substantially large training sets. In other words, GAP and SchNet potentials presented here are technically underfitted compared to standard practices. More information can be found in the ESI.
Fig. 1 shows relative isomer energies with respect to the lowest energy structure for MBE (light color) and mbGDML (dark color) methods. The ESI provides comparable figures for mb-GAP and mbSchNet. Figures showing absolute energy predictions for these structures are also shown in the ESI and help determine where error cancellation comes into play. First, we discuss the inherent errors in MBE data versus supersystem MP2 data (gray). These MBE predictions generally capture the relative energy trends of water, acetonitrile, and methanol isomers. Water predictions showed increasing errors with system size, indicating the importance of higher-order contributions (as expected). Acetonitrile 5mers and 6mers (Fig. 1D-F) show small energy differences that are not monotonically increasing. This is likely due to challenging electrostatics and polarization from dipole-dipole interactions. Methanol isomer MBE predictions showed this same trend as water, but higher-energy structures now exhibit lower MBE errors. This suggests that higher-order contributions are crucial for stabilizing low-energy methanol structures. Incomplete basis sets and basis set superposition errors (BSSE) are known to impact MBE accuracy. [47][48][49][50][51][52] The def2-TZVP basis set was chosen for its balance of cost and accuracy, as the larger aug-cc-pVTZ basis set only improved the energy MAE by 0.15 kcal/mol. BSSE corrections are not included here because the n-body energies and forces would depend on the original supersystem-thereby limiting data set transferability.
We now discuss mbGDML data, which approximates the MBE potential energy surface. In general, mbGDML reasonably mimics MBE data, including innate errors made by the MBE framework, as seen in the acetonitrile 5mer and 6mer data. Note that fortuitous error cancellations of 2-and 3-body mbGDML predictions sometimes give the appearance of higher accuracy than MBE. mb-GAP and mbSchNet potentials occasionally are better or worse than mbGDML; however, as previously mentioned, these methods generally require larger training sets and are likely to underperform. For example, Table 1 shows the MBE, mbGDML, mbGAP, and mbSchNet energy and force mean absolute errors (MAEs) with respect to supersystem MP2/def2-TZVP calculations for all 4-6mer structures considered here. All mbML models perform similarly for water and acetonitrile, but the methanol isomer errors 53 Their 2-body training set included 34 431 structures containing the global dimer minimum, saddle points, artificially compressed geometries, and geometries from path-integral molecular dynamics (PIMD) simulations using HBB2-pol. 54 Their 3-body training set contained 10 001 structures from HBB2-pol MD and PIMD of small water clusters, liquid water, and ice phases. Table 2 shows 2-and 3-body interaction energy MAEs with their models and those calculated here. The PIP, BPNN, and GAP water models trained on large data sets achieved 2-and 3-body interaction energy MAEs on the order of 0.033-0.145 and 0.007-0.123 kcal/mol, respectively. Alternatively, our GDML force fields trained on only 1000 structures achieved MAEs of 0.030-0.047 and 0.041-0.093 kcal/mol for 2and 3-body interaction energies. This shows that GDML models using small training sets can perform similarly to well-trained potentials requiring large training sets. 53 We highlight that the GAP results from Ref. 53 demonstrate substantial accuracy improvements possible with more extensive training sets. We reiterate that ML potential accuracy is intricately linked to reference data sets, but we specifically opted to show the promise of GDML with small training sets. In almost all cases, the water 2body models prepared here outperformed those from Ref. 53 that used larger training sets. Presumably, our smaller data sets may contain structures that enhance the perceived accuracy of these models. The 3-body data exhibit the opposite trend, which can be attributed to data set quality. In general, additional sampling of configurational spaces would improve our mbML models; however, the objective here is to evaluate ML potentials that can be trained with minimal computational cost.
Predictions of medium-sized structures provide a straightforward test of size transferability. There are additional, albeit typically small, higher-order contributions in larger structures. Also, the n-body cutoffs are now in effect to reduce the number of computations. Table 3 shows energies and forces of hexadecamers (16mers) from the literature [55][56][57] computed with RI-MP2/def2-TZVP and compared against mbGDML, mbGAP, and mbSchNet results. The truncated MBE contributes a few kcal/mol errors depending on the system. For example, the MBE prediction of (H 2 O) 16 results in a 3.3 kcal/mol error whereas (MeCN) 16 has only a 0.2 kcal/mol error. Missing higher-order contributions or basis set errors are the most likely causes. All mbML models performed similarly well with (H 2 O) 16 . Most errors originated from
We find that both mbSchNet and mbGAP models trained from smaller data sets have abnormally high errors for (MeCN) 16 and (MeOH) 16 , respectively. In both cases, the 3-body model has substantial error accumulation. Cutoffs are not the issue because only 0.006 of the 16.1 kcal/mol error in mbSchNet's prediction of (MeCN) 16 is from cutoffs implemented in the 2-and 3-body models. Prediction errors contribute the most; a massive 15.2 kcal/mol error comes from the 3-body SchNet model.
Assessing inadequacies of training data is more complicated. If 3-body structures from (MeCN) 16 are substantially different from the data sets, then the model may break down. To investigate this, we used dimensionality reduction to visualize high-dimensional similarity in 2D space. Similar structures in feature space should be clustered together and vice versa. Fig. 2A shows the GDML feature space, a 2D embedding of trained and 3-body structures from (MeCN) 16 using Uniform Manifold Approximation and Projection (UMAP). 58 There is a significant overlap between the GDML training set feature space and the structures from (MeCN) 16 . High overlap suggests that GDML should have low prediction errors, which is the case. SchNet, on the other hand, has several test structures isolated from training data, resulting in higher errors (shown in the ESI). Not all structures with a high error are dissimilar in feature space. Models should have learned similar structures and thus should have performed well. A simple, ad hoc geometry descriptor (discussed in the ESI) applied in Fig. 2B shows that all higherror structures are dissimilar to anything in the data set. SchNet has some difficulty with these structures, which results in a substantial 16.1 kcal/mol error. Many-body GAP's 17.6 kcal/mol error in (MeOH) 16 is likely for the same reason. However, GAP uses a local descriptor, making feature space more complicated to analyze. Models under these circumstances were excluded from further analyses (namely, mbSchNet for acetonitrile and mbGAP for methanol).
Truncated higher-order contributions could be pertinent for accurate absolute energies, as seen in the previous 16mer data. In practice, relative energy accuracy is of primary importance. Yao et al. 59 trained mbML methanol potentials and analyzed their performances on five (MeOH) 20 isomers. They used a Generative Adversarial Network (GAN) trained on RI-MP2/cc-pVTZ energies with the Coulomb matrix descriptor. Training included 80% of their data sets that contained 844 800 monomers, 74 240 dimers, and 36 864 trimers.
Relative isomer energies of their methods are reported in Table 4. Their mbGAN potential accurately captures the isomer rank-
As previously mentioned, many ML potentials use local descriptors for size transferability. Recent developments of ML potentials with local descriptors have involved GNNs. 13,14 NequIP uses equivariant, continuous convolutions where edges connect every atom within a cutoff radius. 13 NequIP has achieved remarkable accuracy and data efficiency on the MD17 data set, bulk water, formate dehydrogenation, and amorphous lithium phosphate. Such models are inherently size transferable, but the accuracy is not typically studied when trained exclusively on small clusters. In theory, these potentials can train on the same data sets, but instead of n-body interactions, they would use total energies and forces. This would eliminate the need for an MBE framework. We trained NequIP on total energies and forces of 1000 trimers for water, acetonitrile, and methanol to assess this approach. Another 2000 trimers were used as a validation set.
We emphasize that this is an edge case of GNN potentials. If energies and forces of larger structures were readily available, these data would improve size transferability if they were included in the validation set. However, mbGDML models were never exposed to these larger clusters during training since the objective was to reproduce the 1-, 2-, and 3-body PES. Thus, training a NequIP on only trimers represents a straightforward comparison to mbGDML. These models were then tested against the identical isomers discussed above, with the results shown in Table 5. While NequIP can expectedly extrapolate to larger clusters, the Training on large clusters or bulk systems is likely more efficient if a lower scaling method is satisfactory. However, mbGDML becomes particularly useful when applications require force fields based on higher scaling methods. Recovering truncated higherorder contributions would also expectedly improve errors, but explicit 4-body interactions are rather challenging due to high demands on precision and combinatorics. 48,49,60 Electrostatic 61 and more general quantum embedding approximations may be a practical route to avoid calculating higher-order contributions, but they are not considered here.
While accurate predictions of static clusters are essential, compelling applications for mbGDML would involve molecular simulations. Low energy and force errors are not conclusive of accurate molecular simulations, 62 but experimentally measurable dynamic properties are an alternative and rigorous way to evaluate ML potentials. For example, the radial distribution function (RDF) is a vital bulk property that quantitatively defines liquid structure. Locations and intensities of peaks and valleys represent the solvation shells and liquid ordering. Accurately reproducing reference RDF curves is crucial for a practical size-transferable potential.
Periodic NVT simulations driven by mbGDML force fields were performed at 298.15 K for 10-30 ps in the atomic simulation environment (ASE). 63 Note that NVT simulations could artificially bias intermolecular distances due to the volume constraint. NPT simulations would be a more rigorous metric, but these are not yet implemented in mbGDML, and this will be the focus of future work. The minimum-image convention was used with cubic boxes with lengths of 16 Å (137 molecules), 18 Å (67 molecules), and 16 Å (61 molecules) for water, acetonitrile, and methanol, respectively. Production trajectories were used to compute all possible RDFs. Some RDF curves are shown in Fig. 3. Wa- ter, 64 acetonitrile, 65 and methanol 66,67 reference RDF curves are from neutron diffraction experiments. Results from classical MD simulations [68][69][70][71][72][73][74] are also shown in Fig. 3. Note that classical references often include some fitting to empirical data, 68,[71][72][73][74] whereas mbGDML and others 69,70 run calculations with no explicit empirical fitting. Individual figures of all RDF curves, along with labeled classical references, are shown in the ESI.
Dispersion and polarization are not always accurately treated with MP2 theory, 75,76 and likewise, the underlying model chemistry (MP2/def2-TZVP) to train mbGDML force fields may not accurately reproduce experimental liquid properties. For example, MP2 yields excellent results for liquid water simulations when appropriate density corrections or basis set error cancellation schemes are employed, 77,78 but these were not used here. To our knowledge, a thorough investigation has not been performed for liquid acetonitrile and methanol with MP2, so the agreement with experiments is more uncertain. Note that molecular simulations using Kohn-Sham density functional theory (DFT) results in comparable differences in RDFs shown in Fig. 3, depending on the exchange-correlation functional and dispersion treatment used. [79][80][81][82][83] The simulated RDFs with mbGDML fairly agree with the reference curves. In particular, the water g OO (r) in Fig. 3A agrees remarkably well with experimental data. This is consistent with fragment-based ab initio MD (AIMD) simulations. 84,85 However, these AIMD simulations include higher-order contributions through electrostatic embedding. Deviations in the g OH (r) and g HH (r) curves are partially due to the neglect of quantum nuclear effects. 86 In all cases, acetonitrile peaks from mbGDML are less intense than the reference curves. This indicates that the predicted liquid structure with mbGDML is less ordered than the deuterated neutron diffraction data. 65 Notably, g NN (r) is wide with two distinct peaks that deviate from the experimental reference. However, classical RDFs from the literature can vary substantially. Some classical potentials 70,73 result in a similar g NN (r) shape while others 69,71,72 better resemble the experimental reference.
Methanol simulations appear more challenging for mbGDML. RDF peaks with respect to experimental data are less intense (same as acetonitrile). The shape of g OO (r), Fig. 3C, agrees well with the digitized experimental data. Classical simulations using GROMOS96 and OPLS/AA potentials have significantly more ordered liquid structure. 74 For instance, their g OH (r) peaks are around 1.24 higher in intensity than mbGDML. While the g OO (r) is in good agreement with the experiment beyond 5 Å, the g OH (r) and g HH (r) curves are missing long-range liquid structure. Even though GDML employs a global descriptor, mbGDML is not capturing these long-range interactions. We suspect this is caused by truncations and cutoffs used in the MBE framework.
To summarize, even though the mbGDML models used here only included up to 3-body contributions, they generally predict the liquid structure of water, acetonitrile, and methanol well. Moreover, these force fields automatically include fully flexible molecules and perform no fitting to experimental properties. Further improvements could be made with more expansive training sets and higher-order contributions. For systems without classical parameters, mbGDML can be rapidly trained on relatively small amounts of data and provide valuable dynamical insights for explicitly solvated systems.
We have introduced a GDML-driven, many-body expansion framework that enables state-of-the-art size transferability toward molecular simulations of solvents. mbGDML force fields trained on only 1000 1-, 2-, and 3-body interactions accurately modeled small and medium isomers of water, acetonitrile, and methanol. Size-extrapolated predictions on static clusters of up to 20 monomers had energy errors of less than 0.38 kcal/mol per monomer for all three solvents. These results outperform NequIP trained on the same trimer data set by up to 34 kcal/mol for 16mers. Dynamic simulations of bulk systems using our mbGDML force fields provide semi-quantitative insights while avoiding ex-pensive training data on bulk systems and fitting to experimental data.
It is important to note that the accuracy of mbGDML is generally limited to that of the underlying MBE framework. More extensive and diverse n-body data sets can help minimize mbGDML deviations from the MBE reference. If further accuracy improvements are desired, explicit 4-body ML force fields, classical models, or hybrid methods like MIM and MTA could be required. While these approaches are certainly possible to implement, we focused on providing a proof-of-concept of mbGDML. We thus anticipate promising applications for complex, explicitly solvated systems where high levels of theory are desired.
1-10
|1-10