Absolute hydration free energy calculations support

[jmichel80]

Currently we only support relative free energy calculations but absolute solvation free energy calculations are routinely carried out in the field with well described protocols.

@lohedges @ppxasjsm how much work would be required to extend the current FreeEnergy module and support setup, simulation and processing of absolute solvation free energy calculations ?

[lohedges]

Hi @jmichel80. Having not run these simulations before I'm not sure what is involved in the setup. Do you have any good references for setting up and running absolute free energy calculations with, e.g., GROMACS?

The CCP-BioSim conference has been useful so far. One attendee is keen on running free energy perturbation simulations using NAMD and would be happy to provide a set of example input files if this is something that we wanted to add support for. (We already support basic MD protocols with NAMD.)

[chryswoods]

The collaboration with syngenta was to port their BEE code for auto-setup of absolute binding free energy calculations in gromacs into BioSimSpace. I can show you the code and ideas next week (when I am finally free from the Tier-2 bid). The main challenge is automatically defining the restraints used to hold the molecule as it vanishes. BEE has some great code to do that, and we agreed with Syngenta that we could create an open source port. The developer of BEE (Davide Sabbadine) since went to take up a postdoc position in Barcelona, so if you want to do this I can get back in touch with him and make sure this would still be ok.

[jmichel80]

Hi @lohedges ,

Absolute hydration free energy calculations are easier to setup than relative because there is no need to align and map atoms. They are also easier to setup than absolute binding free energy calculations because of the restraints. I would not jump too quickly into implementing absolute binding free energy calculations because the restraints are quite code specific.

For an absolute hydration free energy calculation, the calculation involves perturbing from a fully interacting solute to one in an ideal state that either has:

• zero partial charges and null LJ parameters

OR

• no intermolecular interaction energy

What is actually easier to implement depends on the package used.

The progression from interacting to ideal solute typically involves first removing electrostatic interactions, followed by removal of the LJ interactions (to avoid charge penetration issues).

You can see sample input files for absolute hydration free energy calculations for several simulation packages here:

https://github.com/halx/relative-solvation-inputs/tree/master/SIMS/

(package subfolder, then absolute subsubfolder)

For instance with SOMD to do an absolute hydration free energy calculation we first discharge the solute. For ethane the pert file looks like this:

version 1
 molecule LIG
     atom
         name C1
         initial_type    c3
         final_type      c3
         initial_charge -0.09435
         final_charge    0.00000
         initial_LJ      3.39967  0.10940
         final_LJ        3.39967  0.10940
     endatom
     atom
         name C2
         initial_type    c3
         final_type      c3
         initial_charge -0.09435
         final_charge    0.00000
         initial_LJ      3.39967  0.10940
         final_LJ        3.39967  0.10940
     endatom
     atom
         name H3
         initial_type    hc
         final_type      hc
         initial_charge  0.03145
         final_charge    0.00000
         initial_LJ      2.64953  0.01570
         final_LJ        2.64953  0.01570
     endatom
     atom
         name H4
         initial_type    hc
         final_type      hc
         initial_charge  0.03145
         final_charge    0.00000
         initial_LJ      2.64953  0.01570
         final_LJ        2.64953  0.01570
     endatom
     atom
         name H5
         initial_type    hc
         final_type      hc
         initial_charge  0.03145
         final_charge    0.00000
         initial_LJ      2.64953  0.01570
         final_LJ        2.64953  0.01570
     endatom
     atom
         name H6
         initial_type    hc
         final_type      hc
         initial_charge  0.03145
         final_charge    0.00000
         initial_LJ      2.64953  0.01570
         final_LJ        2.64953  0.01570
     endatom
     atom
         name H7
         initial_type    hc
         final_type      hc
         initial_charge  0.03145
         final_charge    0.00000
         initial_LJ      2.64953  0.01570
         final_LJ        2.64953  0.01570
     endatom
     atom
         name H8
         initial_type    hc
         final_type      hc
         initial_charge  0.03145
         final_charge    0.00000
         initial_LJ      2.64953  0.01570
         final_LJ        2.64953  0.01570
     endatom
end molecule

And we also setup a 'vanishing' pert file.

version 1
 molecule LIG
     atom
         name C1
         initial_type    c3
         final_type      du
         initial_charge  0.00000
         final_charge    0.00000
         initial_LJ      3.39967  0.10940
         final_LJ        0.00000  0.00000
     endatom
     atom
         name C2
         initial_type    c3
         final_type      du
         initial_charge  0.00000
         final_charge    0.00000
         initial_LJ      3.39967  0.10940
         final_LJ        0.00000  0.00000
     endatom
     atom
         name H3
         initial_type    hc
         final_type      du
         initial_charge  0.00000
         final_charge    0.00000
         initial_LJ      2.64953  0.01570
         final_LJ        0.00000  0.00000
     endatom
     atom
         name H4
         initial_type    hc
         final_type      du
         initial_charge  0.00000
         final_charge    0.00000
         initial_LJ      2.64953  0.01570
         final_LJ        0.00000  0.00000
     endatom
     atom
         name H5
         initial_type    hc
         final_type      du
         initial_charge  0.00000
         final_charge    0.00000
         initial_LJ      2.64953  0.01570
         final_LJ        0.00000  0.00000
     endatom
     atom
         name H6
         initial_type    hc
         final_type      du
         initial_charge  0.00000
         final_charge    0.00000
         initial_LJ      2.64953  0.01570
         final_LJ        0.00000  0.00000
     endatom
     atom
         name H7
         initial_type    hc
         final_type      du
         initial_charge  0.00000
         final_charge    0.00000
         initial_LJ      2.64953  0.01570
         final_LJ        0.00000  0.00000
     endatom
     atom
         name H8
         initial_type    hc
         final_type      du
         initial_charge  0.00000
         final_charge    0.00000
         initial_LJ      2.64953  0.01570
         final_LJ        0.00000  0.00000
     endatom
end molecule

There is no need to perturb bonded terms. Because this implementation turns off the intramolecular non-bonded interactions it is necessary to carry out a vacuum simulation as well. One has to be a bit careful here that the electrostatic interactions were consistently described between 'free' and 'vacuum' legs. The implementation described here uses a reaction field for the free leg, and nocutoff for the vacuum leg. This requires a trajectory postprocessing calculation to work out a correction term (implemented here )

See the README.md file here as well

For GROMACS see the instructions here

Not saying it has to be done now, but good to discuss how easily this could be implemented in the current framework.

I would implement this before absolute binding free energy calculations since the latter require all the above, plus the definition of restraints.

[lohedges]

Looking at the examples, I think absolute solvation calculations should be easy to implement in the current framework. We just update the BioSimSpace.FreeEnergy.Solvation code to also handle systems containing a single non-solute molecule that is non-perturbable, then setup and run an absolute calculation in that case (using the protocols you describe). Assuming the scripts for the correction term still work, then BioSimSpace can just run those behind the scenes when doing the analysis.