[gmx-developers] rcoulomb, rlist

[David Mobley]

Berk,

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

It seems surprising that this would do so much here. Perhaps there is
a way to test this -- if I remember correctly beta is also controlled
by ewald-rtol, so I could mimic the effects of altering rcoulomb from
0.9 to 1.0 by changing the relative tolerance?


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

Yes, I was in on those discussions, but needed something like this to
refresh my memory.

> 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?

I don't believe so, for the following reasons:
1) qA and qB are both zero for the atoms in acetamide for this calculation.
2)

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

Yes. In fact in this case we disagree with Shirts' value for acetamide
(and more so with rcoulomb=rlist than with rcoulomb<rlist), which is
one of the things that concerns us. I hand-checked our FFAMBER
tip4p-ew implementation and it is OK, and it does give the correct
density, so we're very perplexed here.

Also, for what it's worth, our TIP3P results for acetamide are in
agreement with yours and Michael's and don't appear to exhibit this
large dependence on rcoulomb, which makes matters even more
perplexing. Is there anything water-model-specific in the PME code?

Thanks,
David

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