[gmx-developers] Soliciting feedback on one aspect of free energy calculation in 4.6

[Shirts, Michael (mrs5pt)]

I'd like to solicit a bit of feedback on one free energy feature in 4.6.

When preparing for 4.5, Berk noticed that with some molecules with large
charges and small sigmas on the H's (below 0.3), then when using a softcore
pathway to disappear the entire molecule, there were minima created, where
the interactions with the H were significantly stronger than the
interactions with heteroatom H was bonded to.  These "artificial"
intermediate minima can slow sampling and increase statistical error, and
should be avoided when possible.

The solution put in place for 4.5 at that point was that whenever the sigma
was lower than sc_sigma (which defaults to 0.3, and was previously only
imposed for molecules with zero sigma on C6), it was reset to sc_sigma.
This does not change the endpoints (so the free energies were always still
correct): however, it did change the pathway, which created a difference
from previous versions, making comparisons more difficult.

This change was documented, but it was not necessarily clear what the
implications were. This is equivalent to, when sigma drops below sigma_min,
the effective softcore alpha parameter being increased.  Specifically:

The formula (older code):
sigma^6/(sigma^6*alpha*lambda + r^6)
becomes (4.5):
sigma^6/(sigma_min^6*alpha*lambda + r^6

The first is equal to:
1/(alpha*lambda + (r/sigma)^6
The second is:
1/((sigma_min/sigma)^6*alpha*lambda + (r/sigma)^6

Assuming sigma = 0.25 for the H-atoms in question (usually about right),
then alpha essentially becomes (0.3/0.25)^6*0.5 = 1.5 for these H atoms, and
stays 0.5 for the other atoms, with rapidly increasing alpha for smaller
sigma.  Looking through some of the literature, it seems to indicate that
these higher alpha than 0.5 could decrease the inefficiency of the pathway,
though this could use some more study.  It's also not a standard way of
doing things in the literature or in other codes.

My proposal was to keep this behavior for transformations when soft coring
coul+vdw together, but return to the old, standard, behavior for soft core
when coulomb and vdw are turned off separately. I believe this preserves the
most efficient pathway for sc both with and without charges, as well as
returns the vdw-alone soft core the way to it was for many iterations of
gromacs, as well as being what people understand that soft core does from
the literature.  In the free energy changes scheduled for 4.6, it will be
much easier to turn off first vdw and then coul in the same set of
simulations (one .mpd and .top required), and research indicates it's
generally a bit more efficient than changing both vdw and charge
simultaneously.

Any comments would be great!  Please respond soon (next day or two), as
validation of new free energy code needs to be finished soon to not delay
4.6.

Best, 
~~~~~~~~~~~~
Michael Shirts
Assistant Professor
Department of Chemical Engineering
University of Virginia
michael.shirts at virginia.edu
(434)-243-1821