A Python tool to set up relative free energy calculations in GROMACS

Abstract

Free energy calculations based on molecular dynamics (MD) simulations have seen a tremendous growth in the last decade. However, it is still difficult and tedious to set them up in an automated manner, as the majority of the present-day MD simulation packages lack that functionality. Relative free energy calculations are a particular challenge for several reasons, including the problem of finding a common substructure and mapping the transformation to be applied. Here we present a tool, alchemical-setup.py, that automatically generates all the input files needed to perform relative solvation and binding free energy calculations with the MD package GROMACS. When combined with Lead Optimization Mapper (LOMAP; Liu et al. in J Comput Aided Mol Des 27(9):755-770, 2013), recently developed in our group, alchemical-setup.py allows fully automated setup of relative free energy calculations in GROMACS. Taking a graph of the planned calculations and a mapping, both computed by LOMAP, our tool generates the topology and coordinate files needed to perform relative free energy calculations for a given set of molecules, and provides a set of simulation input parameters. The tool was validated by performing relative hydration free energy calculations for a handful of molecules from the SAMPL4 challenge (Mobley et al. in J Comput Aided Mol Des 28(4):135-150, 2014). Good agreement with previously published results and the straightforward way in which free energy calculations can be conducted make alchemical-setup.py a promising tool for automated setup of relative solvation and binding free energy calculations.

Keywords: Automated setup; Free energy calculation; Hydration free energy; Transfer free energy.

Fig 1

alchemical-setup.py provides the missing link between planning and analysis of relative free energy calculations for GROMACS. Previous work in the group has focused on automated planning and on analysis of these calculations; alchemical-setup.py automates the setup of the calculations which can then be run to arrive at the analysis stage.

Fig 2

Two main types of relative free energy calculations. (a) Single topology, explicit intermediate. The transformation of molecule A into molecule B is split into two steps: each molecule is transformed to an intermediate whose atoms that need to be replaced are turned into dummy atoms. At this point, the two intermediates are equivalent (essentially corresponding to the CSS) except that their dummy atoms differ. However, these differences cancel when computing the total free energy change. (b) Single topology, implicit intermediate. Molecule A (left) is transformed into molecule B (right). Here, the endpoints of the transformation are identical in that the number of the atoms the molecules are comprised of is intact. As described in text, we are dealing with a single entity here, comprised of the CSS region (highlighted in blue) and peripheral atoms that are subject to get appeared and disappeared.

Fig 3

An excerpt from the topology file of the mannitol-to-tetrahydropyran interconversion displaying the content of the [atoms] directive. The entries corresponding to to-be-annihilated and to-be-appeared atoms are highlighted in green and purple, respectively, while those of the CSS are left plain. As discussed in the text, soft-core potentials must be used for both types of non-bonding interactions, as the dummy atoms are present in both end states. The chimeric molecule for this transformation is shown in Fig. 4, and the atom numbering and coloring schemes here correspond to those in Fig. 4 as well.

Fig 4

The chimeric molecule for the mannitol-to-tetrahydropyran transformation. The atoms colored as follows: atoms to be appeared are red, atoms to be disappeared are green, and the common substructure is black (the blue oxygen atoms of the common substructure are those that will be converted in the hydrogen atoms, i.e. will change their atom type). As mentioned in Section 2.2, the -CH₂- fragment (index 2 in mannitol and index 7 in tetrahydropyran) was not included in the common substructure to avoid biasing the conformation of mannitol. A section of the GROMACS topology file describing this transformation is shown in Fig. 3, and uses the same coloring scheme and atom numbering.

Fig 5

A diagram of how input for relative alchemical free energy calculations is constructed in GROMACS. Two molecules (A and B) along with corresponding coordinates and topologies are provided; a mapper is then used to determine a common substructure shared by the two molecules (as in Fig. 2a, for example); the resulting map or bijection is then used by the alchemical-setup.py algorithm to define a chimeric molecule which is a hybrid of the two (Fig. 4, for example) and generate output topology and coordinate files. Thus, alchemical-setup.py covers the final green box shown here. The main text provides a more detailed description of what is involved.

Fig 6

A table of relative hydration free energies found for four pairs of molecules depicted in the panels above. The free energies were computed with a new scheme (left column) and obtained as a difference between the absolute hydration free energies reported earlier (right column). The units are kcal/mol.