[gmx-users] Question about decomposition of Lennard-Jones contributions in free energy calculations with Gromacs

[David van der Spoel]

On 11/11/16 01:22, Joel Jose Montalvo Acosta wrote:
> Dear Dr. van der Spoel,
>
> Thank for your fast answer.
>
> I took time to think and to make tests about why is meaningful to have a
> dispersion/repulsion decomposition in an alchemical calculation. In
> general, the solvation free energy (ΔG-solv) can be decomposed in
> electrostatic (ΔG-elec) and van der Waals (ΔG-vdw) contributions, where
> this last term is also possible to decompose in dispersive (ΔG-disp) and
> repulsive (ΔG-rep) contributions, following the Weeks-Chandler-Andersen

Sorry to be harsh, but this kind of decomposition does not mean anything 
thermodynamically since the distinction between electrostatic forces and 
van der waals forces is arbitrary. It is only relevant algorithmically 
for e.g. convergence of free energy calculations. Hence, unless you are 
trying to develop smarter or more efficient sampling methods, you are 
wasting your time. Any paper that is based on the physical 
interpretation of Delta G-elec and/or Delta G-vdw should be rejected.




> (WCA) approach. Each contribution has different convergence profile in a
> FEP/TI run, thus, ΔG-elec and ΔG-disp can be obtained using a small
> number of equally-sized lambdas because the plots ⟨dUelec(λ)/dλ⟩ and
> ⟨dUdisp(λ)/dλ⟩ vs λ usually present a near linearity. However, ΔG-rep
> will need more and non-equally sized λ since this contribution is
> related to the cavity formation of the solute in the solvent. Then, the
> computational time needed to get converged ΔG-elec and ΔG-disp values
> will be less than the needed for converged ΔG-rep values.  As you
> suggested, the most efficient order to do the decomposition will be
> electrostatic, dispersion and finally repulsion, this last term will
> converge because it's only related to a repulsive potential for
> solute-solvent interactions, and for a large enough box of solvents, the
> pressure will not be significantly affected for the solute-solvent
> perturbation and the run could be done in NPT ensemble.
>
> As one example, I computed the desolvation free energy (-ΔG-solv) of
> fullerene C60 in toluene (OPLS-AA parameters from virtualchemistry.org
> <http://virtualchemistry.org>) with gromacs using 25 λ (production of 5
> ns per λ) and the same lambdas' distribution used in your paper J. Chem.
> Inf. Model., 2015, 55 (6), pp 1192–1201. I used a fullerene model
> without charges, thus ΔG-solv = ΔG-vdw. I attached in this e-mail a plot
> of  ⟨dUvdw(λ)/dλ⟩ vs λ for this example. As you can see, from λ=0 to
> λ=19 are mainly dispersive contributions and from λ=19 to λ=24 are
> repulsive contributions. Most probably 10 λ (and equally sized) are
> enough to obtain converged dispersive contributions (but not 20 since
> I'm overkilled the efficiency) and 15 λ (maybe less) will be required to
> get converged repulsive contributions (but not 5, because it's not
> enough). However, I could not know this information a priori, before to
> run any simulation, (ie, how many λ to use and how to distribute them)
> if I do the decoupling of dispersion and the repulsive contributions
> simultaneously. In general, it will be most efficient if I can decompose
> separately all nonbonding interactions, since I know a priori that
> converged ΔG-elec and ΔG-disp values will be obtained with only few and
> equally sized λ while I will put more and unequally sized λ (ie., more
> effort) to get converged ΔG-rep values.
>
> As you suggested, I tried to make the van der Waals decomposition
> modifying the topology files, however I didn't figure out how to make
> it. Could you provide me more details about this procedure?. I
> understand that to modify the code is not an easy task, it involves a
> lot of effort writing, validating and optimizing the code, however it
> will be more flexible if the three decomposition could be done
> separately from the mdp input file without touching the topology file.
>
> I will post this question about the van der Waals decomposition in the
> gmx-users discussion, thus other users with the same doubt can be
> informed about the possible solutions.
>
> I appreciate your time and comments on this issue.
>
> Joel
>
> 2016-11-07 16:16 GMT+01:00 David van der Spoel <spoel at xray.bmc.uu.se
> <mailto:spoel at xray.bmc.uu.se>>:
>
>     On 07/11/16 15:55, Joel Jose Montalvo Acosta wrote:
>>     Dear Dr. van der Spoel,
>>
>>     My name is Joel Montalvo, I'm a PhD student computational
>>     chemistry in the University of Strasbourg (France). I'm interested
>>     in computing solvation free energy of solutes with apolar
>>     character as fullerenes and graphene in organic solvents, using
>>     free energy perturbation (FEP) and thermodynamic integration (TI)
>>     techniques as implemented in gromacs.
>>
>>     I write because I would like to decouple independently the
>>     dispersion and repulsion components of the Lennard-Jones term in a
>>     FEP/TI run such as it happens with electrostatics and
>>     Lennard-Jones interactions. How can I modify the gromacs source
>>     code in order to make this decouple independenly?. I image that in
>>     the mdp input file, the options couple-lambda0 and couple-lambda1
>>     could accept new values as "disp" for dispersion and "rep" for
>>     repulsion. Also, it is possible to have new
>>     lambda-vector sections, "disp-lambdas" and "rep-lambdas" for
>>     controlling independently both contributions during a FEP/TI run.
>>
>>     I find interesting this decomposition because, first, the
>>     convergence profiles for both contributions are different, so I
>>     could use a different number of lambda for dispersion than for
>>     repulsion, similar to what happens with electrostatics and
>>     Lennard-Jones contributions. Second, it will be valuable to have
>>     the dispersion/repulsion contribution ratio for solvation free
>>     energies of different solutes.
>>
>>     Thank you in advances for any comment or help
>>
>>     Joel Montalvo Acosta
>>
>     This is a regular gromacs question so please continue the discussion
>     on the gromacs user list.
>     I think you can simply define A and B values for C6 and C12 in your
>     topology file and the do the TI in multiple steps. You will have to
>     turn off the C6 first I guess, since the goes to -infinity. However,
>     then I don't see how you could ever get converged values for the
>     next step, since then you are left with a completely repulsive
>     potential (it might work in NVT though).
>
>     However before you even start you should ask yourself what it means.
>     Will you get any meaningful results? Decomposing in enthalpy and
>     entropy is meaningful, but energy components is not.
>
>     --
>     David van der Spoel, Ph.D., Professor of Biology
>     Dept. of Cell & Molec. Biol., Uppsala University.
>     Box 596, 75124 Uppsala, Sweden. Phone:	+46184714205 <tel:%2B46184714205>.
>     spoel at xray.bmc.uu.se <mailto:spoel at xray.bmc.uu.se>    http://folding.bmc.uu.se
>
>


-- 
David van der Spoel, Ph.D., Professor of Biology
Dept. of Cell & Molec. Biol., Uppsala University.
Box 596, 75124 Uppsala, Sweden. Phone:	+46184714205.
spoel at xray.bmc.uu.se    http://folding.bmc.uu.se