[gmx-developers] rcoulomb, rlist

[hessb at mpip-mainz.mpg.de]

Him

There is no water specific free-energy code.
I also don't recall any other bugs than the one I mailed.

The one I mailed would actually be consistent with your observations.
The bug causes atoms with qA=0,qB!=0 to be missing
in the PME mesh part.
If you use 0.9 iso 1.0 nm, more energy goes to the mesh part.
If you are affected by the bug the shorter cut-off would
give a larger discrepancy.

Did you modify any code?
Or did you only do fancy topology stuff.

Berk


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