{'fraction_non_alphanumeric': 0.033007019592402, 'fraction_numerical': 0.017460726941928826, 'mean_word_length': 4.873365231259968, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Regularized by Physics: Graph Neural Network Parametrized Potentials for the Description of Intermolecular Interactions

Simulations with an explicit description of intermolecular forces using electronic structure methods are still not feasible for many systems of interest. As a result, empirical methods such as force fields (FF) have become an established tool for the simulation of large and complex molecular systems. However, the parametrization of FF is time consuming and has traditionally been based largely on experimental data, which is scarce for many functional groups. Recent years have therefore seen increasing efforts to automatize FF parametrization and a move towards FF fitted against quantum-mechanical reference data. Here, we propose an alternative strategy to parametrize intermolecular interactions, which makes use of machine learning and gradient-descent based optimization while retaining a functional form founded in physics. This strategy can be viewed as generalization of existing FF parametrization methods. In the proposed approach, graph neural networks are used in conjunction with automatic differentiation to parametrize physically motivated models to potential-energy surfaces, enabling full automatization and broad applicability in chemical space. As a result, highly accurate FF models are obtained which retain the computational efficiency, interpretability and robustness of classical FF. To showcase the potential of the proposed method, both a fixed-charge model and a polarizable model are parametrized for intermolecular interactions and applied to a wide range of systems including dimer dissociation curves and condensed-phase systems.


Introduction
Computer based simulations are a powerful tool for the investigation of chemical systems [1]- [3]. Performing such simulations requires an accurate description of intermolecular forces [4]. However, due to the computational complexity of ab initio methods [5]- [9] or density functional theory (DFT) [10], [11], an exact description is out of reach for most systems, particularly systems of biological relevance [12], [13]. As a result, more approximate methods have been developed, which can be broadly categorized into three classes: semi-empirical methods, classical force fields (FF), and machine learning (ML) based models.
Semi-empirical methods explicitly describe the electronic structure [14], [15]. However, various approximations are introduced to reduce computational costs. Existing methods attempt to compensate for these approximations by introducing a relatively small number of empirical parameters [16]- [18]. Classical force fields, on the other hand, forego an explicit description of the electronic structure and employ instead a predefined functional form and associated parameters, which together encapsulate aspects of a given interaction [19], [20]. Due to their simplicity, they can be evaluated much more efficiently than ab initio or semi-empirical methods, but generally require a larger number of parameters and extensive parametrization. In recent years, ML-based models have emerged as a third alternative [21], which assume fewer inductive biases, but require an even larger number of parameters compared to the two previous approaches. Even though very promising results have been reported for ML potentials MLFF [21]- [26], their application to condensed-phase systems and the prediction of experimentally measured properties has been fairly limited [27]- [32].
The relative scarcity of application of ML potentials to propagate molecular dynamics (MD) simulations is likely a result of insufficient data efficiency due to a lack of inductive biases as well as difficulties posed by the presence of a large number of relatively weak and long-ranged interactions in condensed-phase systems [33]. However, there are also discussions whether ML potentials describe features of the potential-energy surface (PES), such as its curvature, sufficiently accurate to perform MD simulations [34]. With their predefined functional form, which is physically motivated, FF have become an established tool to simulate large condensed-phase systems such as solvated proteins over long time-scales [1]. While the FF formalism provides a computationally efficient, interpretable and robust way to describe forces in molecular systems, this robustness comes at a price. Most of the commonly used FF do not account for phenomena such as charge anisotropy or polarization. In addition, the development of a FF is still a non-trivial process, despite advances in automation over the past years [35]- [41]. The OpenFF initiative in particular has initiated a grand effort to fully automate this process, including atom-typing [42], [43], data generation [44], parametrization [45], and validation [46].
While FF were historically and continue to be (partly) fitted to experimental data, the advances in computational power and improved scaling of methods based on quantum-mechanical (QM) calculations has opened up new opportunities [19]. As a readily accessible data source, focus has shifted to parametrization with respect to QM reference data such as torsion profiles or interaction components. In addition, there is an increased effort to extract FF parameters directly from electron densities [47]- [53]. As an example, van Vleet et al. [54] developed an ab initio FF based on a parametrization formalism using Slater functions from which certain parameters were directly derived. An alternative approach is to keep a predefined functional form to describe intermolecular potentials but obtain the parameters from a ML model [55], [56]. In a similar fashion, Wang et al. [57] and Harris et al. [58] investigated graph neural networks (GNN) and graph-convolutional neural networks (GCNN) in combination with automatic differentiation as a method to parametrize FF. Focusing on intramolecular interactions, they could show that GNN can be used to predict FF parameters from potential energies and recover human-defined atom types. Finally, we also point out recent efforts which use ML in a complementary fashion to extract symbolic expressions from data [59]- [62].
In this work, we build on these developments to propose a universal framework for the parametrization of FF, focusing on intermolecular interactions. Besides generalization of the parametrization process, we describe a formalism for end-to-end differentiable FF, taking advantage of learned atom types. The proposed approach is applied to the parametrization of a non-polarizable as well as a polarizable FF. Multipoles and monopoles used to describe electrostatic interactions were obtained separately from our previously introduced equivariant GNN model without further modifications [63]. For both the non-polarizbale and polarizable FF, the model is trained on the PES of dimer interaction potentials of the recently published DES5M data set [64] and a data set of intermolecular potentials of molecular crystals, which was generated for this work. Hence, parameters are learned by the model from scratch to reproduce the given PES. We show that a fixed-charge FF parametrized in such a manner can be used to successfully reproduce experimental condensed-phase properties. We find that the resulting models outperform comparable models for a wide range of test cases.
The work is structured as follows: In the Theory section, the proposed formalism and potential energy terms are introduced. In the Methods section, the training and validation procedures is outlined. Finally in the Results and Discussion section, the performance of the models on a wide range of test systems and properties are discussed. MLFF … Figure 1: Overview of the proposed parametrization workflow: I) Atomic environments are encoded as feature vectors with a GNN. II) Parameters θ ij are predicted for tuples of atomic features. III) In conjunction with predefined parametric interaction potentials V , the predicted parameters θ ij are used to evaluate the potential energy of a state. IV) The prediction is evaluated against reference properties. V) Using automatic differentiation, errors with respect to the target properties are backpropagated to improve the quality of the predicted parameters.
θ to a molecular graph G and a functional form where the total potential energy of a state x is decomposed into V θ,i . The function G that assigns parameters θ can be understood as a parametrized, or learnable, function itself. For commonly used FF, G is generally expert-devised and depends only on atomic features such as element types and hybridization states. For FF that are parametrized based on electron densities, G partitions and assigns parameters to an electron density. Given a parametrization function G, a FF can thus be interpreted as a function V θ with parameters θ and a functional form V, V θ : X → V , which maps a PES V ∈ R to the states X ∈ R n of a system, with n denoting the dimensionality of the system. A system can thus be propagated in time by integrating the negative derivative of the potential energy v ∈ V with respect to its current state x ∈ X, where F ∈ R n is the gradient field of the negative potential energy. Propagating the system for an appropriate amount of time, a system property P can be derived as P : X → P with a function P, which assigns a property to a state or set of states of the system. Generally, system properties can be scalar, vectorial, or tensorial. They may be defined for each configuration or ensemble average, and depend on the functional form of the FF and its parametrization V θ (Figure 1). Given a system property obtained from a simulation, P pred and a reference value P ref (e.g., from an experiment or QM reference calculation), a loss L can be defined as L : (P pred , P ref ) → L with a loss function L. Due to its generality, any computable property can be used as a target. Examples include the potential energy, gradients or Hessians from QM reference calculations, but also experimental properties would be possible such as geometrical constraints from NMR or crystallography, vibrational spectra or ensemble properties such as enthalpies of phase transitions.

Parametrization
Given a loss function L and parameters θ, a FF may be optimized to yield more accurate properties. Specifically, the derivative of the loss function with respect to the FF parameters, can be used to update FF parameters successively with gradient descent based optimization methods analogously to backpropagation used in deep learning [65]. With the help of automatic differentiation libraries, these gradients can be evaluated with minimal additional programming effort and computational cost [66], [67].

Graph Neural Networks as Universal Parametrization Functions
GNN are ML models parametrized by artificial neural networks (ANN) that process graph-structured data. In the commonly used form, node, edge and/or global features are iteratively refined based on the current features. GNN models differ mainly by the features used, the way the underlying graph is constructed as well as the updating or feature-refinement process applied [68]- [71]. Considering a molecular graph G = (A, B) with nodes a i ∈ A and edges b ij ∈ B, message passing can be defined as with h l i ∈ R n describing the hidden feature vector of node a i after l graph processing steps, b ij ∈ R n the bond feature of bond b ij between node i and j, N (i) denoting the set of neighbours of a i and m ij the message from node j to node i. φ b and φ h represent ANN parametrized edge and node update functions. The superscript l denotes the current layer or the current iteration in the recurrent realization. After n iterations, the refined node feature h n i is used as an atomic-environment descriptor in subsequent steps. We note that h n i does not necessarily have to be obtained from a GNN. In principle, any other descriptor can be used, for instance atom-based topological fingerprints. However, GNN present themselves as a natural choice to work with graph-structured data, which in turn is ideally suited for a classical description of molecules.

Parameter Prediction
To predict atomic parameters, learned atom features h n i obtained from the GNN serve as descriptors of the atomic environment. Assuming that FF parameters are a function of the atomic environment, such features can be used to distinguish atom types and to assign FF parameters. Hence, atomic parameters are predicted as with φ atom denoting the readout function that assigns the parameter θ i for a given interaction to an atom type. In general, FF also describe how parameters of two particles are combined to describe a given interaction. Standard biomolecular FF use combination rules to derive parameters of two distinct atom types [19]. By design, combination rules must be symmetric under arbitrary permutations of types present in the interaction. Instead of using established combination rules such as arithmetic or geometric means, pairwise parameters are predicted as a function of two atom types. Specifically, the following approach was chosen In this case, φ pair is an ANN parametrized function, which assigns the parameters of a specific interaction to a pair of atomic environments. This approach guarantees symmetry with respect to permutation of two atoms and allows for more complex combination rules.

Force-Field Parametrization Function
Combining the previously introduced concepts with the formalism described in Section 3.1, the FF parametrization G φ can be defined as the combination of the following two components: (i) a typing function to assign atom types to a given system, and (ii) a combination rule, which returns parameters for a set of atoms partaking in a given interaction. In the current work, the function assigning atom types is modelled with a GNN consisting of a node and edge update layer, φ h and φ b , and a combination rule φ Atom or φ Pair for atomic and pairwise parameters, respectively. In the context of the GNN formalism, φ θ can be considered a readout function, which maps the hidden state of node features to a label. In our case, FF parameters are mapped to interactions between given atom types. Evidently, established FF atom type definitions and parametrization procedures can be viewed as a special case of the described formalism. For example, the SMARTS patterns used for chemical perception in Open FF [43] can be cast as graph-based operations that account for features such as the element, its coordination number, bonded neighbours, or sub-graph features. In the case of models, which derive parameters from electron densities, the partition function used to decompose the electron density into atomic contributions takes up the role of the typing functions. In these cases, the combination rules are often derived from first principle considerations or empirically fitted [47], [49].

Models
Besides the aforementioned components, which assign parameters to a given interaction or atom, a FF must further define potential-energy terms and an associated functional form. The chosen functional form reflects the assumptions of the model, and thus determines the accuracy level of the model, its computational cost, and its capabilities.
In this work, we consider two models. The first model is based on three interactions: (i) repulsive, (ii) attractive, and (iii) electrostatic. The first two components follow the functional form of the Mie potential [72] with a repulsive C 9 term and an attractive C 6 term. The electrostatic component is described with Coulomb's law and fixed partial charges. Hence, this model assumes an isotropic-pairwise-additive form, and we will refer to it as "IPA model". The second model considered in this work is based on multipole electrostatics and polarizable atoms, therefore violating isotropy and pairwise additivity. We will refer to this anisotropic-non-additive model as "ANA model".
The IPA model is similar to the functional form found in the most widely used FF [19], [45], [73]- [77]. The functional form of the ANA model, on the other hand, is an attempt to explore the limit of a purely classical model through an implicit description of the electron distribution based on atomic multipoles and induced dipoles. Its functional form is similar to the class of polarizable and QM derived FF [55], [56], [78]- [83]. In the following section, the components used in the respective models are described.
For consistency, the following notation is used: Capitalized letters refer to pairwise parameters, i.e., parameters given for a pair of atoms. Small letters are used for atomic parameters, i.e., parameters assigned to one specific atom. Subscripts are used to further clarify interaction partners. In general, the indices i, j iterate over each unique pair of atoms.

Potential-Energy Terms: IPA Model
The IPA model for the non-bonded potential energy includes three components: attractive-repulsive described with the Mie potential [72], and the electrostatic interaction described with atomic partial charges interacting through Coulomb's law,

Attraction-Repulsion Potential
As the simplest potential with two parameters that can reproduce the qualitative features of the dissociation of uncharged atoms, the Mie potential [72] is used, where C n and C m denote the coefficients used to describe the strength of the repulsion and the attraction, respectively, and r ij the distance between two atoms i and j. For n = 12 and m = 6, the well-known Lennard-Jones potential is obtained [84]. While the attractive part is often set to m = 6 motivated by the asymptotic behaviour of the dispersion interaction, there is no such obvious choice for the repulsive part. In the past, n = 12 was commonly chosen for its computational efficiency. In this work, a softer n = 9 repulsive interaction is used as in Ref. [85]. This choice allows for a more accurate description of the repulsive interaction while still retaining computational efficiency and comparability with the more common C 12 − C 6 formulation.

Electrostatics
Electrostatic interactions in the IPA model are treated on the basis of atomic monopoles (partial charges), with monopoles q i and q j , and the vacuum permittivity 0 . In this work, the monopoles are obtained from our previously introduced equivariant GNN model [63], which was trained on minimal basis iterative Stockholder multipoles (MBIS) [86].

Predicted Parameters
The IPA model predicts a set of C 6 and C 9 parameters for each atom pair, whose features were represented with the permutation-invariant pairwise feature combination shown in Eq. (6). No prior knowledge such as baseline or default parameters is used in the training of the IPA model. As mentioned above, the monopoles q i are obtained from a separate GNN model [63] and not further modified for the present work.
The atomic monopoles and the parameters of the Mie potential are considered fixed parameters, i.e., they remain constant during a simulation and do not change in response to changes in the molecular geometry. This setting was chosen to be comparable with existing FF [19].

Potential-Energy Terms: ANA Model
Unlike the IPA model, the second model includes explicit treatment of polarization effects and anisotropy, resulting in an anisotropic and non-additive model, abbreviated as 'ANA' model. The model is motivated by the desire to develop a fully classical description that performs with an accuracy expected from semi-empirical methods. In addition, we aim to demonstrate the power of the proposed parametrization strategy through its application to a model with several inter-dependent components. At its core, the ANA model is based on an implicit description of the electronic structure through the use of atomic multipoles and a polarization model. Adding additional interaction terms allows for a more detailed decomposition of the total energy. Consistent with the decomposition by symmetry adapted perturbation theory (SAPT) [87], components for the dispersion, electrostatic, induction, and exchange potential energy are used. As a further benefit, it is possible to include SAPT terms in the fitting procedure.
The functional form of the non-bonded potential energy of the ANA model is inspired by previous work on polarizable FF and intermolecular potentials [56], [78], [79], [88].
The dispersion interaction (V disp ) is described with dispersion coefficients in conjunction with the Becke-Johnson damping model [83], [89]. Following the model used in AMOEBA/AMOEBA+, induction is included through the Thole damping model (V ind ) and a charge-transfer potential (V ct ) [78], [79], [90], [91]. Exchange and electrostatic interactions (V ex and V ele ) are based on anisotropic potentials derived for atomic multipoles following the work by Rackers et al. [80]- [82], [91]. Atomic multipoles are obtained with our recently introduced equivariant GNN [63]. This model was not further modified for the present work. A detailed description of the terms is given in the Appendix A2.1.

Predicted Parameters
As in the IPA model, atomic multipoles were obtained using our previously developed equivariant GNN [63]. In total, five atomic parameters and five pairwise parameters were predicted by the ANA model. Atomic parameters include the atomic polarizability (α in Eq. (A10)), the exponential factor for the electrostatic damping function (b in Eqs. (A2) and (A3)), the exponential factor used in the damping function of the exchange potential (b in Eqs. (A16), (A17) and (A18)), the scaling factor used to adjust the strength of the exchange potential (k in Eq. (A20)), and the number of valence electrons (q val ), which is added to the negative atomic monopole to obtain a scalar that replaces the atomic monopole in the anisotropic exchange potential. Further, pairwise C 6 , C 8 and C 10 parameters were independently predicted for each atom pair as well as the exponents for the short-range induction potential (A and B in Eq. (A15)).

Differentiable Force Field
To achieve end-to-end differentiability, all FF terms and parametrization models were implemented in Tensorflow (version 2.5) [66], [92], taking advantage of its automatic differentiation capabilities as well as (batched) GPU accelerated computation. The particle-mesh-Ewald method (PME) [93], [94] implemented in OpenMM (version 7.7) [95] was used to obtain the long-range electrostatic contributions for periodic systems.

Graph Construction
The approach applied in this work is based on graphs constructed from the molecular topology, referred to as 'topological graphs'. Topological graphs do not include information about the Euclidean distance between atoms but only atomic connectivity. Including geometrical information could be advantageous for certain applications but would require frequent recalculation of parameters, which would limit the performance. In addition, model robustness might suffer from insufficient sampling of intramolecular degrees of freedom. On the other hand, topological graphs may be ill-defined for certain cases and are unable to describe phenomena such as bond forming and breaking. For the envisioned application, i.e., the classical description of molecular motion, the shortcomings of topological graphs are acceptable while presenting a robust and efficient solution. We note that some degree of conformational dependence is present in the overall approach due to the GNN model used for the prediction of atomic multipoles [63].
No chemical concepts such as bond types or hybridization states were included in the graphs. Hence, graphs only contained a description of the element type of an atom and the presence of a covalent bond between two atoms. Further information regarding the construction of molecular graphs is given in the Appendix A2.2.

Loss Weighting
Models were optimized by minimizing the mean square error L between the predicted intermolecular potential energy V pot,θ and a reference intermolecular potential energy V pot,ref , which was used as the target property with N denoting the batch size, and i iterating over each sample of the batch. w i is a scalar used to weight the contribution of each sample. The importance of each sample was scaled according to its Boltzmann weight in the following manner Figure 2: Data sets used for training, validation and testing of the IPA and ANA models. The number of molecular systems in the respective data set is given below the set indication. For dimers in the DES5M data set [64], a variable number of samples is found, indicated by 'xN'. Details on the data sets are given in Appendix A2.6.
Samples beyond the equilibrium distance r eq were weighted with w i = 1, and samples closer than the equilibrium distance weighted with the Boltzmann weight of the difference between the potential of the given sample and the potential of the equilibrium sample for a given system. The inverse-temperature β was used to determine the relative importance. To give more importance to near-equilibrium energy samples towards the end of the training procedure, an exponential decay was used to simulate annealing with a decay rate γ and an initial temperature T 0 and a stopping temperature T min .

General Training Strategy
The recently published DES5M data set [64] (see Appendix A2.6.1) was used as main data source for the model training ( Figure 2). The data set includes spin-network-scaled MP2 (SNS-MP2) [9], [96], [97] intermolecular potentials and SAPT0 components [87], [98] for a large number of small-molecule dimers in vacuum. In total, 113 800 dimer sets were included with 100 000 sets randomly selected for training and the remaining 13 800 used for validation. Training was performed over 512 epochs. During each epoch, a total of 1024 batches were presented. Each batch contained one dimer set, i.e., all interaction potentials of one dissociation curve or set of clusters for the same two monomer molecules.

IPA d Model
The IPA model was first optimized against dimer intermolecular potential energies in vacuum from the DES5M data set [64] through minimization of the expression found in Eq. (11), with SNS-MP2 energies serving as V pot, ref . Annealing was performed for the loss weight term w i (V pot, ref (x i )) using the described exponential decay with T 0 = 8 000 K, γ = −8 · 10 −3 , and n being incremented after every epoch. Annealing was stopped at T min = 400 K. Annealing was of particular importance for the IPA model as the attraction-repulsion potential used is not able to describe very short-ranged interactions accurately.
The IPA model was regularized with the following term MLFF During the training on dimers, an additional term L C6 = C −1 6 was added. Both regularization terms were averaged over all interaction pairs of a molecule. The model obtained in this manner is called IPA d in the following

IPA d+c Model
The IPA d model was optimized in a second step with respect to intermolecular potentials of molecular crystals. Intermolecular potentials were calculated with Quantum Espresso [99]- [101] using the B86bPBE functional in conjunction with the XDM dispersion correction [89], [102]- [105]. For this step, a total of 11 489 molecules from the CSD database were used [106]. 10 000 molecules were randomly assigned to the training set and the remaining 1 489 molecules formed the validation set. Each molecule entailed five intermolecular potentials, which served as one batch. These five potentials were obtained from the relaxed crystal structure and through the expansion/contraction of the lattice of the relaxed structure. Further details on this data set are given in the Appendix A2.6.2.
Again, the mean square error between the predicted and the reference intermolecular potentials was minimized. The electrostatic component was independently calculated with the PME method [93] and not optimized. Fixed partial charges were obtained as described in Section 4.7. The loss was weighted with Unlike the weighting function used for dissociation curves in Eq. (11), the weight was only set to one for the equilibrium structure. In addition, the potential-energy terms were scaled by the number of atoms in the respective molecule. Loss weighting is necessary since contracting/expanding the relaxed crystal structures resulted in highly unfavourable structures in certain cases. Training was performed over 512 epochs. During each epoch, 512 randomly sampled batches were presented. Each batch contained five intermolecular potentials of one specific molecule. Throughout, the temperature was set to T 0 = 128 K. A cutoff of 10 Å was used for nonbonded interactions.

Training Strategy: ANA Model
Optimization of the ANA d model followed the same procedure used to optimize the IPA d model based on the dimer data set, except for the following differences. Since the functional form of the ANA model allows for a more accurate description of short-range interactions, the annealing schedule was modified to T 0 = 40 000 K, γ = −7.5·10 −3 , and T min = 2 000 K. Furthermore, the L C 6 and L C 9 terms were replaced with additional loss terms for the energy components as the ANA model permits a decomposition of the total energy into components, which can be related to SAPT components [87], [98].
with λ iterating over all SAPT components, i.e., exchange, induction, electrostatic and dispersion. The weights w i were calculated for the total reference potential. The SAPT loss term was scaled by a factor κ and added to the total loss in Eq. (11). The weight for the total energy component with respect to the SNS-MP2 calculation was kept at 1. κ was initially set to 0.5 and reduced to 0.01 after 256 epochs had passed. We note that the inclusion of the SAPT loss term is a double-edged sword. Preliminary investigations indicated that the SAPT loss term serves on one hand as a regularizer, which also accelerates convergence. On the other hand, larger values for κ limit the degree of error cancellation between potential-energy terms. It is further important to keep in mind that the accuracy of the employed SAPT0 method is lower than SNS-MP2. Specifically, for the considered subset of the DES5M data set, a MAE of 4.19 kJ/mol and 0.77 kJ/mol was found when using the weighting in Eq. (15) with T= 2000 K. With a mean error of −4.11 kJ/mol and −0.71 kJ/mol, SAPT0 overbinds relative to the SNS-MP2 results. This observation is consistent with previous benchmarks [87].

Calculation of Condensed-Phase Properties
The IPA d+c model was also applied to calculate properties of condensed-phase systems. The simulations were performed under periodic boundary conditions using OpenMM [95]. Up to 32 conformers were generated for each molecules with the ETKDG conformation generator [107] as implemented in the RDKit [108]. An RMS pruning threshold of 0.1 Å was used. Monopoles were predicted for each conformation using our previous equivariant GNN [63], and then averaged over the conformational ensemble to yield the fixed partial charges for the simulation. The same partial charges were used for the simulations in the condensed phase and in the gas phase. Bonded terms were parametrized with OpenFF 2.0 [109] since the IPA model handles only intermolecular interactions. Further details regarding simulation setups are given in the Appendix A2.5.

Lattice Energy
In the present work, the lattice energy V lattice is approximated as ignoring the contribution of the intramolecular interactions. Z refers to the number of molecules in the unit cell, and V pot, inter to the total intermolecular potential energy for a unit cell with Z molecules under periodic boundary conditions. V pot, inter was calculated for the experimental geometries without relaxation. A list of CSD codes is given as Supporting Information.

Heat of Vaporization
The heat of vaporization was computed from the difference between the mean potential energy in the gas phase V pot, gas and the mean potential energy per molecule in the condensed phase V pot, liq corrected by a factor of RT, where R is the gas constant, and T the absolute temperature.

Density
The condensed-phase density was calculated as the total mass m box in the simulation box divided by its

Static Dielectric Constant
Static dielectric constants were obtained from the fluctuation of the dipole M of the system as described in Ref. [110],

Isothermal Compressibility
Similarly, isothermal compressibilities κ were obtained for fluctuations of the system volume V box [110], where P is the system pressure.

Thermal Expansion Coefficient
Thermal expansion coefficients α were computed via the following relation [110], where H l is the total enthalpy of the box.

Results and Discussion
The performance of the IPA d , IPA d+c and ANA d models to describe intermolecular interactions were investigated for a wide range of systems and environments. For the IPA d+c model, emphasis is put on the performance in condensed-phase systems, which are generally challenging for non-classical models.
With its more sophisticated functional form, the ANA model is applied to intermolecular potentials of small-molecule dimers, which permits a direction comparison with first-principle and DFT methods. Finally, further explorations of parameters learned by both approaches are presented. We note that the ANA model has so far not been parametrized and applied to condensed-phase systems.

Intermolecular Potentials in Vacuum
Small-molecule dimers are the largest systems that can be treated with highly accurate wave-function methods. As such they present a valuable validation case to probe the accuracy of the description of specific interactions. For this purpose, several established non-covalent interaction benchmarks were taken as test sets from the Biofragment database [111]- [117] and the non-covalent interaction (NCI) atlas [118]- [122] (Figure 2). In addition, the models are tested on the supramolecular S12L [123] data set. For systems in the S12L, binding energies calculated with quantum diffusion Monte Carlo were used as reference [124]. In all cases, only systems consisting of neutral monomers with more than two atoms were included. The performance results on the benchmarking sets with a total of 10 894 unique data points are summarized in Table 1. The full error statistics is given in  [64]. The full error statistic is given in Tables S1 -S3 in the Supporting Information.
The results of the three models investigated in this study are shown in Table 1. Consistent results over these diverse sets are observed. For comparison, reference values are given for the data sets in the Biofragment database where available ( Table 2). Three methods were chosen to represent classical force fields (CHARMM General FF (CGenFF) [125]), semi-empirical models (PM6-DH2 [126]), and DFT methods (PBE0-D3 [127]- [129]). The IPA models perform comparable to previously reported results for empirical (FF) or semi-empirical models. In general, the simpler models are outperformed by the ANA d model, which achieves for some data sets an accuracy comparable to dispersion corrected hybrid functionals like PBE0-D3BJ. Exceptions are the BBI data sets, where the IPA d model performs better than ANA d , and the S12L data set for which IPA d+c outperforms the other models.  Comparison between IPA d+c and IPA d suggests that additional training on crystal data points does not necessarily lead to a strong negative impact on the description accuracy of dimers in vacuum. This is surprising for two reasons: (i) The method used to compute the potential energy of crystal structures is considerably less accurate than the methods used to compute the dimer interaction potentials in the benchmarking data sets and the DES5M training data set. (ii) Polarization effects, for which the IPA models must implicitly account, differ between vacuum and condensed-phase environments. This effect is visible in the case of the S12L data set, where IPA d+c clearly outperforms the IPA d model. Interestingly, training on crystal potential energies appears to provide a regularizing effect for the IPA d+c model, offering not the most accurate but the most consistent results over a wide range of applications beyond dimers (see below).

MAE for interactions in vacuum
The results for the five neutral systems of the revisited S12L data set are worthy of special attention. For these supramolecular systems, a consistent overestimation of the interaction potential is observed for the ANA d model with a mean error equal to the mean absolute error (Table S1 in the Supporting Information). A large contribution might be due to the missing treatment of many-body effects, which was noted by the creator of S12L [123] and further discussed for the revisited values [124]. Estimation of the three-body dispersion contribution would explain approximately half of the overbinding reported here [123]. Previous work observed similar effects with ML models [55]. However, despite using a dispersion correction with many-body effects, Ref. [55] reported comparable results to our models. The performance could only be improved through the inclusion of the same structures in the training set [55]. In a similar vein, a false balance between many-body effects (e.g., between the potential energy due to induced dipole and pairwise interactions such as dispersion and the short-range induction potential) might contribute further to this effect. A weak overestimation of pairwise interactions might not be noticeable for small systems but could be amplified for larger systems. Overbinding of SAPT and its components could further contribute to this observation. Thus, the S12L data set may be exposing a possible limitation of polarizable models and small molecule data sets. In this context, we also note a recent study where the authors observed considerable discrepancies between reference methods for large complexes [132]. Investigating which effects need to be included to permit the generalization to larger systems and the condensed phase is therefore an important open question. Developing additional benchmarking data sets with a broader coverage of the space from medium to large structures could be MLFF highly beneficial to validate such efforts.

Model Parameters
The influence of feature size and the number of GNN layers was investigated for the IPA d+c model. For this, a range of models with varying feature dimensionality and GNN layers was trained on the DES5M data set as described in the Methods section. Shorter epochs were used, i.e., 512 batches instead of 1 024, while all other parameters remained unchanged. The performance was assessed on the basis of absolute weighted (T= 400 K) errors for the DES5M validation set, which contains 150 molecules resulting in 13 800 dimer sets and a total of 471 149 data points. The results shown in Figure 3 support the inclusion of information about the atomic neighbourhood. The first two coordination shells are particularly beneficial, with marginal benefits or even detrimental effects for larger numbers of steps. Purely local models (number of steps = 0), which only include information about the atomic element, perform considerably worse than non-local models. Nevertheless, even for this most simple case, acceptable performance is observed, which might be due to the monopoles already capturing the most important features of atomic interactions. Interestingly, increasing the numbers of parameters (i.e., feature dimensions) of the model itself does not seem to provide benefits beyond a given range. It is important to note that the molecular graphs for the considered DES5M validation set contain comparatively small monomers with a mean graph diameter of 4.2 and a maximum graph diameter of 8. Some of these conclusions might change for larger molecules. It could be particularly illustrative to investigate the potential benefit of additional graph-updating steps for large aromatic systems with non-local effects due to substituents such as nitro groups or long and branched systems.

IPA d+c Model -Condensed-Phase Properties
Reproducing properties of condensed-phase systems is an important validation task for intermolecular potentials. QM methods are in general not feasible to simulate such systems, and the wide range of interactions poses considerable challenges to ML models. In the following, the IPA d+c model is validated on various condensed-phase properties of pure organic liquids.

Intermolecular Potentials in Crystals
Intermolecular potential energies calculated for a wide range of molecular crystals were used to parametrize the IPA d+c model in addition to the dimers in vacuum. Table 3 shows the mean absolute error on the molecular crystals for the IPA d+c and the IPA d models. Only equilibrium structures of the crystals were included. The training set and test set contained 10 000 and 1 507 data points, respectively. As can be expected, the IPA d model performs considerably worse than the IPA d+c model. As the latter model performs also comparatively well on dimers in vacuum, these results point at a general advantage for the training of IPA-type models on both gas-phase and condensed-phase data. For this reason, only the IPA d+c model is used in the following.  To obtain a picture of the range of van der Waals parameters predicted by the model, Figure 4 shows the C 6 and C 9 parameters from the IPA d+c model for all pairwise interactions in the crystal data set. Further information is given in Figures S1 and S2 in the Supporting Information, showing the resulting well-depth and the minimum distance by atom pairs. The nature of the electrostatic interaction is indicated by the color, with blue for attractive interactions and red for repulsive interactions. Interestingly, a large part of the correlation between C 6 and C 9 parameters might be captured by a power law. This observation could potentially be used to construct FF based on a single parameter and appropriate scaling laws, for instance based on a notion of atomic volumes. Furthermore, the predicted C 6 and C 9 parameters form a continuum over large ranges. While there are some distinct islands, in particular interactions with hydrogens in the bottom left, the results in Figure 4 are nevertheless an indication that the model takes advantage of continuous atom types.

Lattice Enthalpies: X23 and G60
Performance on crystal structures was evaluated on the X23 benchmark [133], [134] using the revised values found in Ref. [135] as well as the G60 data set [136] (Table 4). Two settings were considered: For setting (1), labelled with 'Non-relaxed', the lattice energy was approximated with the intermolecular potential for the non-relaxed experimental geometries. In setting (2), labelled with 'Relaxed + intra', structures were relaxed and the potential-energy difference between the gas-phase minimum conformation and the crystalline phase minimum was included. For the second setting, six systems had to be excluded due to problems with the bonded terms from OpenFF (i.e., X23: CO 2 and UREAXX12; G60: CTMTNA03, METNAM08, MTNANL, OCHTET13).
For both data sets, the cohesive energy is systematically underestimated with the IPA d+c model (i.e., positive mean errors in Table 4) and the MAE is above chemical accuracy. Nevertheless, the IPA d+c model reproduces lattice energies more accurately than most existing models reported so far in the literature, e.g. DFTB-D3 with a MAE of 10.38 kJ/mol for the X23 data set [137], except for some of the best performing dispersion corrected DFT functionals such as PBE0-MBD with a MAE of 3.9 kJ/mol on the X23 data set [133], [134]. This shows that the accurate description of the lattice energy of molecular crystals remains a particularly challenging problem, which will continue to serve as an important reality check, specifically for models that are not parametrized on condensed-phase systems.  Table 4: Mean absolute error (MAE) and mean error of the IPA d+c model for the X23 and G60 benchmark sets based on revised values provided by Ref. [135] and reference values collected in Ref. [136]. The full error statistics are provided in Table S4 and S5 in the Supporting Information.

Pure Liquid Properties
Properties of pure organic liquids are commonly used to validate classical FF, serving in many cases also as parametrization targets (e.g., [75], [138], [139]). To further explore the performance of the IPA d+c model beyond vacuum and the crystalline phase, applications to the liquid phase are shown in the following. Three benchmarks covering a wide range of systems and properties are considered for this: (i) 13 sulfur compounds taken from the publication of the OPLS4 release [140], (ii) 29 molecules containing H, C, O taken from a recent investigation of condensed-phase parametrization of OpenFF [139], and (iii) 57 organic compounds from the GROMOS 2016H66 validation [138]. Unlike the referenced FF, parametrization of the IPA d+c model did not include experimental liquid properties such as the density or the heat of vaporization. Therefore, pure liquid properties present an interesting test case for this model. For the liquid simulations with the IPA d+c model, bonded interactions (i.e., bonds, angles, dihedrals, and 1-4 nonbonded interactions) were treated with OpenFF 2.0 [109], while all other interactions (i.e., nonbonded terms) were treated with the IPA d+c model. A detailed description of the simulation protocol is given in the Appendix A2.5.
Sulfur Compounds. Systems containing sulfur are challenging for fixed-charge FF due to the polarizability of sulfur and the presence of higher-order multipole components. Recent work on the OPLS4 FF [140] improved the performance on several challenging motifs, including sulfur interactions and σ-holes. Specifically, OPLS4 improved the RMSE of the heat of vaporization (H vap ) for 13 sulfurcontaining systems by more than 1 kJ/mol (i.e., 0.3 kcal/mol) compared to the previous OPLS3 version [140]. With an RMSE of 2.5 kJ/mol for H vap , the IPA d+c model performs comparable to OPLS4 (Table  5). This result is remarkable for several reasons. First, unlike the IPA model, OPLS4 uses virtual-sites to represent lone pairs and anisotropic Lennard-Jones interactions, which were specifically introduced to improve the description of systems containing sulfur and halogens. Second, the OPLS FF family was specifically parametrized with respect to liquid properties such as H vap [75], whereas the IPA d+c model was only trained on QM potential energies. Third, the intramolecular potential of OPLS4 is jointly optimized with the intermolecular potential, allowing for a larger degree of consistency between the two parts. It is likely that considering the above points in future work on the IPA model could result in further improvements. 26.0 - Table 5: Root-mean-square error (RMSE) for the heat of vaporization (H vap ) and density (ρ) for 13 sulfur-containing compounds investigated in Ref. [140]. Values for OPLS4 from Ref. [140] are given as comparison and were converted from kcal/mol to kJ/mol by a factor of 4.184. Note that densities were not reported for OPLS4 in Ref. [140]. The full error statistics are provided in Table S6 in the Supporting Information. The individual numerical values are given in Table S11.
Test systems from OpenFF. To gain a better understanding for the role of the bonded terms taken from OpenFF, results for 29 pure liquids from a recent benchmark of OpenFF are presented here [139]. The compounds contain only H, C, and O. The referenced work is particularly interesting for its investigation of opposing forces during the parametrization with respect to mixing enthalpies, vaporization enthalpies, and densities. Only pure liquid properties were considered here.
As can be seen in Table 6, similar errors are observed for the IPA d+c model and the standard OpenFF 1.0. Re-optimization of the FF with respect to the pure liquid properties of the training set (also compounds containing only H, C, and O) improved the accuracy of OpenFF 1.0 considerably [139] ( Table  6). The authors observed thereby opposing gradient components for the simultaneous optimization with respect to densities and heats of vaporization. No further liquid properties were considered in Ref. [139] (such as dielectric permittivity, thermal expansion coefficient, etc.). It would thus be interesting to see the performance of the re-optimized OpenFF (termed 'Pure only') on other properties.
The fact that all models shown in Table 6 use the same functional form with very similar bonded terms may indicate that liquid properties cannot be reproduced more accurately without either improving the description of the bonded interactions and/or extending the functional form, for instance through the use of a polarizable model and multipoles. We note that the simulations with the IPA d+c model used bonded terms from OpenFF 2.0 were employed, while OpenFF 1.0 was employed in Ref. [139].   [139]. Values for OpenFF were taken from the referenced publication. The label 'Pure only' refers to a version of OpenFF 1.0, which was re-optimized on the densities and heats of vaporization of the training set in Ref. [139]. The full error statistics are provided in Table S7 in the Supporting Information. The individual numerical values are given in Table S12.
Test systems from GROMOS 2016H66. The 57 pure liquids from the GROMOS 2016H66 [138] release include extended coverage of the chemical space and additional properties such as the isobaric thermal expansion coefficient (α), the static relative dielectric permittivity ( ), and the isothermal compressibility (κ). For the considered properties, the IPA d+c model performs comparable to the 2016H66 FF (Table 7). While H vap is less accurately reproduced by the IPA d+c model, smaller errors are observed for the remaining properties. Note that 2016H66 was parametrized on the H vap and density values of 27 of the considered 57 molecules. The observation that the IPA d+c model outperforms 2016H66 on properties that were not used for its parametrization may be an indication that the 'bottom-up' approach MLFF of the IPA approach, focusing on the reproduction of the PES, is a valuable parametrization strategy. It further demonstrates that there is still room for improvement for the fixed-charge FF model.  Table 7: Root-mean-square error (RMSE) for pure liquid properties of 57 systems used in the calibration and validation of the GROMOS 2016H66 FF [138]. Values for GROMOS 2016H66 were taken from the referenced publication. The full error statistics are provided in Table S8 in the Supporting Information.
The individual numerical values are given in Table S13.

Learned Parameters
As the ANA approach cannot be used yet for condensed-phase simulations, we validated the ANA d model by comparing the predicted molecular polarizabilities and intermolecular C 6 dispersion coefficients to experiment. The data set consists of molecular polarizability values for 87 compounds [141]- [146] and C 6 coefficients for 231 molecular pairs [147]. Since the ANA model predicts atomic parameters, the molecular polarizability is obtained as the sum of all atomic contributions. Dispersion coefficients are summed over all intermolecular atom pairs as in previous work [104], [128], [148]. Figure 5 shows the comparison between predicted and experimental values. Mean absolute relative errors of 19.4 % and 2.1 % and Spearman correlation coefficients of 0.95 and 0.94 are observed for the C 6 coefficients and the molecular polarizabilities, respectively. For the C 6 coefficients, two sets of outliers are found. One set includes interactions with tetrachloromethane (labeled with CCl 4 ) while the second set includes interactions with chloromethane (labeled as CH 3 Cl). Figure 5: Comparison of the C 6 dispersion coefficients (right) and molecular polarizabilities (left) obtained from experiment and predicted by the ANA model. Experimental values for the molecular polarizability of 87 compounds were taken from Refs. [141]- [146], and experimental C 6 coefficients of 231 molecular pairs were taken from Ref. [147].

MLFF
In Figure 5, we can see a systematic overestimation of C 6 coefficients by the ANA d model, which could be a result of two factors: First, the value of the damping function in the Becke-Johnson scheme [83], [89] is related to the dispersion coefficients through Eqs. (A9) and (A8). Since all dispersion coefficients are treated as free parameters, compensation effects may arise. Second, the magnitude for higher-order coefficients might be too weak compared to other components. As a result, the relationship assumed in Eq. (A9) might not hold in the same way as for coefficients obtained from DFT densities. Similarly, a weak overestimation is found for the molecular polarizabilities. As already noted, this may be caused by a false balance between the short-range induction potential (V ct ) and the long-range model. Inclusion of induced dipoles in the exchange potential and the parameter used for the Thole damping might also play a role. In addition, it should be noted that topology-based parameters must compensate for damping effects due to the surrounding environment, which are not accounted for but may play an important role [149].
Overall, the results show that the ANA d model predicts physically meaningful parameters from scratch given a physically motivated functional form. As such, including constraints through the use of predefined functional forms or known (in-)equalities and parameter-relations might be the method of choice to regularize ML models applied to physical problems. Further improvements may take into account the relationship between polarizabilites and dispersion coefficients with atomic volume ratios, which has been shown in several studies [148]- [152]. Using an independent model to handle polarizabilities and dispersion coefficients could not only introduce sensitivity to the surrounding environment and conformational changes but also reduce the number of fitted parameters and the resulting interdependencies.

General Discussion and Outlook
The presented results attest to the power and feasibility of the proposed approach to use ML models to predict the parameters of a physically motivated functional form of a FF. Most importantly, it indicates that a general FF, which performs well on a wide range of applications, might be attainable. To explore whether the ML model can truly learn parameters from the atomic environments alone, no baseline parameters were used in the present work. However, introducing such baseline parameters as well as appropriate scaling of parameters may improve the model in the future. For the IPA model, this could be accomplished by using an existing FF as baseline and tuning the parameters with the ML model. For the ANA model, inclusion of additional non-fitted parameters such as the atomic volumes could not only introduce conformational sensitivity but also help to reduce the number of parameters that have to be fitted. The atomic polarizabilities and dispersion coefficients, in particular, are suitable candidates for such a strategy as they can be obtained from atomic volume ratios given an appropriate model [149]. Including further information such as spectral data from experiments or higher-order derivatives from theoretical calculations could also prove fruitful. This might be particularly interesting for the description of intramolecular interactions.
Model accuracy is mostly limited by two factors: First, an implicit description of electrons is only possible if effects due to changes in the electron distribution can be captured accurately. For this, an accurate description of multipoles is required. Improving multipole prediction can likely improve the model performance further. Second, the presented ML-based approach crucially depends on the accuracy and availability of reference data. Expected improvements in hardware, software, and methodological advances may provide access to highly accurate calculations of condensed-phase properties in the near future.

Conclusion
A formalism based on automatic differentiation for the parametrization of classical FF was introduced. The proposed formalism cannot only be viewed as the generalization of commonly used FF definitions and parametrization procedures, it also describes the regularization of ML models through physics. To showcase its strength, the method was applied to the parametrization of an isotropic-pairwise-additive FF (IPA model) as well as a polarizable FF (ANA model). The performance of the models was first assessed on intermolecular interaction energies of dimers in the gas phase. Both approaches perform on par or MLFF better relative to comparable methods on several common benchmark data sets, while requiring little computational efforts and human intervention. Importantly, the models were exclusively parametrized with respect to the PES of reference QM calculations, without the addition of experimental data. The IPA d+c model was furthermore applied to the calculation of condensed-phase properties (lattice energies and pure liquid properties), while the ANA d model was validated by comparing molecular polarizabilities and dispersion coefficients to experimental values. We found that the IPA d+c model, i.e., a fixed-charge FF parametrized on intermolecular potentials from DFT calculations, can provide consistent results over a wide range of systems, ranging from dimer interaction potentials in vacuum to pure liquid properties. While a completely classical description of molecular forces has clear limitations (e.g., no bond formation or breaking), it might still be highly competitive to semi-empirical methods as shown for the ANA model. In particular, the implicit description of the electron distribution through atomic multipoles offers a very promising alternative to semi-empirical models, providing a comparable level of accuracy but only requiring a fraction of the computational cost.
Automatic differentiation presents a powerful tool for the development of parametrized models, which could also be applied to problems other than molecular interactions. As shown, this approach accelerates and simplifies the parametrization process of classical FF and can take advantage of large data sets. In combination with ML based techniques, such as the presented GNN-based atom-typing, the best of both worlds can be harvested. A universal optimisation toolkit combined with robust and physically-constrained models. In future work, the exploration of additional FF terms and the application of the ANA model to condensed-phase systems will be investigated. In addition, employing the end-to-end differentiable approach to the parametrization of bonded/intramolecular interactions will be explored.

Acknowledgment
This research was supported by the NCCR MARVEL, a National Centre of Competence in Research, funded by the Swiss National Science Foundation (grant number 182892).

Software and Data Availability
The intermolecular potentials for 11 489 molecular crystals which were used during training of the IPA d+c model are available in the ETH research collection (https://doi.org/10.3929/ethz-b-000549359). A GitHub repository including weights and models used to produce the results in this work can be accessed under the following link: https://github.com/rinikerlab/GNNParametrizedFF.

A2.1.1 Electrostatics
The electrostatic interaction is described through the use of multipoles up to quadrupoles. Multipoles were in all cases obtained from our previous equivariant GNN [63] for the prediction of atomic multipoles, which was trained on MBIS reference data [86]. The implementation follows the formalism proposed by Refs. [153]- [155]. Following the formalism introduced in Ref. [153], the total potential energy due to the interaction of point multipoles at site i and site j is obtained as with r ij = r i − r j and r ij = | r ij |. The radial functions B l (r) are defined as B l (r) = (2l − 1)!!/r 2l+1 and the coefficients G l ( r) arise through the interactions between components of two multipole sites. For the considered case of a treatment up to quadrupoles, terms up to l = 4 are included. Coefficients G l ( r) can be thought of as contributions due to the interactions of multipoles of a given order. A list of G l ( r) is given in Ref. [153]. Deficiencies of the multipole description at short ranges are compensated through the use of the charge penetration model introduced in Ref. [80]. In this model, effects of interactions between charge distributions are treated through the use of damping functions. Specifically, a damping function for the interaction between a charge distribution and a point charge and a damping function for the interaction between two charge distributions are introduced. The damping parameter b i describes the extent of an exponentially decaying charge distribution centered on atom i. The above damping functions give rise to damping coefficients of a given order λ l (r) which are given in the Supporting Information of Ref. [80]. Combining the damping coefficients with the radial functions B l (r) gives rise to the damped radial functions In the charge penetration model, the standard radial functions are replaced with the damped radial functions. In addition, core-core, core-multipole, and multipole-core interactions are included. Thus, the complete description of the electrostatic potential is obtained as and G l M M describe the core-multipole, multipole-core, and multipole-multipole interactions, respectively. The core-multipole and multipole-core coefficients are obtained by replacing the monopole with the respective core charge.

A2.1.2 Dispersion
The dispersion interaction is described based on the formalism used in the XDM model proposed by Becke and Johnson [83], with dispersion coefficients C n and a damping function depending on a damping parameter R vdW , which is obtained as follows [89], with two positive parameters a 1 and a 2 and In the present work, a 1 was set to 1 and a 2 to 0, i.e. R vdW = R c .

A2.1.3 Induction
Polarization is treated based on the Applequist model [156] including the modification proposed by Thole [90] and follows the formalism described by Stone [4]. The long-ranged component is treated on the basis of atomic dipoles, which are induced by the external electric field and scaled by the atomic polarizability [4] where F D gathers the electric field components at each atom formed by the static multipoles of the surrounding molecules, i.e., only intermolecular contributions. The polarizability matrix B is formed as with the atomic polarizability α i and the elements T ij of the dipole-dipole interaction matrix. The 3N × 3N polarizability matrix is inverted to obtain the induced dipoles. Given self-consistently induced dipoles, the potential energy due to induction is given as the inner product with the external field [4] To prevent the divergence of induced dipoles ('polarization catastrophe'), elements of the polarizability matrix B are damped based on the modifications proposed by Thole analogously to the damping function used to model charge penetration effects [90].
The exponential damping function as used in AMOEBA was used for this purpose [78] f Thole using a damping factor a and the polarizability-normalized distance As in the original AMOEBA FF, the damping factor a was globally set to 0.39 [78]. In addition, a charge transfer potential was added to improve the treatment of polarization at short ranges. This potential is based on work proposed for the AMOEBA+ FF [79] using an exponential form where A describes the strength of the interaction and C is used to approximate the degree of electron density overlap between the respective atom pair.

A2.1.4 Exchange
The exchange interaction is treated with the anisotropic repulsion model proposed by Rackers et al. [82].
Using atomic multipoles, their work derives a description for the overlap between two atoms analogously to the electrostatic interaction between atomic multipoles leading to the following expression consistent with the expression obtained by Salem [157]. Specifically, the damping function used to construct the damped radial functions B damp l as shown in Section A2.1.1 is replaced with the following damping function for the case b i = b j and leading to the radial function for the exchange potential with higher order radial functions B damp 0 (r) following analogously to the damped radial functions presented in the description of the electrostatic potential.
The overlap defined as is then used to obtain the exchange potential energy contribution, with k i being the relative size of atom i, and S 2 the multipole derived orbital overlap. We note that b used to damp the exchange interaction (Eqs. (A16) and (A17)) and the electrostatic interaction (Eqs. (A2) and (A3)) are treated as independent parameters despite relating to the same underlying feature, i.e., an exponentially decaying charge distribution. Induced dipoles are added to the static dipoles. As in the original work [82], the monopole is replaced with an additional atomic parameter q val , which weights the influence of the multipole interaction coefficients present in G l ( r). Following Ref. [82], q val is set to 1 for all hydrogens and > 2 for all other elements. This parameter is added to the negative monopole, yielding the final q ex parameter, which is used in place of the monopole used to compute the multipole interaction coefficients in G l ( r).

A2.2 Graph Construction
Since topological information is not available for all data sets, graphs were constructed from monomer coordinates. Each graph was built by adding a node for each atom and an edge between bonded nodes. Bonds were added by first assigning hydrogen and halogen atoms to its nearest neighbours. For all other elements, all nearest neighbours within a given cutoff were assigned as bonded neighbours to the respective central atom. For C, N, O and S, a cutoff of 2.0 Å, 1.8 Å, 1.8 Å and 2.25 Å , respectively, was used. Element types were encoded as one-hot vectors serving as node features. Edge features were built by concatenating the node features of the binding partners. We note that no distance information or chemical concepts, such as bond types, were included in the graph except for the aforementioned assumptions regarding the extraction of bonded neighbours.

A2.3 GNN
Following the previously introduced notion, the parametrization model consists of a GNN and a readout layer or combination-rule layer. Node and edge features were initially embedded as 64-dimensional vectors. The GNN consisted of independent graph-updating layers, which were composed of two fully connected feed-forward layer with 64 units, each combined with the Mila non-linearity using β = −1 [158].
Each edge-and node update layer consisted of the following module ([64, Mila, 64, Mila]) for the IPA model and ( [64,Mila]) for the ANA model. For both models, three graph-updating layers were used. The GNN module was followed by a readout module/combination rule parametrized by two fully connected layers and an output layer with n output neurons equivalent to the number of predicted parameters [64, Mila, 64, Mila, n, Softplus + ]. For the ANA model, a small term ( = 10 −3 ) was added to the output of the Softplus activation to avoid numerical instabilities, and two independent readout modules were used. Layer weights were initialized with the method introduced by He [159].

A2.4 Model Optimization
Model parameters were optimized using ADAM and the same exponential decay schedule used for annealing with learning rates (4 · 10 −4 , 4 · 10 −6 ) [165]. Gradients were clipped by their global norm with a clip norm of 1 [166].

A2.5 Condensed-Phase Simulations
Condensed-phase simulations as well as the evaluation of the electrostatic potential for intermolecular potentials of crystals were performed under periodic boundary conditions using OpenMM (7.7) [95]. Bonded terms were parametrized with OpenFF 2.0 [109] as our model does currently not provide bonded terms. The C6-C9 potential was implemented using the CustomNonbondedForce class in OpenMM.
For the IPA d+c model, monopoles were predicted for each conformation of an ensemble generated with the ETKDG conformation generator [107] implemented in the RDKit [108]. A RMS pruning threshold of 0.1 Å was used, and up to 32 conformations were generated. Monopoles were predicted for each conformation using our previously introduced equivariant GNN [63]. Fixed partial charges were then obtained by averaging over all monopoles obtained for the conformational ensemble. Charges remained fixed during the simulation, and the same charges were used for the condensed-phase as well as the vacuum simulations. The 1,4-electrostatic interactions were scaled with the same factor as in OpenFF. The 1,4-Lennard-Jones interactions were described with the C 6 -C 1 2 parameters from OpenFF 2.0, and a scaling factor of 0.5 was used. Long-range electrostatics beyond the cutoff were treated with the smooth particle-mesh-Ewald (PME) method [94]. C 6 -C 9 terms were included up to a distance of 10 Å. No long-range correction or shifting function was used for the van der Waals interactions. Initial configurations were generated using packmol [167] and conformations generated with ETKDG [107], [108]. The number of molecules was chosen such that a cubic box with side lengths 50 Å at the experimental density would be filled. Configurations were sampled from an NPT ensemble at 298.15 K using a Langevin integrator [168] with a time step of 2 fs and a collision frequency of 1 ps −1 . To maintain constant pressure, a Monte Carlo barostat [169] with a target pressure of 1 bar and a trial move every 25th step was applied to the system.. Bonds with hydrogen atoms were fixed at the equilibrium distance with the LINCS algorithm [170].
Each box was equilibrated for 2 ns followed by a 20 ns production run, with system data being saved to disk every 4 ps. System properties were averaged over the whole production run. To obtain the potential energy in the gas phase, a single molecule was simulated in vacuum using the same settings and simulation times as for the condensed-phase simulation.

A2.6 Training Data
In the following sections we provide an overview of the data sets used. Generally, only neutral molecules with elements included in {H, C, N, O, F, S, Cl} were used.

A2.6.2 Data Set II: Crystal Intermolecular Potentials
To fit intermolecular potentials in the crystalline phase, a new data set with DFT energies of molecular crystals was built. Calculations for this data set were performed on the Euler cluster of ETH Zürich. Experimental crystal structures from the CSD, which satisfy the following requirements were selected: 1. A single molecule in the asymmetric unit 2. Up to 100 atoms in the unit cell 4. Unit-cell volume up to 1600 A 3

Elements in {H, C, N, O, F, Cl, S}
A total of 35 577 structures were found to satisfy these requirements. Of which 32 811 were successfully relaxed in less than 192 CPU hours under the following settings. Structures were relaxed with the L-BFGS optimizer using the plane wave code QuantumEspresso (QE, 6.8) under the PBE functional and XDM dispersion correction [89], [99]- [104]. PBE in combination with the XDM dispersion correction was shown to perform well for molecular crystals [105]. Coordinates were relaxed under fixed lattice parameters with default QE settings, i.e. an energy convergence threshold of 10 −4 Ry and a force convergence threshold of 10 −3 a.u. The plane wave cutoff was set to 70 Ry and the charge density cutoff to 560 Ry. A uniformly spaced k-point grid was used with the number of k-points for each dimension chosen such that n k i = 25 |x i | . Cutoffs were chosen such that the lattice energy of the optimized structures of a balanced subset of the X23 database (ACETAC07, ANTCEN13, CYTSIN01, ECARBM01, HXMTAM10, TRIZIN, TROXAN11, URACIL) was converged to less than 10 −2 kJ/mol per atom [133]. Further, the publicly available projector augmented-wave (PAW) pseudopotentials (PP) [171] were used. In comparison with other publicly available PP, they provided the most accurate results during our convergence studies. The accuracy was determined as the MAE of the lattice energy with respect to experimental values. For a subset of 11 666 minimized structures, five additional geometries were generated by scaling the unit cell by factors of (0.95, 0.975, 1.0, 1.05, 1.1) without modifying the asymmetric unit. For each structure, a single-point calculation was performed using the XDM dispersion corrected B86BPBE functional [102], [172]. Previous work has shown that B86BPBE-XDM reproduces energies with very high accuracy [105]. PAW PP for the B86BPBE functional were generated with the pslibrary (1.0) [171]. The same k-points scheme was used, while the plane wave cutoff and the charge density cutoff were set to 80 Ry and 800 Ry, respectively. Monomers were calculated with a single k-point sampled at Γ in a cubic box with lengths chosen such that the minimal distance between atoms of the central cell and its periodic images was larger then 12 Å.
The intermolecular potential energies used to train the IPA d+c model were then obtained as The crystal dataset is available in the ETH research collection (https://doi.org/10.3929/ethz-b-000549359).