[gmx-developers] more on electrostatic decoupling

[Michael Shirts]

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

Berk or Erik-

Could you address this part of David Mobley's question?  It seems like
it would be necessary to send these interactions to a separate
neighborlist that
did not get truncated in the same way.

I thought I'd summarize a little more of what I see the issues are for
some of the people who have not been in on the conversations offline.
What we would LIKE to do calculate is the free energy of transferring
an isolated assembly of partial charges called a "solute" into an
infinite bath of classical water particles.  What we end up
calculating with PME involves an infinite periodic lattices of both
types of particles.

"Charge annihilation"  corresponds to one end state being the periodic solute
in periodic water, and the other end state being an uncharged periodic
solute + periodic water (no solute-water interactions).  The free
energy of uncharged periodic solute -> uncharged isolated solute is
zero, as it takes no work to separate uncharged particles (neglecting
vdW greater than the box length).  Assuming we can calculate the
charging free energy of the molecule in isolation (not necessarily
simple if it is flexible, might require a vacuum simulation that will
be hard to with MD because of ergodicity issues with small systems,
and if charged there are problems with the changing "charge jelly"
used to yield cell neutrality) we then have one end state being the
periodic solute in periodic water, and the other end state being
isolated solute + periodic water.

For "charge decoupling", the end states are the periodic solute in periodic
water, and the other end state is periodic solute + periodic solvent
(no interactions).   It seems it wouldn't be that difficult to compute
the periodic solute -> isolated solute, but there may be complications
(periodic space is not rotationally symmetrical, if they are not
rigid, there is some polarization, etc). Certainly, if we are
interested in ligand binding, with the same size systems in ligand
binding for the complex and solvent decoupling simulations, this
periodic -> isolated term is going to identically cancel.  Otherwise
(only one simulation leg, different size boxes), things get more
difficult . . .

I think the general consensus in the field is that periodic PME pure
water is the best approximation to bulk water that we can get right
now with small simulations, and with a decent size system (even as
small as 1 nm or so, but certainly at 2-3 nm), any difference is
pretty small.

But one thing that is unclear to me is the free energy to go from periodic
solute in periodic water to isolated solute in infinite water.  My
point (though I probably got my facts wrong) is that we don't really
understand this transformation.  The proposition put forward seems to
be (total dipole per box -> 0), therefore the free energy is zero.  But I don't
know that we have really articulated well what this actually does
consist.  It's even a little clear what process this transformation
corresponds to.  Is it -inserting- other water molecules, increasing
the box size further, until the periodic copies of the solute have no
interaction?  The free energy of that process is nonobvious to me. 
But somebody must have studied this at some point -- box size
dependence of HF free energy solvation, perhaps?  Though you would
need to do both the "decoupling" and "annihilation" approaches here,
of course . . . although the same-cell HF Coulombic interaction is
assumed to be zero in molecular mechanics models, the intercellular HF
interaction changes depending on the treatment.

Just brainstorming here, it seems like it would be possible to write a
version of PME where the "mutated" version of the ligand NEVER (at
either endpoint) has interactions with its self copies in other cells
-- we just need to subtract out this term at all points.  It may
require an extra Ewald calculation, but since it will involve a small
number of atoms, it might not be a huge overhead.

What would this imply . . . that we would have is one single solute
floating in a sea of periodic water with an infinite series of small
cavities.  But the solvent would be slightly polarized around these
cavities, especially with a charged solute. . . . hmm.  This has
complications as well. There might to be an endpoint polarization
problem at the -coupled- end this time.  :)

Best,
Michael