[gmx-developers] rcoulomb, rlist

hessb at mpip-mainz.mpg.de hessb at mpip-mainz.mpg.de
Sun Feb 15 22:25:02 CET 2009 

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

You mean qA and qB for the solute 0,
but I assume you still have charge on the water, right?
So the same bug would still affect your results
when you change only the LJ on the solute between A and B,
but you can still have water charge interactions which could
be missing.

Berk

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

More information about the gromacs.org_gmx-developers mailing list