[gmx-developers] rcoulomb, rlist

[David Mobley]

Berk,

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

OK. I should probably also test whether results in 4.0 agree with
those from 3.3 (with PME bugfix) and 3.3.1.

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

Doublechecked my topology. qA=0, qB=0. I'm attaching a tarball with my
top and gro files.

>
> Did you modify any code?

No.

> Or did you only do fancy topology stuff.

Yes, in that I modified the pairs list and nonbond_params sections to
maintain intramolecular vdW interactions.

Please see attached.

Thanks,
David



> 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.
>>
>
>
> _______________________________________________
> 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.
>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: debug.tar.gz
Type: application/x-gzip
Size: 37647 bytes
Desc: not available
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20090115/4d794713/attachment.gz>