Machine Learning in QM/MM Molecular Dynamics Simulations of Condensed-Phase Systems

Introduction

Classical fixed-charge force fields (FF) are readily used to perform molecular dynamics (MD) simulations of condensed-phase systems 1, 2. They consist of a relatively small number of parameters, such as bond-stretching, bond-angle bending, and torsional dihedral terms, partial charges and Lennard-Jones parameters. They are partially fitted against experimental values such as the density, heat of vaporisation, and solvation free energy.3, 4 FF are the gold standard for simulations over long time scales of systems for which long-range interactions are essential. In classical simulations, averaged properties are computed by neglecting electron rearrangement. The parameterization of FF requires the availability of sufficient experimental data. On the other hand, quantum mechanics/molecular mechanics (QM/MM)5, 6, 7, 8 MD simulations can provide a valuable alternative to classical FF simulations, when changes in the electronic structure are important or if reliable force-field parameters are not available.

In the QM/MM scheme, the QM zone simulated with density functional theory (DFT) or ab initio principles is placed into a classical environment (MM zone). This approach permits the simulation of the electronic structure of small systems in more realistic surroundings. Most crucial in QM/MM simulations is the description of the interaction between the QM and MM zones. This interaction can be based either on (i) mechanical constraints (i.e. “mechanical embedding" scheme), or (ii) electronic perturbations (i.e. “electrostatic embedding”). In mechanical embedding, the particles in the QM zone are assigned partial charges, which then interact with the MM zone on a classical level. Thus, the MM particles do not interact with the electron density of the QM solute but rather with point charges on a classical level. Mechanical embedding favours efficient implementations, however, a major drawback is that a set of classical parameters have to be determined for the QM zone, which is not trivial and sometimes not possible with sufficient accuracy.9 In electrostatic embedding, the MM environment is incorporated into the Hamiltonian operator of the QM zone as electron operators. Thus, the electron density of the QM zone is perturbed by the MM environment, and the MM zone in turn interacts with the perturbed electronic structure of the QM solute. In other words, the MM particles feel a force from the MM-polarized QM solute. The advantages of electrostatic embedding are that no partial charges have to be assigned for the computation of the interactions between the MM and the QM zones, and that the description of the interaction is physically better motivated compared to mechanical embedding. A well known limitation is the neglect of polarization of the MM environment – unless polarizable FF are used10. Furthermore, computation is more expensive compared to mechanical embedding, and it is not clear whether the partial charges of the MM zone are suitable for inclusion as one-electron Hamiltionans in the QM zone.9, 11 Generally, an electrostatic embedding scheme is to be preferred as it has been shown to be more accurate, and it reflects thus the current standard protocol for QM/MM simulations.9, 7, 12, 13, 11

In QM/MM, the description of the QM zone becomes the computational bottle-neck as it requires an expensive self-consistent field (SCF) procedure, and explicit treatment of all valence electrons14, 15. To partially circumvent these issues, semi empirical methods12, 16 can be used to describe the QM zone17. This extends the accessible time scales but also reduces the accuracy.12 An alternative is to employ machine-learned (ML) potentials.

In recent years, a lot of research effort has been invested into the development of ML models trained on the potential energy surface (PES) of QM systems18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30. This enables ML-MD simulations at an accuracy level close to that of the electronic structure method chosen to generate the training set. The costs of the resulting ML-MD simulations are reduced drastically compared to a normal DFT or ab initio MD simulation as a SCF procedure is no longer necessary and the valence electrons do not have to be treated explicitly. A large amount of approaches have been reported in the literature to achieve this task. These approaches differ in:

Use case: e.g. gas-phase19, 20 or periodic box18, small compounds20 or larger compounds,29 system specific27, 28 or system unspecific29

Descriptor used to encode the chemical structure: e.g. Coulomb matrices31, distance matrices27, 28, low-order polynomials,19, 20 or so-called symmetry functions18, 23, 29

Target output: e.g. energies18, 19 or gradients27, 28, 23

Scope: e.g. global27, 28 or local18, 21, 22, 29, 23

ML method: e.g. neural networks18, 19 or Kernel methods27, 28, 32

Applying a restriction of the final form of the PES32

Convolution of the chemical graph: Using a set of hard-coded functions,18, 27, 19 or using other techniques such as graph (convolutional) architectures33

Other ML models do not target the PES but rather the electron density34 or the wave function itself.35 However, these will not be subject of the current study. All of the approaches listed above partially generate the correct parity behaviour of the system. In other words, swapping two elements of the same type and rotating or translating the system leaves the output unchanged or changes it equivariantly. The usage of the ML models heavily depends on the specific use case. In the following, selected approaches are discussed in more detail.

Permutational-invariant polynomial neural networks (PIP)19 are suitable for accurately describing reactions of small systems in the gas phase (less than 5 atoms). Here, low-order polynomials, which depend on the internuclear distance of the compounds, are used as a descriptor for rigorously introducing the correct parity behaviour for the system studied. The neural networks are trained on the total energy of the system and are global. Symmetric gradient domain machine learning (sGDML)27, 28 targets systems up to 30 atoms in the gas phase, which do not undergo a change in topology. In this case, the compound is encoded as an inverse distance matrix. The latter is fed into a Kernel ridge regression model and trained in the gradient domain, i.e. the target output is the gradient of the system. The energy can be obtained by integrating the forces. The definition of the kernel used ensures that the integrated forces give rise to a continuously-differentiable PES. The models can be employed for spectroscopic studies28.

Larger compounds and box simulations can be tackled with high-dimensional neural network potentials (HDNNP)18, 21, 22. The environment of each atom is encoded using symmetry functions. Simple examples of symmetry functions are Gaussians that depend on the interatomic distance, or trigonometric functions. They map the input to a higher dimension, which is subsequently summed up. The summation operator and the mapping to higher dimensions using symmetry functions introduces the correct parity behaviour,36 while still allowing the HDNNP to learn the PES accurately. Using symmetry functions, Smith et al.29 showed that a universal force-field (ANI-1) can be constructed, which is capable of describing unseen compounds containing the same elements. The model was trained in a system unspecific manner on gas-phase data. The models can be employed in the drug-discovery process.37 For periodic-box simulations, HDNNP can be trained system-specific on previously sampled data points.18, 21, 22 The sampling of the data points can be done on a lower level of theory, as long as the trajectory samples all structures important for learning the PES. The data points are then re-computed on the desired level, and the HDNNP is trained to interpolate between the training data points. These models employ a cutoff for the symmetry functions and are thus local. They are trained on the energy and use optionally a gradient correction.22 HDNNP based on symmetry functions motivated other groups to apply a similar approach for different target outputs such as gradients, or with different ML models such as Kernel methods.23. The uniqueness of the descriptor has been discussed in theoretical work.38

Another important contribution to the field is the concept of $\Delta$-learning.24 Here, the ML model is trained to reproduce the difference between two QM calculations, an expensive, higher-level and a cheaper, lower-level method. Different ML methods and descriptors can be used. Subsequently, the output on the higher level can be partially recovered by performing the calculation with the cheap method and applying the trained ML model. A good example of the $\Delta$-learning scheme was proposed by Shen et al.,39 where a modified HDNNP was trained on the energy difference between a semi-empirical and an expensive QM method. It was used to reweight the free-energy profiles of reactions obtained with a QM/MM approach. More recently, the same authors introduced a model, which can be used to perform MD simulations. The approach requires as input for the ML model a reaction coordinate and the partial charges from the lower level method.40

In our opinion and that of others,30 HDNNP have been proven to be currently the most suitable method for performing system-specific periodic-box simulations. However, their usage in the simulation of condensed-phase systems with biological relevance is still hampered because such systems typically involve:41, 11, 26, 30, 32

Large number of element types: HDNNP input descriptors scale exponentially with the number of element types.

Many rotatable bonds and/or no symmetry, large system size: The quality of HDNNP model depends on a densely sampled training set. The number of different possible system configuration scales exponentially with the number of rotatable bonds, which renders sufficient sampling for the training set increasingly challenging. In principle, if the ML model is able to learn the underlying physics, it could extrapolate to unseen system configurations. However, to the best of our knowledge, this was only partially achieved yet. The models, which achieve this task partially, employ a small cutoff for the many-body (0.32 nm) and two-body terms (0.52 nm).29, 37

Important long-range interactions: HDNNP make use of a locality ansatz42, which fails to describe long-range electrostatics. Long-range interactions can partially be introduced via a charge partition scheme and a charge summation scheme (e.g. Ewald summation).22 In QM/MM, this would correspond to a mechanical embedding scheme.43 ML approaches have been proposed to specifically target this issue,42, 38 but these are to the best of our knowledge not yet applicable to simulations of large condensed-phase systems.

Large cutoff: Classical force fields typically work with a cutoff of 1.0 - 1.4 nm for pairwise nonbonded interactions.3, 1, 2 Such large cutoffs are necessary to achieve the desired accuracy for bulk properties such as heat of vaporization or solvation free energy.43

Long time scales: Although a prediction with the ML model is faster than a QM calculation, a small integration step (i.e. 0.5 fs) is still required for simulations. Furthermore, obtaining the gradients from a ML potential is usually significantly slower than calculating the gradients with a classical force field.37

In addition, HDNNP and other ML models are generally plagued by two issues: (i) Generating training sets can imply a significant time investment, and (ii) frequently used fully connected neural networks are prone to overfitting44, especially if a high number of weights is necessary to converge to chemical accuracy.45

A promising alternative for condensed-phase systems is therefore the combination of HDNNP and classical FF in a QM/MM-type approach, where the HDNNP is used for simulating the QM particles and to compute the interactions between the MM and the QM zone with electrostatic embedding. By using HDNNP for the QM zone, longer time scales will become accessible than with standard QM/MM simulations. Such a hybrid approach requires the training data set to be generated with a QM/MM scheme. In other words, the MM environment is incorporated as one-electron Hamiltonians in the QM reference calculation. With such an approach, the complexity of the generation of the training set is reduced drastically as not all atoms are treated on a QM level. Furthermore, it permits the use of a larger cutoff. Including long-range interactions directly in the ML model is therefore possible.

In this work, we assess the concept of (QM)ML/MM MD simulations for condensed-phase systems. To enable electrostatic embedding, we integrate the MM environment directly as an additional element type in the HDNNP. The symmetry functions including MM particles are weighted with their respective partial charge. Two-body and three-body symmetry functions are included for the description of the MM environment. The symmetry functions used follow the proposal of Smith et al.29, 37 and are adapted for larger cutoffs. In the first part, we explore the effect of different parameters on the accuracy of (QM)ML/MM calculations of small molecules in water. The parameters include the application of gradient correction, the inclusion of three-body terms, and the cutoff distance for the descriptors. In the second part, we extend the methodology to a $\Delta$-learning scheme24 with a semi-empirical method. Semi-empirical methods are known to recover long-range interactions well11, an issue ML models with a locality ansatz usually struggle with. The approach is validated using different molecular systems. Finally, MD simulations are performed for retinoic acid in water (50 atoms in the QM zone and approximately 2500 classical partial charges), and the interaction of SAM with cytosin (63 atoms in the QM zone and approximately 3000 partial charges). The results are compared to standard QM/MM simulations.

Theory

In non-relativistic Born-Oppenheimer approximated,46 time-independent wave function mechanics47, the electronic energy including nuclear repulsion $E_{\text{QM}}$ can be computed by solving the Eigen function,48

where $Z_{j}$ is the charge of the nuclei $j$. The indices run over the total number of electrons $N_{\text{el}}$ and the total number of nuclei $N_{\text{QM}}$. The quantum Born-Oppenheimer approximated potential energy $E_{\text{QM}}$ is obtained by computing the expectation value,

According to Newton’s second law,49 the gradients on the nuclei of the system can be computed by taking the derivative with respect to the coordinates $\vec{R}$,47

In Eq. (4), $\vec{F}_{i}$ is the force of particle $i$, $m_{i}$ its mass, $\vec{v}_{i}$ its velocity, and $t$ is the time. The system can then be propagated in time $t$ using an integration algorithm such as the leap-frog algorithm.50 Propagating the system using the quantum Born-Oppenheimer approximated potential-energy surface is computationally unfeasible for large systems. Hence, approximated solutions to Eqs. (1) and (2) were developed, which can be grouped into four classes:51 (1) Ab initio methods use the correct Hamiltonian and an approximated wave function $\psi(\vec{r})$. (2) Density functional theory52 (DFT) computes the electron density $\rho(\vec{r})$, which is used to calculate the energy of the molecular system, i.e. $E[\rho({\vec{r}})]$. Modern DFT approaches achieve an accuracy comparable to $ab$ $initio$ methods (or above), while still being able to simulate larger systems. However, DFT requires a treatment of all valence electrons and an SCF procedure, which makes it unfeasible to perform long MD simulations of large systems. (3) Semi-empirical methods use an approximated Hamiltonian $\hat{H}$ and correct for the error made by introducing a set of fitted parameters. Most prominent example is the PM7 method17. Another attractive alternative is density functional tight binding (DFTB),53 which can be considered a semi-empirical variant of DFT. In DFTB, the electron density $\rho$ is expanded in a series around a reference density $\rho_{0}$ and truncated. The error made is corrected with fitted parameters, similarly to other semi-empirical methods. (4) The last class includes empirical models with effective parameters such as classical force fields.1, 16, 2

1.2 Classical Force Fields

In classical fixed-charge MD simulations, the potential energy of the system is calculated with a force field,1, 16, 2 which is the sum of bonded and nonbonded interactions terms,

where $E^{\text{bond}}$ is the contribution of all covalent bonds, $E^{\text{angle}}$ that of all covalent angles, and $E^{\text{dihedral}}$ that of all covalent dihedrals. The nonbonded terms consist of electrostatic ($E^{\text{el}}$) and van der Waals ($E^{\text{vdW}}$) interactions. In many force fields, the van der Waals interactions are only calculated up to a cutoff and neglected afterwards (straight truncation).1 To indicate that only short-range van der Waals interactions are included, we will use the notation $E^{\text{vdW,SR}}$ in the following. The parameters of these interactions terms are typically fitted to reproduce experimental and/or QM reference data. Eq. (5) is inexpensive to solve, and depends only on the number of nuclei (atoms). The inclusion of experimental data in the fitting procedure allows to incorporate long range effects in an implicit manner. Thus, they are capable of reproducing bulk properties, such as the solvation free energy accurately (see e.g. Refs. 54, 55, 56). However, fixed-charge force fields treat electronic effects only in an averaged field, and can thus not be used to simulate chemical reactions. The simplicity of Eq. (5) allows to use relatively large cutoffs, and methods are available to include long-range interactions beyond a given cutoff.

1.3 QM/MM Scheme

QM/MM is a hybrid approach that combines a QM subsystem with an MM environment,5, 6, 7, 8 thus striking a balance between cost and accuracy. The challenge is to describe the interactions between the parts appropriately. Here, the question is how to obtain the energy $E_{\text{QM/MM}}$ of the system, which is now a combination of the QM subsystem $E_{\text{QM}}$ and the MM surrounding $E_{\text{MM}}$. A first approach to compute the total energy of the system is a subtractive scheme,7

where $E_{\text{QM}}(\vec{R})$ is the energy of the QM subsystem, $E_{\text{MM}}(\vec{R})$ the energy of the complete system calculated with the classical force field, and $E_{\text{MM}}(\vec{R}_{\text{QM}})$ the energy of the QM zone calculated with the force field. Note the distinction between $\vec{R}$, which refers to all nuclei in the system, and $\vec{R}_{\text{QM}}$, and $\vec{R}_{\text{MM}}$, which refer to the nuclei treated quantum-mechanically (QM zone) and the atoms treated classically (MM zone), respectively. The more prominent alternative is the additive scheme,7 which is also used in this work, i.e.,

where $N_{\text{MM}}$ is the number of partial charges, and $q_{j}$ the partial charge of MM atom $j$. The QM subsystem is thus directly influenced by the MM partial charges, while the latter particles feel a force of the perturbed QM zone. $\hat{H}^{\text{SR}}_{\text{QM/MM}}$ deals with the short-range van der Waals interaction, and is treated classically. It reads as follows,

where $\epsilon_{ij}$ and $\sigma_{ij}$ are fitted parameters. Thus, in electrostatic embedding Eq. (7) becomes,

Typically, not all MM partial charges are included in the summation in Eq. (8), but only those within a cutoff radius $R_{c}$.7 This leads to a PES, which is non-continuous at the cutoff. This issue can be (partially) resolved by adaptive resolution schemes.57 It has been found that the cutoff radius $R_{c}$ needs to be chosen relatively large (i.e. 1.4 nm) to converge to accurate results.11

2 Machine Learning of Potential-Energy Surfaces

Feed-forward neural networks (FNN) consist of a number of sequentially chained (non)-linear functions $\sigma^{l}$, weight matrices $W^{l}$ and biases $b^{l}$, where $l$ is the number of layers. A FNN can be defined recursively as,58

where $x$ is the input vector. The weight matrices and biases correspond to the parameters $P$. In addition, there are user-defined hyperparameters such as the number of layers $l$, the number of neurons per layer, and the type of functions $\sigma^{l}$.

The FNN architecture does not contain any information about the parity of the system. In chemistry, this is an issue because the rotation and translation of the system or swapping of atoms of the same type need to leave the energy of the system unchanged. It is possible to augment the training data set in order to force the FNN to learn the correct symmetry. However, this increases the training set by multiple orders of magnitude. Behler et al.18 identified this issue and developed a FNN architecture and an input-coordinate transformation $h$, which guarantees the correct parity behaviour, the so-called high-dimensional neural network potentials (HDNNP). Instead of using only one FNN, each element type in the system is described by its own FNN. Each FNN is called an NNP because it describes the atomic energy of the respective atom type. The total energy of the system $E$ is thus separated, i.e., $E=\sum_{i}^{N_{\text{H}}}E_{i}^{\text{H}}+\sum_{i}^{N_{\text{C}}}E_{i}^{\text{C}}+...$, where $N_{\text{H}}$ is the number of hydrogen atoms, $N_{\text{C}}$ the number of carbon atoms, and so on. Note that this separation ansatz is an approximation, it has no quantum-mechanical foundation. It guarantees that the total energy of a system is invariant with respect to swapping the positions of two atoms of the same type. However, rotation and translation still change the total energy. To address this issue, Behler et al. used the concept of symmetry functions.18 The Cartesian coordinates are transformed using a set of symmetry functions, which are invariant with respect to rotation and translation,18, 29

Applying symmetry functions to periodic boundary conditions requires a cutoff. For this, a cutoff function $f_{C}$ can be added with the cutoff parameters $R_{\text{sym,Rad}}$ (for radial symmetry functions) and $R_{\text{sym,Ang}}$ (for angular symmetry functions), i.e.,

with the cutoff function $f_{C}$ defined as,

Note that the symmetry functions $S$ take over the role of the function $h$ introduced in the previous section. To mimic long-range interactions, a second ML model can be trained on QM-derived partial charges and multipole moments to be used in an Ewald summation scheme. For covalently bonded or metallic systems, the cutoff can be chosen relatively small because the contributions of the long-range interactions are minor. However, as discussed in the Introduction, much larger cutoffs are required for condensed-phase systems.

HDNNP present a powerful technique as these models can learn the complete Born-Oppenheimer approximated PES (given an appropriate training set), i.e.,

Applying $\Delta$-learning and using the definition of $\Delta E$ in Eq. (11), we can write

However, HDNNP are limited in the following ways: First, the number of angle terms of the symmetry functions scale exponentially with the number of element types as all possible combinations must be considered. For example, a system consisting of hydrogen and carbon atoms involves five unique angle combinations. Introducing an additional element type (e.g. oxygen) adds five more unique combinations. This issue can be resolved partially by introducing weighted symmetry functions. Instead of introducing a new element type, the symmetry functions are weighted with respect to their atomic number and/or charge.26 Other issues with HDNNP are that the full correct permuting behaviour is only partially recovered,63 and that the descriptors used are non-unique.42 However, the latter issue is in our opinion more of theoretical nature, and should not limit the usage of HDNNPs in practical applications. Lastly, one has to keep in mind that HDNNPs do not resolve the inherent general limitations of ML approaches. For example, the requirement for large training sets, danger of overfitting, and limited extrapolation capabilities apply also to HDNNPs.

3 (QM)ML/MM MD Simulations Using HDNNP

In the original work, HDNNP are used to learn the complete Born-Oppenheimer approximated PES. In this work, we propose to combine the HDNNP with classical force fields in a QM/MM-type approach to reduce the complexity of learning. Hence, instead of Eq. 18, our target is

or in the case that $\Delta$-learning is applied it is,

$E_{\text{MM}}(\vec{R}_{\text{MM}})$ and $E^{\text{vdW,SR}}_{\text{QM--MM}}(\vec{R})$ are provided by the classical force field. Thus, it is expected that the complexity of the model can be reduced drastically. Furthermore, long-range interactions are now included via point charges, allowing to use a larger cutoffs in the construction of the HDNNP training set. For this approach, it is necessary to incorporate the MM environment into the descriptor of the ML model. The MM environment can be encoded in the HDNNP by introducing the MM particles as an additional element type, which is weighted with the respective partial charge. For the QM particles, we obtain

MM particles are introduced as a new element type by

In Eqs. (24-26), the summation operator iterates over all MM particles $N_{\text{MM}}$ and QM particles $N(t)_{\text{QM}}$ of type $t$. The function $Z(i)$ returns the partial charge of the MM particle. For SPC/E water, we have $Z(\text{O})=-0.8476$ and $Z(\text{H})=0.4238$. $R_{\text{sym},\text{Rad}}^{\text{QM-QM}}$ and $R_{\text{sym},\text{Rad}}^{\text{QM-MM}}$ are cutoff parameters for the radial symmetry functions, the first one influencing the QM-QM interaction and the second one targeting the inclusion of the MM particles.

Methods

Two types of systems were investigated: (i) single molecules in water, including benzene, uracil and retionic acid, and (ii) chemical reactions in water (constrained close to the transition state), namely the SN2 reaction of CH3Cl with Cl- and the reaction of S-adenosylmethionate (SAM) with cytosine. The systems were chosen to represent different system sizes, use cases, and difficulties. The systems were either treated in a static or dynamic manner. Static refers to a retrospective evaluation of sampled configurations, whereas dynamic indicates that an actual MD simulation was performed using the gradients provided by the ML model.

2 General Computational Details

All MD simulations were performed using the GROMOS software package64, 65 interfaced to DFTB+/19.253 and ORCA/4.2.014. The functionals used in this study are the gradient corrected Becke-Perdew functional (BP86), Head-Gordon’s range-separated functional $\omega$B97X-D3 66, and Grimme’s double hybrid functional B2-PLYP67. The basis set used is Ahlrichs and Weigend’s def2-TZVP68.

All QM calculations used the resolution of identity69, Weigend’s auxiliary basis70 functions, and Grimme’s dispersion correction with Becke-Johnson damping71, 72 (except for $\omega$B97X-D3). DFT calculations, which require the computation of a Hartree-Fock73 reference wave function, were further accelerated using the so-called COSX74 approximation. The ORCA computations used TightSCF convergence criteria and the integration Grid5. The grid was changed to Grid6 in the final iteration. Otherwise, standard parameters were used. DFTB computations were performed with Grimme’s D3 dispersion correction with Becke-Johnson damping.71, 72

All point charges within the cutoff radius $R_{c}$ were included in an electrostatic embedding scheme in the QM computation. Selected solvent atoms beyond the cutoff were included to avoid bond-breaking and creation of artificial charges. The convergence criteria was set to $10^{-8}$ eV for all systems. Structures, which did not converge in the maximum number of steps, were discarded. MD simulations were performed with a convergence criteria of $10^{-6}$ eV and a Broyden mixer75 with a mixing parameter of $\alpha=0.3$ to ensure numerical stability.

Newton’s equations of motion were integrated with a time step of 0.5 fs. The temperature was kept constant in the MD simulations with the Nosé-Hover chain76, 77 thermostat with a coupling constant of 0.1 ps and two baths (one coupled to the internal motion and the rotation of the solute, and the other one to the solvent). A weak-coupling78 barostat was used for constant pressure simulations with a coupling constant of 0.5 ps and a isothermal compressibility of $4.575\cdot 10^{-4}$ (kJ mol-1 nm-3)-1. Long-range electrostatic interactions beyond the cutoff of 1.4 nm were included using a reaction-field method.79 Note that the reaction field acts only on the MM particles. The water model used in this study was SPC/E80. All bonds between MM particles were constrained using the SHAKE algorithm81 and a relative tolerance of 10-4. The motion of the center of mass was removed every 1000 steps. A topology file of the QM solute is required as input in GROMOS. The Lennard-Jones parameters are needed to evaluate $E_{\text{QM--MM}}^{\text{vdW,SR}}$, since this interaction between the QM zone and MM zone is treated on the classical level. All other parameters are discarded for the QM solute. The topology files were obtained from the ATB82 server.

The neural networks were implemented using Tensorflow/keras.83 We exported the computational graphs trained in Python to the C++ GROMOS code.

3 Initial Structures

Initial structures of the individual solutes were generated using the ATB82 server. Initial structures for the chemical reactions were generated by placing the reactants manually close to the transition state. During sampling, a biasing potential was applied of the form,

where $k$ is the force constant, and $R_{ij}$, $R_{kl}$ are the distances between the atoms $i$ and $j$, and $k$ and $l$, respectively. $R_{0}$ is the ideal distance, and $d\in\{-1,+1\}$ determines, whether the distances are subtracted or added. For both systems, we set $d=-1$ and $R_{0}=0.00$. For the CH3Cl/Cl- system, the indexing is shown in the following,

The structure was minimized using a force constant for the biasing potential $k$ = 2’000 kJ mol-1 nm-2.

For the SAM/cytosine system, atom $i$ was the carbon atom participating in the reaction, atom $j$ the sulfur atom, atom $k$ was the same carbon atom as atom $i$, and atom $l$ corresponds to the carbon atom in $\alpha$-position to the amine group.

For all systems, the solute(s) was solvated in a periodic box of SPC/E water using the GROMOS++ package of programs,84 with a minimum solute-wall distance of 1.6 nm and a minimum solute-solvent distance of 0.25 nm. The size of the simulation boxes were between 3.3 nm and 5.0 nm. All systems were initially relaxed at 0 K using a gradient descent algorithm. A configuration was considered to be converged when the predicted change in energy was less than 0.1 kJ mol-1 averaged over all particles. The thermostat temperature was set to 298 K or 400 K with a coupling constant of 0.1 ps.

4 Sampling of Configurations

All systems were simulated for 20’000 steps at $T=400$ K, $p=1$ bar. The first 10’000 steps were discarded as equilibration. The remaining 10’000 frames were split into 70/20/10 training/validation/test sets. Thus, the last 1’000 steps were considered as the external test set. The performed splitting mimics the envisioned use case and gives a validation and test set, which do not resemble the training set as much as if a random partitioning would have been performed.

For benzene in water, a cutoff for the partial charges in the QM/MM Hamiltonian of $R_{c}=0.6$ nm was used. For all other systems, the cutoff for the partial charges was set to $R_{c}=1.4$ nm. For CH3Cl/Cl-1 and SAM, a biasing potential with $k$ = 2’000 kJ mol-1 nm-1 was employed. Note that for the training of the HDNNP, the biasing potential is removed. Commonly used sampling strategies such as normal-mode sampling are difficult to employ for these systems as they do not account for the solvent environment.

5 MD Simulations Using the Fitted ML Models

For the large systems, (QM)ML/MM MD simulations were performed using the fitted ML models. The time step was set for both system to 0.5 fs, the temperature was set to $T=298$ K, the pressure was set to 1 bar. For the SAM/cytosine system, the biasing force constant was set to $k$ = 2’000 kJ mol-1 nm-2.

6 Neural Networks

We implemented the HDNNP in Tensorflow/keras with a float64 precision. Each HDNNP used at most two hidden layer and thus followed the recommendation in Ref. 61. The activation function mila85 was used with a coefficient $\beta=-0.25$, because it is continuously-differentiable and is known to outperform commonly used activation functions such as tanh. The ML models were trained using the Adam optimizer86 with varying learning rates. The number of neurons per hidden layer ranged for each NNP from 10 to 80 neurons. As a loss function, we used

where $N_{\text{MM}}$ is the number of MM particles, $N_{\text{QM}}$ the number of QM particles, and $\omega_{0}$ and $\omega_{1}$ are weight parameters for the gradient contribution. If not noted otherwise, $\omega_{0}=1$ and $\omega_{1}=1$. Note that $L$ is a non-gradient corrected loss function and $L^{\prime}$ is a gradient corrected loss function. All models were initially trained for 50 steps with a learning rate of $3.5\cdot 10^{-4}$. The SAM/cytosine system was trained for 20 epochs. We monitored the validation loss during the training process. After training, the model with the lowest loss on the validation set was recovered.

7 Symmetry Functions

The symmetry functions were extracted from Refs. 29, 37. In contrast to Smith et al., we did not apply any prefactors. In addition, we allowed for a longer cutoff and extended the spacing of the radial functions. Table 1 lists the parameters for the symmetry functions in Eqs. (13) - (26), which encode the environment. The cutoffs for the radial symmetry functions were varied in this study and are provided in the Results and Discussion section for all models.

For the MM particles, we used either the same symmetry functions as for the QM subsystems (two-body terms and three-body terms), or omitted the three-body terms (only two-body terms). All symmetry functions that include MM particles were weighted with the corresponding partial charge. Note that symmetry functions could be varied to optimize the results. However, Smith et al.29, 37 have shown that their set of symmetry functions covers important local interactions and torsional space in organic compounds, while still being computationally efficient and less prone to overfitting.

In the following, we will refer to three cutoffs: $R_{c}$ is the cutoff chosen in the QM/MM calculation to generate the data sets. The summation in Eq. (8) goes over all partial charges within $R_{c}$. Similarly, the pairwise van der Waals interactions in Eq. (9) are calculated within $R_{c}$. $R_{\text{sym,Rad}}^{\text{QM-QM}}$ is the cutoff used for the radial symmetry functions encoding the QM environment (see Eq. (22)), whereas $R_{\text{sym,Rad}}^{\text{QM-MM}}$ is the cutoff used for the radial symmetry functions encoding the MM environment (see Eq. (24)). $R_{\text{sym,Rad}}^{\text{QM-QM}}$ and $R_{\text{sym,Rad}}^{\text{QM-MM}}$ correspond to the cutoffs used in the HDNNP.

Results and Discussion

The loss function of neural networks has multiple local minima. The training procedure usually finds a configuration of weight parameters close to a local minimum, but not necessarily close to the global minimum. To assess the impact of gradients correction, many-body terms, and the cutoff size, we chose test systems with a conformationally rigid solute in water where the conformational space can be sampled completely. Three different solutes were investigated, which differ in the polarity: (i) benzene (apolar), uracil (polar), and the transition state of CH3Cl/Cl- (charged). The trajectories were split into 70/20/10 training/validation/test sets (see Methods). As discussed in the Theory section, the ML models learn the quantity $E_{\text{QM}}(\vec{R}_{\text{QM}})+E^{\text{el}}_{\text{QM--MM}}(\vec{R})$. The performance of the ML models are thus assessed based on the error on this quantity (energy), and on the respective derivatives (forces), i.e. the gradients on the QM particles

and the gradients on the MM particles resulting from the QM zone,

The first test system consists of benzene in water. Benzene is an apolar solute. Thus, long range interactions do not contribute significantly to the potential energy compared to polar solutes. Point charges further away than 0.6 nm from the solute atom were not included in the evaluation of Eq. (8). The test system is suitable to study the effect of including gradient correction and many-body terms (Eqs. (24) - (26)) in the description of the PES for the MM part. Note that the descriptor always contains many-body information for the QM particles (Eq. (16)). Table 2 summarizes the settings of four ML models used for the comparison. All models use the parameters proposed by Ref. 29. However, we account for the interactions of the QM subsystem with the MM particles by introducing them as a new element type (see Theory section). Model 1 (M1) is the baseline model, including only two-body information and no gradient correction. M2 includes three-body terms for the MM point charges. M3 only includes gradient information. M4 additionally weights the gradients of the MM particles in the loss function. Note that the cutoffs $R_{\text{sym,Rad}}^{\text{QM-QM}}$ and $R_{\text{sym,Rad}}^{\text{QM-MM}}$ for M1-M4 are larger than $R_{c}$. This means that the models “know” about all particles contributing to the potential energy.

We first analyzed the model performance as a function of the number of neurons in the hidden layers for M1 and M2. The accuracy on the energies of the validation/test set is far above chemical accuracy for both models, and also the performance on the QM gradients (mid panel) and the MM gradients (right panel) is poor. Thus, including three-body terms in the descriptor for the MM particles (M2) does not improve the accuracy. Interestingly, the performance decreases even further with an increasing number of neurons in the hidden layers. This is in line with the observation in ML that large numbers of weight parameters should be avoided, when regularization techniques are not applied.

M3 includes gradient information of the QM and MM particles in the training procedure, and thus regularizes the weights of the models. The accuracy on the energies of the validation set is now clearly below chemical accuracy, and also the performance on the gradients is significantly improved. We observe that the performance is stabilized, meaning that a larger number of neurons per layer slightly improves the performance instead of worsen it.

The influence of regularizing the models via gradients can be further highlighted by varying the parameter $\omega_{1}$. In Eq. (30), $\omega_{1}$ scales the contribution of the error on the MM gradients. The gradients on the MM particles are much smaller than the QM gradients, since the quantity only incorporates the electrostatic interaction with the QM zone, and not the interaction of the MM particles with themselves. The latter interactions as well as the short-range van der Waals interactions between the QM and MM particles are modulated with the classical force field. Thus, the contribution of the MM particles to the loss function is substantially less than the contribution of the energy and forces of the QM particles if $\omega_{1}=1$ is used. The performance on the training set measured by the loss function is thus significantly increased at the cost of a worse description of the interaction with the MM particles. This behaviour is clearly non-physical and leads to a HDNNP, which performs worse on the validation (and test) set. Scaling the MM gradients via $\omega_{1}$ solves this issue. Figure 2 shows the performance of M4 on the energies and gradients as a function of the parameter $\omega_{1}$. By increasing $\omega_{1}$, the prediction accuracy can be further improved.

1.2 ΔΔ\Delta-Learning

As ML models use a locality ansatz, it is difficult for them to describe long-range interactions. Semi-empirical methods, on the other hand, are known to be able to account for long-range interactions. Thus, a $\Delta$-learning24 scheme based on a semi-empirical method might be a valuable approach. In such a scheme, the ML model introduces a correction to the semi-empirical method to recover the accuracy of the higher-level method. Although this decreases naturally the computational efficiency of the production run, since an additional computation is required at each time step $t$, we will demonstrate in the following that the performance is increased significantly, and that the number of weight parameters can be reduced. As a result, the amount of training data points needed is reduced, the training procedure is accelerated, and the extrapolation capabilities of the model are possibly increased.

Table 3 summarizes the settings. The $\Delta$-learning model M5 uses the same settings as M1, i.e. only two-body terms for the QM - MM interactions, no gradient correction ($\omega_{0}$=0, $\omega_{1}$=0), and a cutoff for the symmetry functions of 1.0 nm. M6 uses the same settings as M3, i.e. gradient information are included in the training procedure ($\omega_{0}$=1, $\omega_{1}$=1). M7 uses stronger gradient regularization ($\omega_{0}=1$, $\omega_{1}=200$), and is thus comparable to M4. In Figure 3, the $\Delta$-learning models are compared to M3.

The results show that using gradient correction is crucial for performance as M3 (gradient correction but no $\Delta$-learning) clearly outperforms M5 ($\Delta$-learning but no gradient correction). At the same time, the prediction accuracy can be substantially increased by using a $\Delta$-learning in combination with gradient correction (M6). However, while the MAE with M6 is below 4.18 kJ mol-1 for the validation set, it is above it for the test set. By increasing the regularization (M7, $\omega_{1}=200$), the performance can be further improved for the test set. M7 presents thus the overall best performing model. In general, it is important to carefully monitor the performance of neural networks as a function of their parameters.

Using $\Delta$-learning and gradient correction improved the accuracy of all three properties, i.e. energies, QM gradients, and MM gradients, significantly already at a small number of neurons per hidden layer. This indicates that only a reduced number of training data points is necessary to achieve the accuracy of the reference method (see discussion below). It also suggests that the extrapolation capabilities of M6/M7 might be better than those of the previous models. This makes the model particularly interesting for use cases, where the ensemble of conformations/configurations of the system in the training set is not complete. To the best of our knowledge, the lowest level of theory available for such a $\Delta$-learning scheme is DFTB. Other semi-empirical methods do not incorporate the partial charges in an SCF manner in the Hamiltonian. Rather, they describe the interaction based on parameters. This means that the gradients on the MM particles usually do not correlate with those computed by the reference method (DFT).

1.3 Size of the Training Set

A major issue for learning of complex systems is that not all relevant configurations can be enumerated. Thus, ML approaches that can be trained on smaller data sets are especially valuable. As shown above, $\Delta$-learning requires less weight parameters to achieve chemical accuracy on the validation and test set of benzene in water, which indicates that less training points are needed compared to the ML models trained directly on the full energies and gradients. To quantify this observation, the M4 ($\omega_{1}=200$, 80 neurons per layer) and M7 ($\omega_{1}=200$, 10 neurons per layer) models were trained on only 10% of the data set. This was done by changing the training/validation/test splitting from 70/20/10 to 10/20/10 (i.e. discarding the last 60% of the training set). More precisely, the ML models with a reduced data set set are fitted on the first 1000 data points of the training set, while the validation and test set remain the same. The results are presented in Table 4. It is striking that for M4 the prediction accuracy on the energies and gradients of the same validation and test sets drops by up to 50%, when trained only on 10% of the data set. For M7, on the other hand, the error on the energies remains relatively constant when only trained on 10% of the data set. The performance on the gradients worsens as well but to a smaller extent.

1.4 Cutoff Size

Fitting a ML model for the benzene in water system did not require the inclusion of long-range interactions as the partial charges beyond 0.6 nm were truncated. To assess the effect of the cutoff size and thus the improved accounting for long-range interactions, the test system of uracil in water is used in the following. The cutoff $R_{c}$ is increased up to 1.4 nm, which corresponds to 1500-2000 MM partial charges in the cutoff sphere. A cutoff of 1.4 nm is used in the GROMOS force field,3 of which the SPC/E water model is part of. Table 5 lists the settings M8-M13. All models use gradient correction, M11-M13 follow the $\Delta$-learning scheme with DFTB. M11 serves as a “worst-case” model for $\Delta$-learning, where the MM particles are simulated solely with the semi-empirical method, and the ML model introduces a local correction to the QM subsystem.

The results with the models M8-M13 and the DFTB baseline for uracil in water are shown in Figure 4 and summarized in Table 6. We will first discuss the models without $\Delta$-learning, M8-M10. For these, long-range interactions need to be learned completely by the ML model. M8 has no information about the surrounding environment, and thus the performance on the energies and the QM gradients is poor. Increasing the cutoff of the symmetry function for the description of the MM environment to 0.52 nm (M9) increases the performance significantly. The MAE on the energies on the validation set is reduced by nearly a factor 4, and the description of the gradients also becomes better. Interestingly, increasing the cutoff further worsens the results. The reason is likely that the training data set is relatively sparse (7000 data points). Increasing the cutoff in the symmetry functions also increases the number of weight parameters, which need to be fitted. This complicates the training procedure, making it harder to find meaningful minima. The situation is different for the $\Delta$-learning approach (M11-M13). Here, the semi-empirical method can be used to incorporate important long-range interactions, while the ML models serves as a local correction. A striking example provides M11. M11 like M8 has no information about the surrounding partial charges. Nevertheless, the model almost reaches chemical accuracy on the validation set. Compared to M8, the gradients on the QM particles are improved by a factor three, and the gradients on the partial charges by a factor 7. Including information about the partial charges in the surrounding up to 0.52 nm improves the results further. Again, going from a cutoff of 0.52 nm to 1.4 nm worsens the results, likely due to the same reason. However, given a sufficiently large training set (larger than in the current cases), a $R_{\text{sym,Rad}}>0.52$ nm may be desirable.

Lastly, we compare DFTB to the ML models. All ML models outperform DFTB on predicting the QM gradients. However, DFTB outperforms the ML models without $\Delta$-learning on the MM gradients. This highlights again the importance of applying $\Delta$-learning for systems, where long-range interactions are important and only sparse data sets are available. Further, it is interesting to observe that DFTB outperforms M8, which has no information about the surrounding partial charges, on the QM energies. This finding suggests that well-fitted semi-empirical methods may be preferred over ML models, if the latter do not account appropriately for the surrounding partial charges.

To assess the transferability of the conclusions from the previous test systems, we analyzed the performance of the settings M11-M13 on a third test system, the (close to) transition state of the SN2 reaction between CH3Cl and Cl-. For this system, larger changes in the electronic structure within the configurational ensemble are expected. Furthermore, as it is a charged system, it is interesting to see, how well the long-range interactions are described by the $\Delta$-learning models. The results for M11-M13 for the CH3Cl/Cl- test system are summarized in Table 7. In general, we see the same behaviour as for uracil in water, i.e. a cutoff of 0.52 nm performs best.

2 Application of ML Models in MD Simulations

So far we have assessed the performance of the ML models in terms of MAE (and RMSE) on a validation and test set. However, these results can be to some degree misleading because rare outliers get averaged. Such rare outliers might be less tolerable in the context of an actual MD simulations, where the results of the next step depend directly on the results of the previous step.

In order to judge the usefulness of HDNNP models in practice, it is important to test their performance in actual MD simulations. For this purpose, we have chosen two test systems with larger conformational flexibility: (i) retinoic acid (a vitamin A derivative) in water, and (ii) the (close to) transition state of the reaction of SAM with cytosine in water (Figure 5). Sampling to generate the training and validation sets was performed at elevated temperature (400 K) with a time step of 0.5 fs for both systems. Instead of an external test set, MD simulations were carried out at $T=298$ K. Our main objective in this section is to show that given a sparse training set the proposed $\Delta$-learning scheme is able to produce stable trajectories, which describe the PES sufficiently well even for unseen data points. We therefore use as above only 10’000 data points for training, even though the complexity of the systems is much higher compared to the previously studied ones. For real case studies, e.g. the computation of free-energy profiles, a more sophisticated sampling procedure for the training set is recommended.

The ML models investigated employ the settings M12, i.e. $\Delta$-learning scheme with DFTB as baseline, 10 neurons per hidden layer, only two-body terms for the QM - MM interactions, gradient correction (with $\omega_{0}=1$ and $\omega_{1}=200$), and $R_{\text{sym,Rad}}$ of 0.52 nm. The QM/MM reference calculations were computed with a cutoff $R_{c}$ of 1.40 nm.

First, we assessed the performance of the ML model on a validation set, as done in the previous sections. Table 8 summarizes the results in terms of MAE and RMSE, Figure 6 shows them graphically. Using the $\Delta$-learning model, the performance is improved compared to the DFTB baseline for all three properties, including the MM gradients.

Next, we tested the trained ML model by performing an (QM)ML/MM MD simulation for 5000 steps (top panel in Figure 7). In order to compare the ML-corrected energies with those of the DFT reference and the DFTB baseline simulation, single point calculations were performed for each configuration in the ML + DFTB trajectory. As shown in Figure 7, the energies of ML + DFTB agree better with the DFT reference than the DFTB baseline. Although the latter performs relatively well, one should keep in mind that the deviations from the reference method accumulate along an MD trajectory. Thus, even relatively small deviations might result in largely different trajectories when starting from the same coordinates. We note that the MAE on the trajectory is with 5.8 kJ mol-1 larger than the MAE of 3.9 kJ mol-1 and 4.4 kJ mol-1 for the training and validation sets, respectively. The reason for this is that the starting point of the trajectory is the relaxed structure, a configuration never seen by the HDNNP. To achieve a higher accuracy in practical MD simulations for complex molecules such as retinoic acid, more training data points will be needed for the HDNNP. To assess the stability of the ML model, we extended the simulation to 200’000 steps (bottom panel in Figure 7).

2.2 SAM/Cytosine Transition State in Water

In biochemistry, S-adenosylmethionate (SAM) is a co-factor for the transfer of a methyl group by enzymes (methyltransferases). Here, we investigated the transition state of the chemical reaction between SAM and cytosine. Again, we first assessed the performance on a validation set as done above (Table 9 and Figure 8). The $\Delta$-learning model clearly outperforms the DFTB baseline.

Next, we performed a (QM)ML/MM MD simulation for 2000 steps using the fitted model (top panel in Figure 9). As for the test system with retinoic acid in water, the energies with DFTB + ML agreed well with the DFT reference, and outperformed the DFTB baseline as indicated by the MAE and RMSE over the trajectory. To assess the stability of the model, we extended the simulation to 50’000 steps (bottom panel in Figure 9). For both test systems, it was possible to carry out a stable (QM)ML/MM MD simulation for the selected number of steps (up to 200’000 steps). Of course, the possibility cannot be excluded at this point that a configuration, which is not well represented by the ML model, might be encountered in even longer simulations.

3 General Discussion

The $\Delta$-learning scheme requires an explicit computation with the lower-level method at each time step. In this work, we have used DFTB as the lower-level method. For the SAM/cytosine system with 63 QM atoms and up to 3500 MM particles in the cutoff sphere, one evaluation with DFTB required less than a second (single core), while the corresponding DFT computation takes 60 to 80 minutes (on 4 cores, CPU time), which is more than three orders of magnitude longer. Thus, the additional time needed for the DFTB calculation at each step of the $\Delta$-learning scheme is less expensive than e.g. the pairlist creation in the MD engine. An obvious limitation of the $\Delta$-learning approach are configurations, for which the DFTB PES is substantially different from the reference DFT PES. In such cases, the DFTB computation will converge very slowly (or not at all), limiting the usage. However, we did not encounter such an issue in the presented MD simulations.

Summary and Conclusion

In this work, we investigated the use of HDNNP in (QM)ML/MM MD simulations of condensed-phase systems with DFT (or $ab$ $initio$) accuracy for the QM subsystem. Standard QM/MM MD simulations at this level of accuracy are very expensive and only applicable to small systems. Using semi-empirical methods to describe the QM parts is much faster but also less accurate. In the HDNNP, the MM partial charges up to a certain cutoff are integrated as additional element type.

We assessed the influence of different parameters on the prediction accuracy of energies, gradients of the QM particles, and gradients of the MM particles. First, we used a simple test system of an apolar, rigid solute in water. Three-body terms to describe the QM - MM interactions could be neglected, however, the inclusion of a gradient correction in the training procedure of the ML model was crucial to obtain sufficient accuracy on the gradients and energies. In general, we found that strong regularization was necessary during the training procedure to predict these gradients with decent accuracy. This may no longer apply when larger training sets are used. As an alternative, we employed a $\Delta$-learning scheme, where the lower-level method is DFTB. Semi-empirical methods are known to handle long-range interactions well. For uracil in water, we showed that $\Delta$-learning allows one to decrease the cutoff of the symmetry functions substantially, while still outperforming classical HDNNP. The removal of symmetry functions reduces the number of weight parameters and the complexity of the training procedure. Smaller models usually require less training data points and display improved extrapolation capabilities. This may be a key factor for large biomolecular systems, where it is not feasible to enumerate all conformations/configurations for the training set.

The final $\Delta$-learning model was tested further by performing actual (QM)ML/MM MD simulations of relatively large systems, i.e. retinoic acid in water and the transition state of the chemical reaction between SAM and cytosine. Stable trajectories were obtained with an accuracy close to the DFT reference. The performance can be further boosted by a more sophisticated sampling strategy for the training and validation set. We envision that the results and findings presented in this study will enable the use of (QM)ML/MM MD simulations in practical applications.

Acknowledgements

The authors thank Patrick Bleiziffer and Annick Renevey for helpful discussions. S.R. gratefully acknowledges financial support by the Swiss National Science Foundation (Grant Number 200021-178762) and by ETH Zurich (ETH-34 17-2).

References