[gmx-developers] rcoulomb, rlist

[Berk Hess]

Hi,

As far as I know, rcoulomb < rlist with PME will only change the Ewald beta
and nothing else.

We have had long mail discussions (including you in the cc?)
in which I explain why we can not have rcoulomb < rlist with PME.
The reason is that Gromacs uses spline interpolated tables without if 
statements,
therefore we can not have abrupt cut-off's (which would give 
(near-)infinite forces).
I don't know how other packages solve this. I guess just by using a 
if-statement,
which would mean you have a jump in the energy and the energy is no longer
the integral of the force.

But your free energy difference seem to indicate a larger error than a small
change in beta can explain.
There was a bug with PME and free energy in pre 3.3.3 versions for atoms 
with qA=0, qB!=0,
which I fixed for 3.3.3.
Could this be what you are seeing?

That is the only bug in pre 3.3.3 version I know of.
I got perfect agreement with Michael Shirst results and also with RF results
and GROMOS package results with RF.

Berk


David Mobley wrote:
> All,
>
> I need some help understanding the behavior of GROMACS 3.3. (
> Apologies on the reference to the old version, but I published some
> data that now others are having trouble reproducing, which
> necessitates me using the version I originally ran the calculations
> in).
>
> I am computing the hydration free energy of acetamide in TIP4P-Ew
> water using the FFAMBER ports, and I am having a problem in the
> component of the calculation where the A state has acetamide in water
> with no partial charges on it, and the B state has acetamide in water
> with all of its interactions with water turned off (but intramolecular
> interactions on) which is accomplished through some topology file
> trickery (modifications to the pairs section of the topology and
> explicit specification of intramolecular nonbonded interactions).
>
> Anyway, the expected value for this component calculation is something
> around -1.8 kcal/mol. My published result (in GROMACS 3.3) was -2.1
> kcal/mol (+/- 0.05 or so). I tried to reproduce the result of -2.1
> kcal/mol to confirm that the difference was legitimate, and found that
> I instead got -4.1 kcal/mol roughly.
>
> After some testing, I established that the difference was due to a
> single setting: In my original calculations, I had used PME, with
> rcoulomb = 0.9 nm and rlist = 1.0 nm. When I tried to reproduce these,
> I used PME, with rcoulomb=1.0 nm and rlist=1.0 nm.
>
> Can someone help me understand what exactly is going on here with this
> (seemingly small) change in the mdp files making such a large
> difference in the computed hydration free energies? I realize that in
> later (3.3.1 and after) versions of gromacs, the code insists that
> rcoulomb=rlist with pme, essentially (which still seems strange to
> me), but I still am not clear on why this is (and what exactly is
> going wrong in the calculations).
>
> For what it's worth, I still get the correct density of water with
> both sets of parameters. Yet this parameter change seems to cause a
> *systematic* difference of some sort -- for example, if I look at the
> plot of dV/dlambda values
> (http://www.dillgroup.ucsf.edu/~dmobley/uploads/novdw_comparison.pdf)
> the differences between the two cases are entirely within a specific
> range of lambda values, and are in a consistent direction. Yet I find
> no systematic differences in system volume or density at these lambda
> values.
>
> That's my question -- any help will be appreciated.
>
> The astute will notice I now also have a separate problem: If the
> "correct" input parameters to use are rcoulomb=rlist=1.0 nm, then the
> correct result is -4.1 kcal/mol, which is even further from the
> expected answer (-1.8 kcal/mol) than I started off with. This is the
> expected answer, in that it is the answer I arrive at using other
> simulation packages (DESMOND, for example) with the same input
> parameters, and the answer my collaborators (using several other
> packages) are getting with the same input parameters. Obviously this
> is a cause for some concern and suggests there may be some sort of an
> underlying problem here. On my list of things to do is to see whether
> I still get the same results in GROMACS 4.
>
> Thanks for your help,
> David Mobley
> _______________________________________________
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the 
> www interface or send it to gmx-developers-request at gromacs.org.
>
>