michellab
Absolute hydration free energy calculations support
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 ?
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.)
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.
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.
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.