[gmx-developers] Non Equilibrium alchemy implementation in gromacs

[David van der Spoel]

Den 2018-12-14 kl. 09:35, skrev 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.

The main novelty of the paper in my opinion is the usage of statistical 
efficiency to compare different methods. It could well be that more 
efficient methods than MBAR can be devised but without efficiency 
measure, how are we to judge which method is more efficient?
I will study your paper, in particular the 2014 JCTC, for calculations 
of kOW.


> 
> 
> 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


-- 
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