GROMACS Implementation of Free Energy Calculations with Non-Pairwise Variationally Derived Intermediates

Gradients in free energies are the driving forces of physical and biochemical systems. To predict free energy differences with high accuracy, Molecular Dynamics (MD) and other methods based on atomistic Hamiltonians conduct sampling simulations in intermediate thermodynamic states that bridge the configuration space densities between two states of interest ('alchemical transformations'). For uncorrelated sampling, the recent Variationally derived Intermediates (VI) method yields optimal accuracy. The form of the VI intermediates differs fundamentally from conventional ones in that they are non-pairwise, i.e., the total force on a particle in an intermediate states cannot be split into additive contributions from the surrounding particles. In this work, we describe the implementation of VI into the widely used GROMACS MD software package (2020, version 1). Furthermore, a variant of VI is developed that avoids numerical instabilities for vanishing particles. The implementation allows the use of previous non-pairwise potential forms in the literature, which have so far not been available in GROMACS. Example cases on the calculation of solvation free energies, and accuracy assessments thereof, are provided.


Introduction
Thermodynamic systems are driven by free energy gradients. Hence, knowledge thereof is key to the molecular understanding of a wide range of biophysical and chemical processes, as well as to applications in the pharmaceutical [1,2,3] and material sciences [4,5,6]. Consequently, in silico calculations of free energies are popular in providing complementary insights to experiments or assisting the selection of chemical compounds in the early stages of drug discovery projects.
The microscopic calculation of the free energy, requires integration over all positions x of all particles in the system, where Z denotes the partition sum, β = 1/(k B T ) the thermodynamic β, k B the Boltzmann constant, T the temperature and H(x) the Hamiltonian. As an exact integration is not feasible for high-dimensional x in case of many particles, sampling based approaches such as Monte-Carlo (MC) or Molecular Dynamics (MD) simulations are commonly used. Furthermore, in practice, it oftentimes suffices to know only the free energy difference between two states, which can be calculated much more accurately. The most basic approach, rests on the Zwanzig formula [7]. The brackets A indicate an ensemble average over A is calculated. More recent methods with close relations to Eq. (3) that use samples from both A and B are the Bennett Acceptance Ratio (BAR) and multistate BAR (MBAR) method [8,9] methods.
For sampling based approaches, the accuracy of a free energy difference estimate between two states A and B generally improves when sampling is not only conducted in A and B, but also in intermediate states. Commonly, a mostly linear interpolation between the end state Hamiltonians H A (x) and H B (x) is used, where λ ∈ [0, 1] denotes the path variable. The λ dependence of the end state Hamiltonians enables the use of soft-core potentials [10,11,12] that avoid divergences in case of vanishing particle for, e.g., the calculation of solvation free energies (where the molecules "vanishes" from solution). A step-wise summation, yields the total free energy difference, where N denotes the total number of states. In the sum of Eq. (5), i = 1 corresponds to state A and i = N to state B, respectively. Alternatively, for many steps the difference can be calculated with Thermodynamic Integration (TI) [13], Importantly, advantageous definitions of intermediate states exist that go beyond the definition of Eq. (4). For example, variationally derived intermediates (VI) [14,15] minimize the mean squared error (MSE) of free energy estimates using FEP and BAR. An easily parallelizable approximation for a small number of states is where, similar to BAR, the free energy difference estimate is optimal if C ≈ ∆G. It is similar in shape to the minimum variance path (MVP) [16,17,18] for TI (2 vs 1/2 in the exponents). Enveloping Distribution Sampling (EDS) [19,20], and extensions such as Accelerated EDS [21,22] use a reference potential similar in shape to Eq. (7) to calculate the free energy difference between two or more end states.
Note a particular characteristic of the VI sequence and related methods, which is illustrated in Fig. 1: Its Hamiltonians cannot be formulated as the pair-wise sum of interaction potentials for all particles. To see this, consider the force on particle j (blue), obtained through the derivative of Eq. (7). It still depends on the full Hamiltonians of the end states. The consequence can be understood by considering a particle i (red), with λ dependent parameters, positioned at a distance r ij so large such that all direct interactions between i and j are negligible. However, when particle i changes its position with respect to its neighboring particles, the end states Hamiltonians also change, and, therefore, so does the force on particle j.
In this work, we, firstly, describe our implementation of the VI approach, and, by extension, also the MVP and basic principles of the EDS methods for two end states, into GROMACS [23,24,25]. It is among the most widely used MD software packages; however, none of the above approaches are available so far in GROMACS. Secondly, we introduce an approach to avoid singularities for vanishing particles with VI.

Avoiding End State Singularities
Interestingly, the VI sequence, Eq.(7), already exhibits soft-core characteristics for vanishing particles, as shown in Fig. (2)(a) on the example of a two-particle Lennard-Jones (LJ) potential. However, divergences can still occur when configurations from the decoupled states are evaluated at foreign states, i.e., the ones that no sampling is conducted in, but that the Hamiltonian is evaluated at such as, e.g., state B in Eq. (3). Furthermore, when two particles start to overlap, very small changes in their separation r lead to large changes in force, which causes instabilities due to finite integration steps.
To avoid these divergences, a dependence of the end state Hamiltonians on λ analogous to common soft-core potentials [10] is introduced, i.e., , respectively. For two particle i and j with distance r ij , the Coulomb Figure 1: Non-pairwise potentials and forces in VI intermediates. Two particles i and j (red and blue, respectively) are considered that are λ dependent, i.e., their interaction potential differs between A and B. It is assumed that direct interactions between i and j in both A and B are negligible. If particle i changes its position, then H A (x) and / or H B (x) change accordingly, and so does H V I (x, λ). Due to the form of the VI sequence, the derivative, and therefore, the force on particle j changes.
and Lennard-Jones interactions in state A and B are calculated based on the modified distances r A and r B , respectively, that are defined as where α and p are soft-core parameters to be specified by the user, and σ ij the Lennard-Jones parameter in the coupled state. For a system of two Lennard-Jones particles, Fig. (2) shows the resulting VI states without (a) and with (b) the use of λ dependent end states. As can be seen, the transition to the overlap region becomes markedly smoother.
Secondly, for increasingly complex molecules, the likelihood of barriers between the relevant parts of configuration space of the end states rises. Aside of additional techniques such as replica exchange, or meta-dynamics, the factor 2 in the exponent can be replaced by a user specific smoothing factor s introduced in the EDS [19,20] method. In the limit of small s, a series expansion of the exponential terms yields the conventional pathway, i.e., Eq. (4). The modified VI sequence thus reads as The force on particle i, depends on the derivatives ∂H A (x, λ)/∂λ and ∂H B (x, λ)/∂λ in the end states. Equation (13) is used for TI.
Due to the dependence of Eq. (10) on C, where the accuracy is optimal if C ≈ ∆G AB , the free energy difference has to be determined in an iterative process, where C n denotes the free energy guess at iteration step n. The free energy difference ∆G AB is obtained from simulations between state A and B , where the latter denotes the end state shifted by the constant C, i.e., that is governed by H B (x, λ) = H B (x, λ) − C. The difference ∆G AB converges to zero, such that the desired quantity ∆G AB = ∆G AB + C n ≈ C n at the end of the iteration process.

Program Structure and Usage
The end states Hamiltonians, can be split into the λ-dependent energy contributions H λ A (x, λ) and H λ B (x, λ), respectively, and the common contributions summarized by H c (x) that are equal in both end states, such as water-water interactions. To calculate H A (x, λ) and H B (x, λ), GROMACS only evaluates the λ-dependent contributions separately for the end states, whereas H c (x) is calculated only once. Note that, due to the λ dependence of the end states, H λ A (x, λ) and H λ B (x, λ) differ for different intermediates for α > 0.
The same holds for the VI sequence, Eq. (10). Inserting Eqs. 15 and 16, yields where H λ V I (x, λ) is described by Eq. (10), where the end states Hamiltonians H A (x, λ) and H A (x, λ) have been replaced by the parts H λ A (x, λ) and H λ B (x, λ), respectively, that only sum over λ-dependent interactions. The same principle applies to the calculation of the forces and λ-derivatives. Therefore, the computational effort of VI is very close to the using conventional intermediates.
However, in the current GROMACS implementation structure, all force and energy contributions from different interaction types are interpolated between the end states right after they have been calculated, i.e., the overall calculation has the form, Whereas this has the least memory requirement, for VI, the full Hamiltonians and forces in the end states need to be known before the individual forces can be calculated. Therefore, the end states Hamiltonians and forces are stored separately. After all λ-dependent contributions have been collected, first the Hamiltonian and subsequently the forces are calculated.
The implementation was built based on the GROMACS 2020 version 1 (forked on October 19th, 2019 from the master branch of the developer's repository). VI can be used with the new following entries in the mdp (i.e., input parameter) file: By nature of Eq. 10, the transformation only takes place along a single λ variable, to be specified by the mdp parameter fep-lambdas. As such, it is not possible to decouple several interactions simultaneously with different λ spacing for each type. It is, of course, possible to decouple electrostatic and LJ interactions in a sequence, that can be defined via coul-lambdas and vdw-lambdas, respectively, whereas the other is set to either zero (full interaction) or one (no interaction) for all intermediate states.

Example and test cases
When VI is switched off, all interactions are calculated as in Eqs. (18), (19) and (13). To test that VI collects all contributions correctly, for the following options in the mdp file, and likewise, for the forces and λ derivatives. Setting the seed to a fixed value such as, ld−s e e d = 1 it can be validated that all energies required for the free energy calculation that are stored in the dhdl.xvg file match between the implementation of the VI and the conventional sequence.

Equilibrium States
As an example case, the solvation free energy of nitrocyclohexane in water was calculated (structure shown in Fig. 3). The topologies of the solvation toolkit package [26] created with the Generalized AMBER Force Field [27] were used. Upon energy minimization, 2 ns NVT (constant volume and temperature) and 4 ns NPT (constant pressure and temperature) equilibration were conducted, followed by 100 ns production runs.
To asses whether the VI implementation yields accurate results consistent with the ones from conventional intermediates, first, through extensive sampling with 101 states (i.e., λ steps of 0.01), a reference value value of (9.85 ± 0.02) kJ/mol was obtained. It can be divided into (10.46 ± 0.01) kJ/mol electrostatic, and (-0.61 ± 0.02) kJ/mol LJ contributions. Next, a set of simulations with 5 states, i.e., λ steps of 0.25, were conducted.
The distribution of the free energy estimates between the different states is shown for Coulomb and LJ interactions in Fig. (4) and differs considerably between the two methods. The bars denote the free energy difference between the states denoted at the bottom. Again, A represents the coupled, and B the decoupled state. The plots shown for VI were created based on the runs where C was set to the respective reference value, and, as such, sum up to about zero. When decoupling Coulomb interactions with a conventional linear interpolation method, shown in panel (a), the largest differences between the states occur in the first steps and gradually decreases. For VI (b), the free energy path along the intermediates has be become very small (note the differing unis on the axis). In contrast, for LJ interactions, the differences for VI (d) become larger than for the linear interpolation (c). The reason is, most likely, that the differences in the contributions from the attractive and the repulsive part of the LJ potential don't cancel for all intermediates.
To compare the accuracy of both methods, Fig. 5 shows the MSEs with total simulation time, distributed equally over all five states. The MSEs were obtained by dividing the trajectories of the production runs into smaller ones, and comparing the resulting free energy difference to the reference value. For VI, two different smoothing values were considered (blue and green lines), as well as an exact initial estimate (solid line) and one that is 1 kJ/mol too low (dashed lines).
For electrostatic interactions, the MSEs in Fig. 5(a) are significantly better for VI with s = 2 and an estimate close to the exact one than the MSE obtained with linear intermediates, thereby validating the result of Ref. 14. However, in this case the MSEs are quite sensitive to the initial guess. For Lennard-Jones interactions, Fig. 5(b), VI and linear intermediates yield similar MSEs, but the VI estimates are less sensitive to the initial guess. In  both cases, the MSEs corresponding to VI with a smoothing factor of 0.1 are close to the linear ones and insensitive to the initial guess for most of the trajectory lengths in Fig. 5. As such, it is advantageous to start the iteration process with a smaller smoothing factor that is gradually increased with an improved estimate for C.

Summary
We have implemented the VI sequence of states into the GROMACS MD software package. For Coulomb interactions, our implementations yields significantly smaller MSEs and, in this sense, higher accuracy as compared to linearly interpolated intermediates. This results requires a sufficiently accurate initial estimate, which for the test cases presented here requires only a few percent of the overall simulation time. Furthermore, using the λ dependence of the end states added to VI, for LJ interactions, similar MSEs as for conventional soft-core approaches are achieved. Given the many stepwise improvements that eventually led to the accuracy of current softcore protocols, the fact the VI approach achieves similar accuracy already in the first attempt suggests that future refinements, e.g., of the lambda dependency on the end states, will push the accuracy even further.

Code and Data Availability
The source code is available at https:// www.mpibpc.mpg.de/ gromacs-vi-extension or https:// gitlab.gwdg.de/ martin.reinhardt/ gromacs-vi-extension. Documentation, topologies and input parameter files of the above test cases are also available on the website and the repository. In the gitlab repository, all changes with respect to the official underlying GROMACS code can be retraced.
As installation is identical to that of GROMACS 2020, refer to http:// manual.gromacs.org/ documentation/ 2020/ install-guide/ index.html for detailed instructions.

Acknowledgments
The authors thank Dr. Carsten Kutzner for help, discussions and advice on GROMACS code development.