Dear Developers, <br>
<br>
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.<br>
<br>
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.<br>
<br>
Berk had previously suggested this:<br>
>For decoupling you want a decoupled state where the decoupled<br>
>molecule has plain Coulomb interactions without cut-off.<br>
>This is not possible with the normal neighborlists, as there is<br>
>no guarantee that the molecule is smaller than the cut-off length.<br>
>I think the simplest way to accomplish decoupling is with pair<br>
>interactions.<br>
<br>
>My proposition is:<br>
<br>
>Add a pairs interaction type 2 without parameters that uses<br>
>the normal LJ parameters and a FudgeQQ of 1.<br>
>Add pairs type 3 and 4 that are equivalent to 1 and 2,<br>
>except that they always use the A state parameters.<br>
<br>
>Add a decouple option to the mdp file where the user can<br>
>select a moleculetype.<br>
>For this moleculetype grompp should then change all pairs<br>
>type 1 to type 3 and for all non-excluded atom pairs add<br>
>exclusions and pairs type 4.<br>
<br>
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).<br>
<br>
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. <br>
<br>
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)? <br>
<br>
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).<br>
<br>
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: <br>
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).
<br>
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.<br>
<br>
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. <br>
<br>
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?<br>
<br>
Thanks very much for your time,<br>
David Mobley<br>
Dill group<br>
UCSF<br>
<br>
<br>