[gmx-developers] Non Equilibrium alchemy implementation in gromacs

[Piero Procacci]

I raed with interest your JCP work on free energy o hydration.
However, before coming to the rather apodyctic conclusion such as
"traditional Multistate Bennett Acceptance Ratio (MBAR) is more
efficient than non-equilibrium methods such as Jarzynsky and Crooks
for hydration free energies", maybe it would be better to test the
method on a more challenging example than the one you have used. I
agree in fact, that for simple rigid molecules, with inexistent or
little conformational disorder, the HFE can be computed more
efficiently (i.e. using less computational resources) with standard
FEP.  Actually even FEP is unnecessary for most of the molecules in your
data set.  As you certainly knwow, the HFE of ethanol, for example,
can be recovered in a *single* slow growth simulation of about 2
ns. In figure 8 (left) of the 2014 JCTC paper 10.1021/ct500142c, the
2.6 ns growth and annhilation of ethanol are depicted, evidencing that
this time for this system basically corresponds to a reversible
tranformation. This would be true for 90\% of the molecules of your
data set. So you could have invested your 10+10=20 ns (!) of
equilibrium simulation of the end points in few back and forth slow
(for these system in this "fast" solvent) growth/annhiliation to get
good HFE with a good confidence interval, overperforming both FEP and
NE. Things, of course, get much tougher than in this nearly ideal
systems where nither FEP nor NE method are necessary. See for example
the growth an annhilation of the same etOH but in 1-octanol on the
right of Fig 8 in 10.1021/ct500142c. If the ligand has a complex
conformational space and/or the environment reorganized much less
rapidly than water, the NE method (coupled with h-REM) can be the only
viable alternative. Of course the NE method does not come for free in
terms of computer resources, but it can be parallelized ideally with
no loss of efficiency no matter how large is ncores.


P.S. Does the official version of gromacs already support NE simulation
in a single parallel job or you somehow patched up your 50 NE
trajectories using the current FEP implementation in discrete lambda
steps?


Piero

Il Venerdì 14/12/2018 07:39 David van der Spoel ha scritto:
> Den 2018-12-13 kl. 14:31, skrev Piero Procacci:
>> Dear colleagues
>> 
>> I would like to implement non equilibrium work (NEW) alchemical
>> methods for solvation and binding free energy calculations in gromacs
>> The idea is to transpose the ORAC hybrid OpenMP/MPI implementation for
>> NUMA architectures (which I wrote) to gromacs.  A NEW calculation is
>> done in two steps:
> 
> For your information, I am not familiar with your method but we have
> worked with some other NE free energy methods and found that
> "traditional" Multistate Bennett Acceptance Ratio (MBAR) is more
> efficient than non-equilibrium methods such as Jarzynsky and Crooks
> for hydration free energies.
> 
> https://aip.scitation.org/doi/10.1063/1.5041835
> 
> Cheers, David.
>> 
>> i) a H-REM simulation where the initial thermodynamic state is 
>> canonically
>> sampled (e.g. fully coupled solute/ligand or fully decoupled
>> solute/ligand)
>> 
>> ii) starting from a representative sample of these canonical 
>> configurations, a
>> swarm of fast annihilation or growth NE trajectories (typically few
>> hundreds) are launched in parallel on the MPI layer. Each trajectory
>> can be n-threaded on the OpenMP layer.
>> 
>> As I have no experience (as contributor) in the gromacs code, before
>> embarking in this task I would like to get some info regarding
>> the current HREM and alchemical implementation in gromacs.
>> 
>> REM question: how does gromacs implement H-REM? namely, are 
>> topology/parameter
>> arrays exchanged or the replica exchange concerns (like in ORAC) only 
>> specific
>> scaling factors (e.g .factors for torsional scaling, 14 scaling, intra
>> scaling, inter scaling etc. etc.) with invariant topology throughout
>> the GE?
>> 
>> ALCHEMY question: when electrostatic of a solute molecule are turned
>> off no longer interacting with the solvent, are the *intramolecular*
>> electrostatic interactions also scaled?  namely, if (e.g.)
>> lambda_q=0.5, do two atomic charges on the solute interact with a 0.25 
>> scaling
>> factor?  And if so, how does this get corrected when calcuting 
>> solvation
>> energies?
>> 
>> Thank you in advance for your help
>> 
>> Piero Procacci
> 
> 
> --
> David van der Spoel, Ph.D., Professor of Biology
> Head of Department, Cell & Molecular Biology, Uppsala University.
> Box 596, SE-75124 Uppsala, Sweden. Phone: +46184714205.
> http://www.icm.uu.se