[gmx-developers] more on electrostatic decoupling

[David Mobley]

Dear Developers,

John Chodera and I have been working with Michael Shirts to try and figure
out how to do free energy calculations where the ligand's electrostatic
interactions with the protein are turned off in the course of a simulation,
but the ligand-ligand interactions are left on. There was a brief exchange
on this list about this a little while ago in response to an e-mail from
Michael, and I wanted to follow up on it, as I'm not sure we've totally
resolved all of the issues. I'll be excerpting some from the previous
exchange below.

First, my main question is this: Is it easily possible with the current way
things are implemented in GROMACS to do as Berk had previously suggested
(just below) and make the ligand decoupled state still have full
ligand-ligand electrostatic interactions with no cutoff? For reasons I'll
explain below, it appears the alternate solution (recomputing ligand-ligand
interactions with PME in the decoupled state) is undesirable.

Berk had previously suggested this:
>For decoupling you want a decoupled state where the decoupled
>molecule has plain Coulomb interactions without cut-off.
>This is not possible with the normal neighborlists, as there is
>no guarantee that the molecule is smaller than the cut-off length.
>I think the simplest way to accomplish decoupling is with pair
>interactions.

>My proposition is:

>Add a pairs interaction type 2 without parameters that uses
>the normal LJ parameters and a FudgeQQ of 1.
>Add pairs type 3 and 4 that are equivalent to 1 and 2,
>except that they always use the A state parameters.

>Add a decouple option to the mdp file where the user can
>select a moleculetype.
>For this moleculetype grompp should then change all pairs
>type 1 to type 3 and for all non-excluded atom pairs add
>exclusions and pairs type 4.

Michael was unclear that this is possible with the way interactions are
currently handled in GROMACS, but I took from Berk's suggestion that the [
pairs ] section must be handled separately from the neighborlist. Can
someone more familiar with the code answer this definitively? That is, if I
specify in the [ pairs ] section an interaction between two atoms which are
separated by MORE than the electrostatic or vdW cutoff, will that
interaction still be evaluated? If so, it sounds like Berk's solution is the
way to go, as this really seems to be the best option. And if so, I could
use some tips on where to go in the code to modify how this pairs section is
handled (that is, to add these pair types).

I do object to the suggestion about grompp. One thing that's essential for
calculating  *standard* free energies of binding is applying restraints
between the ligand and the protein in the decoupled state (both to ensure
convergence and to get the standard volume into the calculation somehow).
See Boresch and Karplus, Gilson, Hermans, and others for this. Anyway, the
point is, it's necessary to restrain the ligand to the protein to calculate
absolute binding free energies. And currently, GROMACS requires that
restraints be *within* the same "molecule"; that is, the ligand and protein
must be defined in the same molecule section in order to restrain the ligand
to the protein. But Berk's suggestion about grompp would make it so the
ligand and the protein would have to be in *different* molecules in order
for grompp to automatically fill in the parameters. If there isn't a
workaround, this would make it impossible to restrain the ligand to the
protein.

Presumably there is some other way to have grompp fill in the parameters?
Perhaps make it so grompp can be provided with an .ndx file containing the
atom list for the atoms to be decoupled (the ligand)?

Aside from the grompp issue, though, I think this makes the most sense in
terms of how to deal with this, so I just want to make sure it's feasible to
implement in the current code (and am willing to work on it in the near
future myself).

Michael's suggestion had been, essentially, to evaluate PME separately for
the ligand in the decoupled state. But we've discussed this with him off
list, and I think that Berk was correct in pointing out that this means that
the ligand will have long range interactions with copies of itself in
vacuum. This appears very problematic, actually: Here are two specific
problems we've thought of based on Berk's line of reasoning:
1) If the ligand has a dipole moment, there will be long range dipole-dipole
interactions between the ligand and copies of itself which will result in
both forces and probably significant contributions to the energy. (The most
extreme case, as pointed out by John Chodera, is a molecule like HF in a
cubic box aligned with one of the box axes. It will be interacting with
copies of itself which are aligned in the same direction, producing a pretty
big net electric field and force keeping dipoles aligned in that same
direction. There should be fairly big energies too).
2) Problems will be box-size dependent. In binding free energy calculations,
you might hope that any mistake you make on the ligand-in-protein part of
the calculation would cancel with the mistake for the ligand-in-water
component. But this won't be the case since the distance between the ligand
and copies of itself in this PME vacuum component are set by the box size,
which will in general be different for the two systems.

The magnitude of the error that these two things would introduce seems
unclear and system-dependent. It may not be large for some systems (for
example, for molecules without a dipole moment, like methane, I have some
results which indicate it's not problematic), but it doesn't seem like it
can be a very good thing. On the other hand, since we're free to choose any
path between the bound state and the unbound state we want as long as we
close it properly, Berk's solution seems perfectly valid.

The only other concern that came up in our e-mail discussion on this is that
Berk's solution might not work very well for any calculation involving
changing more than one molecule at once. Any thoughts on this aspect?

Thanks very much for your time,
David Mobley
Dill group
UCSF
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20060116/3d355f84/attachment.html>