Insights and Challenges in Correcting Force Field Based Solvation Free Energies Using a Neural Network Potential

We present a comprehensive study investigating the potential gain in accuracy for calculating absolute solvation free energies (ASFE) using a neural network potential to describe the intramolecular energy of the solute. We calculated the ASFE for most compounds from the FreeSolv database using the Open Force Field (OpenFF) and compared them to earlier results obtained with the CHARMM General Force Field (CGenFF). By applying a nonequilibrium (NEQ) switching approach between the molecular mechanics (MM) description (either OpenFF or CGenFF) and the neural net potential (NNP)/MM level of theory (using ANI-2x as the NNP potential), we attempted to improve the accuracy of the calculated ASFEs. The predictive performance of the results did not change when this approach was applied to all 589 small molecules in the FreeSolv database that ANI-2x can describe. When selecting a subset of 156 molecules, focusing on compounds where the force fields performed poorly, we saw a slight improvement in the root-mean-square error (RMSE) and mean absolute error (MAE). The majority of our calculations utilized unidirectional NEQ protocols based on Jarzynski’s equation. Additionally, we conducted bidirectional NEQ switching for a subset of 156 solutes. Notably, only a small fraction (10 out of 156) exhibited statistically significant discrepancies between unidirectional and bidirectional NEQ switching free energy estimates.


Additional details for the EXS and UVIE protocols Protocol UVIE-CGenFF as the classical force field ASFE calculations with transformato
4][5][6][7] The solutes were placed in cubic simulation boxes with a side length of ≥ 26 Å, which is sufficiently large to be commensurate with the default CHARMM cut-off radius of 12 Å for Lennard-Jones interactions.
Coulomb interactions were calculated using the particle-mesh Ewald (PME) method 8 (Ewald κ = 0.34Å −1 , 32 × 32 × 32 FFT grid); Lennard-Jones (LJ) interactions were calculated using a cutoff of 12 Å and the default OpenMM switching function between 10 Å and 12 Å.These steps were automated with the Python package macha, * which uses CHARMM scripts generated by CHARMM-GUI 9,10 as templates.When used to compute ASFEs, transformato provides an automated workflow to set up the serial atom insertion (SAI) approach. 11Since electrostatic and Lennard-Jones intramolecular interactions are annihilated, the gas phase correction step (see left side of Figure 1 in the main manuscript) is required to restore both intramolecular interactions.Each ASFE calculation was repeated four times with different initial random velocities.The standard deviation was calculated from these four runs.

Treatment of virtual particles on chlorine atoms
Recent versions of the CGenFF force field 5,6 use virtual particles to better describe the halogen bond, a specific non-covalent interaction between a halogen atom and another electronegative atom, driven by a localized positive electrostatic region known as the "sigma hole". 7While these virtual particles can be handled straightforwardly in the transformato workflow, they must not be present during the NEQ switching simulations from the MM to the ANI-2x level of theory.For the 24 molecules that contain a chlorine atom, we, therefore, had to ensure the closing of the alchemical cycle by introducing a suitable intermediate state.
We calculated the free energy difference between the original state (fully charged virtual particle attached to each chlorine present) and an intermediate state, where the virtual particles are present, but their charges have been transferred to the corresponding chlorine atom.In this intermediate state, the massless virtual particles do not contribute to the system's free energy.Hence, the virtual particles can be safely removed when carrying out the endstate correction from MM to ANI-2x.

Nonequilibrium switching simulations
Setup of the forward (MM→NNP/MM) switching simulations: Four independent 5 ns MD simulations were carried out, using the Langevin integrator with a target temperature of 303K and a friction coefficient of 1/ps.Pooling the trajectories of the 4 simulations of the physical endstate for each of the systems and excluding the first 30% of each run as equilibration, we collected 28,000 coordinate sets from which the forward switching protocol (MM → NNP/MM) was initialized.
Setup of the backward (NNP/MM→MM) switching simulations: At the NNP/MM endstate, a single 10 ns trajectory was generated for each system (all other simulation parameters as just described).Again, the first 30% were discarded, resulting in 7,000 coordinate frames.
Switching simulations: From the pool of initial configurations, 300 coordinate sets were selected randomly (with replacement) as the starting point of the forward and backward NEQ switching simulations, both in the gas phase and in solution.

Protocol EXS-OpenFF as the classical force field ASFE with openmmtools
Molecules were parametrized using Biosimspace 12 and the OpenFF-2.0force field. 13The solvated systems were created by solvating the ligand in a 40Å 3 cubic water box.In solution, electrostatic interactions were computed using PME 8 with a short-range cutoff of 10 Å; no cutoff was used in the gas phase.Bond lengths involving hydrogens were constrained to their parameter value in all simulations.OpenMM 8.0 14 was used to perform all simulations.Before starting an MD simulation, both the gas phase and solvated structure were minimized.Simulations were performed using the LangevinMiddleIntegrator implemented in openmmtools 15 with a target temperature of 300 K and a friction coefficient of 1/ps, using a time step of 2 fs.Each system in the gas phase was equilibrated in NVT for 1 ns.Each solvated system was equilibrated for 1 ns in the NPT ensemble with a Monte Carlo barostat at a target pressure of 1 atm.The ASFE calculation consisted of two phases: In solution, interactions were decoupled/annihilated using 26 λ-windows.First, intra-and intermolecular electrostatic interactions were scaled to zero using an equidistant λ-spacing of 0.1 (i.e., electrostatic interactions were annihilated in 11 steps).Then, the van der Waals interactions were decoupled with a non-equidistant scheduling (using 0.05 spacing until λ = 0.5, then continuing with a 0.1 spacing).Soft core potentials were used as described in Equation 13of Ref. 16.Because of the mixed annihilation/decoupling scheme, the gas phase corrections only had to restore intramolecular electrostatic interactions.Each λ-window was simulated for 10 ns.Errors were estimated from the single repeat equilibrium simulations via the bootstrapping functionality of MBAR as implemented in pymbar. 17

Nonequilibrium switching simulations
Simulations of the physical endstates were repeated, treating the solute as fully flexible (i.e., only water molecules were kept rigid by constraints).During these simulations, 5,000 frames were saved.From these coordinates, 300 were selected randomly (with replacement) and used to initialize 5 ps forward switching simulations (MM → NNP/MM).

An attempt to go beyond mechanical embedding
As described in the main manuscript, correcting hydration free energies to the NNP/MM level of theory led only to statistically insignificant changes compared to the MM results; i.e., the magnitude of the endstate correction was small in most cases.We strongly suspect that this is a consequence of the mechanical embedding 18 currently used for the NNP/MM coupling.Since the performance of ANI-2x is sufficiently fast, we attempted to go beyond mechanical embedding for a handful of compounds as follows.All calculations described below are based on protocol UVIE; the starting point are the trajectories saved during NNP/MM simulations (mechanical embedding) needed for the bidirectional calculation (Crooks' equation) of the end state corrections.Given a set of coordinates, the solute-solvent interaction energy can always be obtained as the difference Here, U(full) is the total energy, U(water) the energy of just the waters, and U(solute) the energy of the solute with the waters removed.For each frame saved during the NNP/MM simulations, the necessary energies were recomputed (1) using the NNP/MM coupling and mechanical embedding as during the production calculations, and (2) using ANI-2x for the full system.We, thus, extracted interaction energies U MM inter and U NNP inter at the MM and NNP levels of theory, respectively.We then estimated a correction based on the NNP interaction energy according to The correction ∆G inter was then added to the endstate corrected hydration free energy obtained with mechanical embedding.
The results are summarized in Table S1.While, with one exception, the magnitude of the endstate corrections resulting from mechanical embedding (column EC) was small (< 0.5 kcal/mol), that of the endstate corrections including interaction energies at the NNP level of theory is much larger.However, in most cases the resulting hydration free energy (column ASFE) is in worse agreement with experiment than the MM and the NNP/MM result obtained with mechanical embedding.resulting from replacing the MM interaction energy by the NNP interaction energy (see Eq. 2); h standard deviation of ∆U inter = U NNP inter − U MM inter ; i hydration free energy after applying the correction for the interaction energy as described.
Nevertheless, some conclusions can be drawn from the data summarized in Table S1.
First, endstate corrections obtained with solute-solvent interactions at the NNP level of

Figure S2 :Figure S3 :
FigureS2: Unidirectional NNP correction (∆G corr MM→NNP/MM ) as a function of the number of rotatable bonds for the 156 compound subset using protocol EXS.The molecules are numbered, starting with one (1) for the compound having the highest correction value.Outliers observed for both force fields (CGenFF and OpenFF) are highlighted in purple.See Figure5in the main manuscript for the analogous plot using the CGenFF force field.

Table S1
also includes the standard deviation σ(∆U inter ) of the difference in interaction energies ∆U inter = U NNP inter − U MM inter .As one sees, σ(∆U inter ) > 4 kcal/mol in all cases.In Ref. 19 it was shown that results obtained from Zwanzig's or Jarzinsky's equation become unreliable if the standard deviation of the energy differences / work values exceeds 4 k B T ≈ 2.4 kcal/mol at room temperature.This is the case for all compounds considered; hence, ∆G inter is likely subject to large systematic errors.

Table S1 :
Results obtained by correcting the solute-solvent interaction energy from MM to NNP.All free energies are in kcal/mol.ID used in FreeSolv;20,21 bcommon name; c experimental hydration free energy reported in FreeSolv; d MM hydration free energy obtained using protocol UVIE; e end state correction as obtained using protocol UVIE (mechanical embedding only); f hydration free energy after reweighting to ANI-2x/MM, mechanical embedding, protocol UVIE ; g correction ∆G inter a