Machine learning of solvent effects on molecular spectra and reactions

I Introduction

The presence of a solvent can dramatically change molecular properties as well as the outcome of reactions . Hence, a profound understanding of the interactions between molecules and their environments is tantamount not only for rationalizing experimental results, but also guiding the way towards controlling, or even designing, reactions and properties . While computational chemistry has described such phenomena with reasonable success, obtaining accurate predictions which can be related to experiment is a highly non-trivial task. Due to the large number of species involved, treating a system and its environment entirely with electronic structure theory is prohibitively expensive, especially for accurate high-level methods. Approximate schemes, on the other hand, are often unable to capture important physical aspects, such as structural features of the environment in the case of continuum models or chemical reactions in the case of classical force-fields. To overcome these issues, we propose a deep neural network potential that includes the influence of external fields, to capture the interactions of the chemical system with the environment.

Recently, machine learning (ML) methods have emerged as a powerful strategy to overcome this trade-off between accuracy and efficiency inherent to computational chemistry approaches. Highly efficient ML models now provide access not only to interatomic potential energy surfaces , but also to a growing range of molecular properties . Specialized ML architectures have been developed for the prediction of vectorial and tensorial quantities , such as dipole moments , polarizabilities and non-adiabatic coupling vectors . Such advances pave the way for using ML models in practical applications, with the simulation of infrared , Raman , ultraviolet and nuclear magnetic resonance (NMR) spectra being only a few examples. At the same time, there is an ongoing effort to incorporate more physical knowledge into ML algorithms, giving rise to semi-empirical ML schemes and even models based on electron densities and wavefunctions .

Most of these approaches operate on closed systems, where the molecule is not subjected to any external environment. Christensen et al. proposed a general framework for modeling response properties with kernel methods. In this context, and in order to model electric field dependent properties, the FCHL representation was extended by an electric field model based on rough estimates of atomic partial charges. This makes it possible to predict various response properties beyond atomic forces, such as dipole moments across compositional space as well as static infrared spectra. Note that this approach inherently relies on molecular representations which are able to capture the required perturbations of the energy.

In this work, we propose the FieldSchNet deep learning framework, which models the interactions with arbitrary environments in the form of vector fields. As a consequence, our model is able to describe various implicit and explicit interactions with the molecular environment using external fields as a physically motivated interface. Moreover, the field-based formalism automatically grants access to first and higher order response functions of the potential energy, with polarizabilities and nuclear magnetic shielding tensors being only a few examples. This enables the prediction of a wide range of molecular spectra (e.g. infrared, Raman and NMR) without the need for additional, specialized ML models.

FieldSchNet can be used as a polarizable continuum model for solvation or to interact with an electrostatic field generated by explicit external charges in an QM/MM-like setup – an approach we refer to as ML/MM. As the model offers speed-ups of up to four orders of magnitude compared to the original electronic structure reference, we can employ FieldSchNet to investigate the influence of solvent effects on the Claisen rearrangement reaction of allyl-p-tolyl ether – a study out of reach for conventional high-level electronic structure or ML approaches. While these simulations would take 18 years with the original electronic structure reference, they could be performed within 5 hours with FieldSchNet.

FieldSchNet goes well beyond the scope of previous ML approaches by providing an analytic description of a chemical system in its environment. We exploit this feature by designing an external field in order to minimize the height of the reaction barrier in the above Claisen reaction. Hence, FieldSchNet constitutes a unified framework that enables not only the prediction of spectroscopic properties of molecules in solution, but even the inverse chemical design of catalytic environments.

II Results

starting from an initial embedding depending only on the respective atom type. Here, $\mathbf{w}_{i}^{l}$ is the standard SchNet interaction (Fig. 1a, left block), while the added terms $\mathbf{u}_{i}^{l}$ and $\mathbf{v}_{i}^{l}$ correspond to dipole-field and dipole-dipole interactions (Fig. 1a, right block). The SchNet interaction update then takes the form

The radial interaction functions $\mathbf{W}_{q}^{l}$ depend on the interatomic distance $r_{ij}$ and are learned from reference data. A fully-connected neural network $\operatorname{NN}$ is applied afterwards performing a non-linear transformation.

In terms of rotational symmetry, the SchNet feature refinements can be interpreted as charge-charge interactions, as the features of $\mathbf{x}_{i}$ are scalars. Hence, the internal structure of SchNet, as well as the generated representation, is invariant with respect to rotations and translations of the molecule. This is an important requirement for machine learning potentials, as the energy of an atomic system exhibits the same symmetry . However, this invariance breaks down in the presence of external fields. In this case, models also need to be constructed such that they are able to resolve rotations and translations relative to an external frame of reference.

Based on these features, FieldSchNet models the interaction between molecule and external fields $\mathbf{u}_{i}^{l}$ with the term

where $\bm{\epsilon}_{\text{ext}}(\mathbf{R}_{i})$ is the field acting on atom $i$. The dot product corresponds to the physical expression for first order approximation to the energy of a dipole in an external field. This ensures that the orientation with respect to the external field is captured correctly and reduces to a constant in the absence of a field.

In addition to the dipole-field interaction, FieldSchNet introduces an update $\mathbf{v}_{i}^{l}$ based on the interaction between neighboring dipoles

where $\mathbf{T}(\mathbf{R}_{ij})$ is a 3-by-3 Cartesian interaction tensor inspired by the classical dipole-dipole interaction

which guarantees the right geometrical behavior. Similar to the SchNet interaction in Eq. 2, $\mathbf{W}(r_{ij})_{\mu}^{l}$ is a learnable radial filter allowing the network to modulate the interaction strength.

For each external field, FieldSchNet adds a set of dipole features and expands the update in Eq. 1 by the corresponding dipole-field and dipole-dipole terms. As the scalar features $\mathbf{x}_{i}$ of the next layer depend on the dipole and field interactions of the current layer, those in turn are coupled with the preceding scalar features via the dipoles (Eq. 3). This allows FieldSchNet to capture interactions with external fields far beyond the linear regime.

As visualized in Fig. 1d at the example of an ethanol molecule, the FieldSchNet descriptor exhibits the same symmetry as the molecule (top). Upon introducing an electric field in the z-axis of the molecule, the descriptor adapts to the changed environment and the original symmetry is broken (middle). This effect can also be observed in the response of the descriptor with respect to the field (bottom). At the same time, multiple successive SchNet and dipole-dipole updates enable the efficient construction of higher-order features of the molecular structure (see Fig. 1c), going beyond the purely radial dependence of charge-like features. Although inspired by electric dipoles $\bm{\mu}$, the dipole-like representations in FieldSchNet are auxiliary constructs and may also represent other quantities, e.g. the nuclear magnetic moments $\mathbf{I}_{i}$ when modeling magnetic fields. An appropriate form of the dipoles corresponding to each field is learned purely data-driven.

Once the FieldSchNet representation $\mathbf{x}_{i}$ has been constructed, the potential energy is predicted via the atomic energy contributions typical for atomistic networks (Fig. 1b)

Since the features $\mathbf{x}_{i}$ depend on the interactions between dipoles and fields of each previous layer and the dipoles are in turn constructed from the features, the potential energy is now a function of the atomic positions, nuclear charges and all external fields as well as initial atomic dipoles. This makes it possible to obtain response properties from the energy model by taking the corresponding derivatives.

II.2 Molecular spectroscopy

FieldSchNet is particular promising for the prediction of molecular spectra, as a single model provides access to a wide range of response properties. A variety of spectra can be simulated in this manner, ranging from vibrational spectra such as infrared and Raman to nuclear magnetic resonance spectra. The availability of molecular forces makes it possible to go beyond static approximations and even derive vibrational spectra from molecular dynamics simulations . Moreover, the high computational efficiency of the machine learning model renders otherwise costly path integral molecular dynamics simulations feasible, which are able to account for nuclear quantum effects and yield predicted spectra close to experiment .

We train a single FieldSchNet model on reference data generated with the PBE0 functional for an ethanol molecule in vacuum. A combined loss function incorporates the energy ($E$), atomic forces ($\mathbf{F}$), dipole moment ($\bm{\mu}$), polarizability ($\bm{\alpha}$) and nuclear shielding tensors ($\bm{\sigma}$) as target properties (see Supplementary Text 1 for details). Excellent fits were obtained for all quantities, as can be seen based on the test accuracy reported in Supplementary Tab. 3. Energies and force predictions fall well within chemical accuracy, with mean absolute errors (MAEs) of 0.017 kcal/mol and 0.128 kcal/mol/Å. The response properties of the electric field exhibit MAEs as low as 0.04 D ($\bm{\mu}$) and 0.008 Bohr3 ($\bm{\alpha}$) compared to value ranges of 4.56 D and 13.58 Bohr3 present in the reference data. In a similar manner, FieldSchNet yields low MAEs for the shielding tensor $\bm{\sigma}$, exhibiting an error of 0.123 ppm for hydrogen (range of 29.43 ppm) and 0.194 ppm for carbon atoms (range of 153.94 ppm).

In the following, we simulate a range of spectra using the multi-property FieldSchNet model. Fig. 2a shows infrared spectra obtained from static calculations and the dipole time-autocorrelation functions simulated by molecular dynamics. In addition, a static spectrum obtained with the reference method, as well as the gas-phase experimental spectrum are provided. While the FieldSchNet model is able to reproduce the static reference almost exactly, comparison of both spectra to experiment demonstrates the drawbacks of relying on purely static calculations for the prediction of vibrational spectra in general.

The potential of the machine learning model becomes apparent when going beyond the static picture. In order to predict vibrational spectra from molecular dynamics simulations, a large number of successive computations are necessary, which quickly become prohibitive when relying on electronic structure methods. However, due to the computational efficiency of FieldSchNet, these simulations can be carried out with little effort. A single evaluation of all response properties takes 220 ms on a Nvidia Tesla P100 GPU compared to a computation time of 207 s with the original electronic structure method on a Intel Xeon E5-2690 CPU, a speedup by almost three orders of magnitude. As a consequence, a simulation which would take 240 days with conventional approaches can now be performed in approximately 6 hours. This huge step in efficiency grows even stronger for larger systems.

The benefits of performing molecular dynamics simulations can be observed in the spectrum recovered in this manner, as it exhibits a much better agreement with experiment in the low frequency regions. Nevertheless, this approach still fails to reproduce positions and intensities of the bands associated with the stretching vibrations of the C-H bonds ($\sim$ 3000 cm-1) and the O-H bond ($\sim$ 3600 cm-1). These differences are primarily due to the neglect of anharmonic and nuclear quantum effects. One way to include these effects is via path integral molecular dynamics, where multiple coupled replicas of the molecule are simulated. While this approach is computationally more demanding than classical molecular dynamics due to the additional replicas, the simulations can still be carried out efficiently with FieldSchNet. As can be seen in Fig. 2a, accounting for anharmonic effects does indeed shift the C-H stretching vibrations to the experimental wavelengths and improves the position the O-H band, yielding a predicted spectrum close to experiment. The remaining difference in the O-H band can most likely be attributed to a general limitation of the reference method.

In addition to infrared spectra, FieldSchNet enables the simulation of Raman spectra, which rely on molecular polarizabilities $\bm{\alpha}$. Since the full polarizability tensor is predicted by the model, it is possible to compute polarized as well as depolarized Raman spectra. Fig. 2b depicts both types of spectra as obtained with FieldSchNet via path integral molecular dynamics, as well as their experimental counterparts . In both cases, very good agreement with experiment is observed, with the most prominent difference being once again the vibrations in the O-H stretching regions. The quality of the predicted spectra is particularly striking, when comparing to the statically computed polarized Raman spectrum, which fails to reproduce the shapes and magnitudes of several peaks.

Beyond vibrational spectra, FieldSchNet can be used to obtain NMR chemical shifts via the nuclear shielding tensors $\bm{\sigma}_{i}$. In this manner, chemical shifts can be obtained for all NMR active isotopes, which in the case of ethanol are 1H, 13C and 17O. The predicted and reference chemical shifts for the equilibrium configuration of ethanol are provided in Fig. 2c and d. In addition, the distribution of shifts sampled during path integral molecular dynamics are shown. 17O shifts are omitted from the analysis, as only one oxygen nucleus is present. However, the shifts are still reproduced accurately and the associated error is given in Supplementary Tab. 3. FieldSchNet predictions agree closely with the reference method for the 1H and 13C isotopes. The 1H chemical shifts of the hydrogens in the CH2 and CH3 groups are close to their expected experimental values of 1.2 ppm and 3.8 ppm, respectively. The peak shifted to 1 ppm is associated with the hydrogen atom in plane with the O-H bond. The resulting band structure vanishes in the molecular dynamics simulation due to rotations of the methyl group. A major disagreement with experiment is the shift of the hydrogen in the O-H group which shows uncharacteristically low values of 1 ppm. This can be attributed to a shortcoming of the reference method, which exhibits the same behavior. The 13C shifts of the CH2 and CH3 carbon atoms agree almost perfectly with their experimental values of $\sim$60 ppm and $\sim$20 ppm.

II.3 Implicit environments with polarizable continuum models

Accounting for effects of the molecular environment and solvent effects in particular is crucial for a wide range of chemical applications. These effects can critically influence the properties of compounds and the outcome of chemical reactions. Due to the large number of species involved, treating a molecule and surrounding environment entirely with electronic structure methods is impractical. In the case of solutions, approximating the solvent by a polarizable continuum model (PCM) with specific dielectric constant $\varepsilon$ has proven as a powerful tool .

By using an expression for the external field adapted from the reaction field approach of Onsager , FieldSchNet can operate as a machine learning model for polarizable continuum solvents with an explicit dependency on $\varepsilon$ (see Sec IV.3). To study this mode of operation, we train such a polarizable continuum FieldSchNet (pc-FieldSchNet) on a reference data set composed of computations for a ethanol molecule in the gas phase, as well as ethanol ($\varepsilon=24.3$) and water ($\varepsilon=80.4$) continuum solvents. Again, the model is trained to predict the potential energy, the atomic forces and the response properties for the molecular spectra. The polarizable continuum FieldSchNet (pc-FieldSchNet) reproduces all quantities with high accuracy, comparable with that of the model for response properties in vacuum (see Supplementary Tab. 3). This demonstrates that pc-FieldSchNet is able to learn the correct dependence on the specific dielectric constant $\varepsilon$.

Supplementary Tab. 3 furthermore shows results for two benchmark datasets of methanol ($\varepsilon=32.63$) and toluene ($\varepsilon=10.3$), solvents that have not been included in training. We find that pc-FieldSchNet generalizes well to these solvents with errors comparable to the prediction of solutions used for training. The accuracy for toluene is slightly lower with the polarizability showing particularly high mean absolute errors of 0.243 Bohr3. However, this is can be attributed to an insufficient sampling of the regions of low polarity, as the most similar solvent included in training is vacuum. The methanol dataset is reproduced with high accuracy, demonstrating that the model is able to generalize across unseen continuum solvents.

II.4 Explicit environments with ML/MM

Although continuum models are powerful tools for the description of solutions, they break down in situations where direct interactions between molecule and environment need to be considered, e.g. solute-solvent hydrogen bonds. In these cases, the solvent has to be treated explicitly, e.g. using QM/MM schemes . These retain the full atomistic structure of the environment but instead treat it with more affordable classical force fields, while the molecule itself is described with electronic structure methods. By expressing the external electric field as the field generated by the point charges of the MM region, FieldSchNet can operate in a similar manner and replace the quantum region in such a simulation, essentially yielding an ML/MM approach. We use generalized atomic polar tensor charges computed with FieldSchNet to model the electrostatic interactions between both regions (see Sec IV.4 for details). In the following, we study the impact of using an ML/MM solvent model on the simulated infrared spectrum of liquid ethanol.

We train a FieldSchNet model on a set of PBE0 reference data of ethanol configurations polarized by external charge distributions that have been sampled from a classical force field. The test set performance of the resulting model is provided in Supplementary Tab. 3. We observe slightly increased errors for the ML/MM model in energy and forces compared to the vacuum and continuum models due to the more complex environment. Still, the model is able to reach high accuracy.

ML/MM molecular dynamics simulations are carried out for a single machine learning ethanol surrounded by a solvent box of 1250 ethanol molecules treated at force field level. The resulting infrared spectrum is depicted in Fig. 3, alongside spectra computed in the same manner using the vacuum model (Sec. II.2) and continuum model ($\varepsilon=24.3$), as well as an experimental spectrum of liquid ethanol. Comparing the gas phase and continuum models, we find that in this case implicit solvent effects lead to no major improvements with respect to experiment. The spectrum simulated via the ML/MM approach on the other hand, yields significantly better predictions. The low wavelength regions in particular show excellent agreement with experiment. The high wavelength regions are shifted to higher wavelengths since anharmonic effects are neglected in the classical ML/MM simulation. However, we still observe the broadening and red-shift of the O-H stretching vibration present in the experimental spectrum. This effect is caused by hydrogen bonding between the O-H groups in the machine learning region and the surrounding ethanols (see inset Fig. 3). Continuum models fail to account for these kind of interactions, as they neglect the structure of the solvent.

II.5 Modeling organic reactions

Effects of the environment and solvents play a central role in many molecular reactions. Certain arrangements and combinations of molecules in the environment can promote or inhibit reactions in a dramatic fashion. As such, a good understanding of these effects is crucial for the development of new catalysts and drugs. However, accurate computational simulations of such systems are hard to obtain since intense sampling procedures are required in order to obtain reliable free energy profiles of the studied reaction. This is further complicated by the need to account for the large number of molecules in the environment, which in many cases cannot be described by more affordable continuum models of solvation due to explicit interactions between solute and environment.

An example for such a reaction is the Claisen rearrangement of allyl-p-tolyl ether (Fig. 4a). The presence of water as a solvent accelerates this reaction by a factor of 300 compared to reaction rates in the gas-phase . Computational and experimental studies have determined that the main reason for this acceleration is explicit hydrogen bonding between the transition state and the water molecules of the solvent. These lead to a lowered barrier and promote the reaction . The combination of computational efficiency and accuracy with the ability to perform ML/MM simulations makes FieldSchNet well suited for modeling such reactions. Beyond that, the model provides access to a range of properties which can be used to characterize the different species formed during reaction.

We train two FieldSchNet models to simulate the first step in the Claisen rearrangement of allyl-p-tolyl ether. The first model is trained on PBE0 reference configurations sampled from a metadynamics trajectory of the reaction simulated at a lower level of theory. A second model is a ML/MM model based on the reference configurations determined above augmented by different MM charge distributions sampled via the TIP3P force field for water (see Sec IV.1 for details). Errors for both models are reported in Supplementary Tab. 4.

The resulting free energy barriers are depicted in Fig. 4b, along with potential energy barriers computed in vacuum using the reference method and vacuum model. The FieldSchNet ML/MM model correctly predicts a lower activation barrier (30.08 kcal/mol) for the aqueous environment compared to the gas phase reaction (33.35 kcal/mol). The overall difference in the barrier height $\Delta\Delta G=3.28$ kcal/mol is close to the experimental value of $\Delta\Delta G=4$ kcal/mol. Analyzing the configurations sampled during the ML/MM simulation, we observe the hydrogen bonding between the ether oxygen and water molecules in the solvent responsible for the acceleration of the reaction . Fig. 4c shows the radial distribution function between hydrogens in the solvent and the oxygen of the transition state. A pronounced peak at a distance of 2 Å indicates that hydrogen bonds between solvent and solute form frequently at this stage. Note that the absolute predicted height of the barriers is underestimated compared to experiment (38.4 and 34.4 kcal/mol for vacuum and water, respectively). An analysis of the gas phase reaction barrier reveals that this is due to the reference and not an artifact of the machine learning model (Fig 4b).

The NMR chemical shifts provided by FieldSchNet can be used to trace structural changes occurring during the rearrangement and allow to connect theoretical predictions to experiment. Fig. 4d depicts the 13C chemical shifts predicted for different stages of the reaction. For example, it shows how the first (C12) and third (C16) carbon atom in the allyl ether exchange chemical environments during reaction. The former moves from typical shifts of allyl ether groups ($\sim$68 ppm) to shifts associated with terminal carbons in conventional allyl groups ($\sim$130 ppm). At the same time, the C16 undergoes this process in reverse, ending at shifts of $\sim$50 ppm characteristic for this position in allyl substituents. Another change of interest are the shifts of the carbon C2 in the aromatic ring where a new bond forms. This atom starts at typical values for aromatic carbons ($\sim$120 ppm), staying there during the transitions state. Upon formation of the product, the aromaticity of the ring is lost and the shift moves to regions more indicative for carbon atoms in cyclic ketones ($\sim$25-50 ppm).

II.6 Designing molecular environments

The analytic nature of neural networks allows to establish a direct relation between molecular structure and properties which can be exploited in inverse chemical design applications . FieldSchNet is well suited for such tasks, as it provides access to a wide range of response properties as a function not only of molecular structure but also the external environment. This offers the possibility to manipulate external fields and molecular environments in order to optimize certain properties or to control reaction rates. In the following, we apply FieldSchNet to design a chemical environment promoting the Claisen rearrangement reaction studied above.

The external field in the ML/MM FieldSchNet model used to describe the reaction depends on the charge distribution of the environment. Thus, the external charges can be optimized to lower the reaction barrier by minimizing

Fig. 5a shows the evolution of the barrier during various stages of the optimization procedure. By designing an optimal environment, the activation barrier can be reduced by $\sim 25$ kcal/mol, lowering it from an initial 35 kcal/mol to 10 kcal/mol. These findings correspond to a rate acceleration by a factor of $\sim 2\times 10^{18}$ at 300 K. Fig. 5b and c show the optimized environment in presence of the transition state, visualizing regions of negative (b) and positive (c) charge. Motifs such as the strong negative density close to the oxygen atom and the neighboring carbons involved in the forming bond are in strong agreement with experimental studies, where it was found that electron donating groups in these regions promote the reaction .

We outline how such an optimized field can guide design endeavors by translating these motifs into an atomic environment of amino acids as would be present in an enzymatic active site. Negatively charged aspartic acid residues (Asp) are placed close to regions exhibiting high negative density, while regions of high positive charge are populated with lysine (Lys) molecules. Using Eq. 8, we optimize the placement of the amino acids based on the electrostatic field derived from their Hirshfeld atomic partial charges. The resulting environment is shown in Fig 5d. Finally, we recompute the reaction barrier in presence of the amino acid charge distribution with the original electronic structure reference (Fig. 5e). Although the optimal environment cannot be reconstructed perfectly due to the constraint imposed by the structure of the molecules, this relatively straightforward approach leads to a significant reduction of the barrier. The activation energy is lowered from 35 kcal/mol to 22 kcal/mol, which still corresponds to an acceleration by a factor of $\sim 3\times 10^{9}$. This demonstrates that the theoretically optimized environment Fig. 5b and d can be mimicked retaining the same trend on the barrier and reaction speed.

III Discussion

FieldSchNet enables modelling the interactions of molecules with arbitrary external fields. This offers access to a wealth of molecular response properties such as polarizabilities and nuclear shielding tensors without the need for introducing specialized models. Leveraging external fields as an interface, the model can further operate as a polarizable continuum model for solvents, as well as a surrogate for the quantum mechanics region in quantum mechanics/molecular mechanics (QM/MM) schemes, yielding a ML/MM approach. As a consequence, FieldSchNet paves the way to many promising applications out of reach for previous techniques.

Computational spectroscopy can benefit greatly from the FieldSchNet framework, as it provides efficient means for computing high quality spectra as we have demonstrated for IR and Raman spectra as well as NMR chemical shifts. Combining the latter with nuclear spin-spin coupling tensors, i.e. response properties of the magnetic moments, enables the accurate prediction of nuclear magnetic resonance spectra. In principle, all other spectroscopic quantities, derived via the response formalism, can be modeled by FieldSchNet as well.

The ability of our neural network to operate as a continuum model for solvation is not only highly attractive for applications in drug design, where efficient models of solvent effects are much sought after, but serves as a starting point for developing more powerful models, which could for example consider structural aspects of the solvent. When treating interactions with the environment explicitly, the introduced ML/MM procedure proves to be a powerful tool combining the accuracy of machine learning potentials with the even higher speed of force fields. While the presented study of solvent effects on the Claisen rearrangement reaction of allyl-p-tolyl ether would require 18 CPU years with the electronic structure reference, it was performed within 5 hours with FieldSchNet on a single GPU. This greatly expands the range of application of ML models and brings the simulation of large, biologically relevant systems, such as enzymatic reactions, within reach.

Moreover, the fully analytic nature of FieldSchNet enables inverse design as we have illustrated by minimizing the reaction barrier of the Claisen rearrangement through interactions with an optimized environment. Coupling this to a generative model of molecular structure , the charge distributions found using FieldSchNet may be populated in a fully automated fashion. Possible application of FieldSchNet to inverse design tasks include the targeted functionalization of compounds or the creation of enzyme cavities promoting reactions.

FieldSchNet constitutes a unified framework for describing reactions and spectroscopic properties in solution. Beyond that, it provides insights on how these quantities can be controlled via the molecular environment. As this opens up new avenues for designing workflows tightly integrated with experiment, we expect FieldSchNet to become a valuable tool for chemical research and discovery.

IV Methods

All electronic structure reference computations were carried out at the PBE0/def2-TZVP level of theory using the ORCA quantum chemistry package . SCF convergence was set to tight and integration grid levels of 4 and 5 were employed during SCF iterations and the final computation of properties, respectively. In the case of allyl-p-tolyl ether, the RIJK approximation was used to accelerate computations . Nuclear shielding tensors were computed with the Gauge Including Atomic Orbitals approach implemented in ORCA, while continuum solvents calculations were performed with the in package conductor-like polarizable continuum model .

The reference data for ethanol was generated by selecting 10 000 random configurations from the MD17 database and recomputing them at the above level of theory. In addition, continuum solvent calculations for the four studied continuum solvents (toluene, ethanol, methanol, water) were carried out for the structures selected in this manner. A training set for continuum models containing 30 000 ethanol configurations were constructed by merging the vacuum, ethanol and water data. Reference data for the ML/MM simulations was generated in a two-step approach. Initially, a periodic box of 1250 ethanol molecules was equilibrated with the NAMD molecular dynamics package using the CHARMM General Force Field for 1 $\mu$s. Using the native NAMD interface to ORCA, electrostatic embedding QM/MM simulations were carried out, where one of the ethanols was described at the PBE/def2-SVP level of theory. CHELPG charges were used as partial charges for the quantum regions and simulations were run for 50 ps using 0.5 fs time steps. For all simulations, temperatures were kept at 300 K using a Langevin thermostat and pressures at 1 atm using a Langevin piston barostat . From this trajectory, 30 000 QM ethanol configurations and the associated charge distributions of the environment were sampled at random and recomputed at the PBE0/def2-TZVP level.

Reference data for the allyl-p-tolyl ether Claisen rearrangement reaction in vacuum was obtained via metadynamics at the PBE/def2-SVP level of theory. The two bonds involved in the reaction (see Fig. 4) were selected as collective coordinates and Gaussians with a height of 1 kcal/mol and a width of 0.529 Å were deposited each 100 simulation steps. The system was simulated for a total of 50 ps using 0.5 fs time steps. Temperature was kept constant at 500 K by means of a Nose-Hoover chain thermostat . We then selected 61 000 configurations from this and recalculated them with the reference level of theory. Data for MM/ML simulations was generated by suspending 20334 configurations sampled during metadynamics in periodic solvation boxes with 9260 TIP3P waters . Keeping the allyl-p-tolyl ether coodinates frozen, the water box was then optimized and simulated for 50 ps with NAMD. For temperature and pressure control, the same setup as in the ethanol box was used. From each of these boxes, 3 ether configurations and associated charge distributions were drawn and recomputed at the PBE0/def2-TZVP level, yielding 61 002 reference data points.

IV.2 Training setting

IV.3 Model for continuum solvation

A FieldSchNet-based machine learning potential for continuum solvent effects can be derived by adapting the Onsager expression for the reactive field . Modeling the external electric field $\bm{\epsilon}_{\textrm{solv}}$ experienced by each atom $i$ due to a continuum solvent with dielectric constant $\varepsilon$ as

allows for a direct coupling between the molecular potential energy and the solvent. The solvent field is adapted in each layer $l$ based on the atomic dipole features $\bm{\mu}_{i}^{l}$. The term $a^{l}_{i}$ models the effective volume of the nucleus accessible due to its environment and is modeled by a neural network.

IV.4 Model for hybrid machine learning/molecular mechanics (ML/MM)

We adopt a QM/MM-like approach for FieldSchNet by coupling the classical and quantum region with electrostatic embedding, resulting in a ML/MM approach. In this case, the quantum region is polarized by the charges of the classical region, while the electrostatic energy of the classical region is in turn modified by the point charges computed for the quantum region. The influence of the external charges on the molecule can be modeled via the electrostatic field exerted by a collection of point charges

where $q$ are the external charges and $\mathbf{R}_{k}$ the associated positions. A suitable set of partial charges for the machine learning region can be obtained in the form of generalized atomic polar tensor charges , which are the response property

These charges are fully polarized charges, depending not only on the molecular structure but also the charge distribution of the environment.

IV.5 Computational Details

Unless stated otherwise, the velocity Verlet algorithm and a time step of 0.5 fs were used to integrate the equations of motion. All simulations not using NAMD were carried out with the molecular dynamics module implemented in SchNetPack .

Classical molecular dynamics simulations for ethanol in vacuum and continuum solvents were carried for 50 ps at a temperature of 300 K controlled via Nose-Hoover chain thermostat with a chain length of 3 and time constant of 100 fs. The first 10 ps of these trajectories were then discarded. Ring polymer molecular dynamics were performed for 20 ps, using a time step of 0.2 fs and a specially adapted global Nose-Hoover chain as introduced in Ref. 81 to keep the temperature at 300 K. Once again, a chain length of 3 and time constant of 100 fs were chosen for the thermostat.

Simulations for the ethanol ML/MM model were carried out using a custom interface between NAMD and our machine learning code. First, a periodic box of 1250 ethanol molecules was equilibrated with NAMD for 1 $\mu$, using the CHARMM General Force Field . Bonds to hydrogens were kept frozen with the SHAKE algorithm and a time step of 1 fs was used. One ethanol was then selected for modeling via the FieldSchNet ML/MM model and ML/MM simulations were carried out using a custom interface between NAMD and our machine learning code for a total of 50 ps. For both simulations, temperatures were kept at 300 K with a Langevin thermostat and pressures at 1 atm using Langevin piston barostat . The first 10 ps of the trajectory were discarded.

Umbrella sampling simulations for the Claisen rearrangement set up according to the following protocol. Using the difference between the bonds formed and broken as the reaction coordinate, we determined the centers for the harmonic bias potentials by choosing 50 equidistant points along the reaction coordinate. The centers ranged from values of -4.15 Å to 5.18 Å with an increment of 0.19 Å. For each center, we selected the closest lying structure in the metadynamics trajectory used for generating the reference data as a starting configuration for the umbrella sampling run. All simulations used a force constant of 112.04 kcal/mol/Å2. Umbrella sampling in vacuum was carried out for each window by first equilibrating the system for 25 ps using a Berendsen thermostat at 300 K (time constant of 100 fs) followed by 25 ps production simulation at the same temperature with a Nose-Hoover chain (chain length of 3 and time constant of 100 fs). For the ML/MM model umbrella simulations, the starting configurations were first solvated in a periodic box of 9260 water molecules treated with the TIP3P force field. Keeping the allyl-p-tolyl ether structures frozen, the water box was first minimized and the equilibrated for 200 ps to a temperature of 300 K and pressure of 1 atm with a Langevin thermostat and Langevin piston barostat using NAMD. Bonds involving water hydrogens were kept frozen with the SHAKE algorithm and a time step of 1 fs was used. Starting from the systems prepared in this manner, ML/MM simulations were performed for 25 ps using the same pressure and temperature control as above.

Free energy profiles were constructed from the umbrella sampling data using the WHAM code with convergence set to 1e-9 and a temperature of 300 K .

Infrared and polarized as well as depolarized Raman spectra were computed from the time-autocorrelation functions of the dipole moment and polarizability time derivatives according to the relations given in Ref. 47. Autocorrelation functions were computed using the Wiener-Khinchin theorem and a autocorrelation depth of 2048 fs. In order to enhance the quality of the spectra, a Hann window function and zero-padding were applied to the autocorrelation functions before computing the spectra. A laser frequency of 514 nm and temperature of 300 K were used for calculating the Raman spectra.

Data availability

All datasets used in this work will be made available on http://www.quantum-machine.org/datasets after publishing.

Code availability

All code developed in this work will be made available upon request.

Author Contributions

MG and KTS conceived the research. MG developed the method and carried out the reference computations and simulations. KTS, MG, KRM designed the experiments and analyses. MG and KTS wrote the paper. KTS, MG and KRM discussed results and contributed to the final version of the manuscript.

Competing Interests

The authors declare no competing interests.

References

Appendix A Supplementary Information: Machine learning of solvent effects on molecular spectra and reactions

The energy predicted by FieldSchNet is an analytic function of the coordinates $\mathbf{R}$, as well as the external fields and their associated atomic dipole moments. This makes it possible to access so-called response properties, which are partial derivatives of the potential energy . Assuming the presence of an external electric $\bm{\epsilon}$ field and a magnetic field $\mathbf{B}$ with its corresponding nuclear magnetic moments $\{\mathbf{I}_{i}\}$, a general response property $\bm{\Pi}$ takes the form

The power of FieldSchNet (and field-based models in general) lies in the fact, that a single energy function provides access to a wide range of quantum chemical properties in a highly systematic manner. Moreover, the expression in Eq. 14 above guarantees the correct geometric transformations of the property tensors with respect to rotations and translations of the molecule in the external field without the need of explicitly encoding the corresponding symmetries. As is the practice with molecular forces, response properties can also be incorporated during training of the FieldSchNet model by including the appropriate squared errors into the loss function

Here, the trade-offs $\eta$ weight the importance of a property in the loss and $N$ is the total number of atoms. The properties predicted by FieldSchNet according to Eq. 14 are indicated with a tilde.