[gmx-developers] suggestion for PME for 3.3.1 for free_energy=yes

Erik Lindahl lindahl at sbc.su.se
Sat Apr 1 21:37:09 CEST 2006 

Hi,

OK, if you have nodes available it would be great if you could test  
the pme_order = 4 vs 6 cases in double, or if you just send me a  
tarball with your tpr files I can do it.

I'll still read through pme.c in the meantime :-)

Cheers,

Erik


On Apr 1, 2006, at 8:39 PM, David Mobley wrote:

> Erik,
>
> This was in single precision.
>
> Incidentally, Michael Shirts pointed out that I hadn't tried *exactly*
> the default, which has ewald_rtol=1e-05, so I have that running now. I
> can update later today or Monday with that information.
>
> Let me know what you figure out.
>
> Thanks,
> David
>
>
> On 4/1/06, Erik Lindahl <lindahl at sbc.su.se> wrote:
>> Hi David,
>>
>> Great work. Was this too with double precision?  In that case I can
>> rule out any accuracy-related bug in pme.c, and will add a warning.
>> If it was single I'll spend an hour hunting for accuracy-sensitive
>> operations in pme.c first ;-)
>>
>> Cheers,
>>
>> Erik
>>
>> PS: Only two bugs to before 3.3.1...
>>
>>
>> On Apr 1, 2006, at 12:53 AM, David Mobley wrote:
>>
>>> All,
>>>
>>> Just reporting back on my tests of various combinations of PME
>>> parameters and rcoulomb. Here's a little table of what I've done.
>>> These are computed charging free energies of phenol with AMSOL  
>>> charges
>>> (which I know is a little strange) in TIP3P water using five lambda
>>> values, Bennett Acceptance Ratio, and GAFF parameters for the  
>>> ligand.
>>> Using this protocol the error in these calculations should be quite
>>> low (0.02-0.1 kcal/mol, in line with my computer uncertainties). And
>>> sorry, my units are in kcal/mol. Also, these are in a dodecahedral
>>> box.
>>>
>>> rcoul  PME order       ewald rtol fourier spacing      delta G    
>>> node hrs/ns
>>> 0.9    4       1e-08   0.12    19.51+/-0.04
>>> 0.9   6        1e-06   0.10    18.77+/- 0.04   2.8
>>> 1.0   6        1e-06   0.10    18.83+/- 0.03   3.8
>>> 0.9    8       1e-06   0.08    18.75+/-0.03    5.7
>>> 1.0    4       1e-08   0.12    18.86+/-0.03    3.36
>>> 1.1    4       1e-08   0.12    18.80+/- 0.02   4.06
>>>
>>> OK, so anyway, I think the moral of all this is that everything  
>>> works
>>> fine *except* rcoul=0.9 (or, presumably, anything less than 0.9)  
>>> with
>>> the GROMACS default PME parameters. Maybe issuing a warning if these
>>> parameters are used?
>>>
>>> Also, I would point out that, at least for the system, it *IS*  
>>> faster
>>> to use rcoul=0.9 with PME order 6 and fourier spacing 0.10 than  
>>> to go
>>> to rcoul=1.0 with the original PME parameters. This also seems to be
>>> the case in a protein-ligand system (26,000 atoms) I've looked at.
>>>
>>> I would point out that while the difference from the first line  
>>> above
>>> to the other lines in the computed free energy is relatively small,
>>> this is a substantial fraction of the hydration free energy for a
>>> typical small molecule, and thus could mess free energy calculations
>>> up a whole lot (in fact, we already had to rerun a bunch when we
>>> realized this). So I suggest at least issuing a warning (or  
>>> putting a
>>> warning in the free energy section of the documentation?).
>>>
>>> It is possible, of course, that this affects other observables as
>>> well. I suppose I can compute water densities and things too if  
>>> people
>>> aren't convinced this is a problem.
>>>
>>> Thanks,
>>> David
>>>
>>>
>>>
>>> On 3/29/06, David Mobley <dmobley at gmail.com> wrote:
>>>> Erik,
>>>>
>>>> These are in dodecahedral boxes. So far I prefer them because they
>>>> save some water molecules for things that are fairly spherical,
>>>> although if there is some technical reason these might be worse,  
>>>> I am
>>>> happy to switch.
>>>>
>>>> Thanks,
>>>> David
>>>>
>>>> On 3/29/06, Erik Lindahl <lindahl at sbc.su.se> wrote:
>>>>> Hi David,
>>>>>
>>>>> Great. Right now 3.3.1 is our priority, but then I'll write some
>>>>> code
>>>>> to actually check the dG/dl values in the PME routine.
>>>>>
>>>>> One more thing - are you using rectangular or triclinic (truncated
>>>>> octahedron/rhombic dodecahedron) boxes?
>>>>>
>>>>> Cheers,
>>>>>
>>>>> Erik
>>>>>
>>>>> On Mar 30, 2006, at 12:07 AM, David Mobley wrote:
>>>>>
>>>>>> Erik,
>>>>>>
>>>>>>> I assume that your values converge if you use even tighter  
>>>>>>> limits
>>>>>>> (e.g. order 8 or spacing 0.8)?
>>>>>>
>>>>>> It looked like it, although I must admit I didn't go to order  
>>>>>> 8, so
>>>>>> I'll start that running now (should be done tomorrow).
>>>>>>
>>>>>>> However, instead of increasing the amount of PME work I would
>>>>>>> recommend checking what happens with longer direct-space
>>>>>>> cutoffs. The
>>>>>>> default PME values are pretty similar to Tom Darden's default
>>>>>>> ones,
>>>>>>> and since Gromacs is much faster on direct-space interactions  
>>>>>>> than
>>>>>>> his code we should move more interactions there.
>>>>>>
>>>>>> Benchmarking here and by Michael suggested it was cheaper to
>>>>>> increase
>>>>>> the pme_order and fourier spacing than increase rcoulomb,
>>>>>> although I
>>>>>> think Michael's suggestion would be to do both  
>>>>>> (rcoulomb>=1.0). I'm
>>>>>> going to start a couple runs going to check the case of rcoulomb
>>>>>> larger (1.0 or 1.1) but PME order 4 still with the 0.12 A
>>>>>> spacing, as
>>>>>> well as rcoul=0.9 (my original) with order 8 and spacing 0.08.
>>>>>> Let me
>>>>>> know if you'd like me to check other combinations, although  
>>>>>> not too
>>>>>> many (I run these for about a day on 5 processors each).
>>>>>>
>>>>>> Thanks,
>>>>>> David
>>>>>>
>>>>>>>
>>>>>>> Cheers,
>>>>>>>
>>>>>>> Erik
>>>>>>>
>>>>>>> On Mar 29, 2006, at 11:33 PM, David Mobley wrote:
>>>>>>>
>>>>>>>> Dear all,
>>>>>>>>
>>>>>>>> I wanted to write and suggest that the default PME parameters
>>>>>>>> be set
>>>>>>>> differently when free_energy=yes in the 3.3.1 release. In
>>>>>>>> particular,
>>>>>>>> I think the default pme_order needs to be 6 with a fourier
>>>>>>>> spacing of
>>>>>>>> 0.10 or better. ewald_rtol probably isn't as important but I
>>>>>>>> would
>>>>>>>> suggest 1e-06.
>>>>>>>>
>>>>>>>> Currently, in 3.3, the default pme_order=4 with
>>>>>>>> ewald_rtol=1e-05 and
>>>>>>>> fourierspacing=0.12. Michael Shirts and I have done some
>>>>>>>> testing (he
>>>>>>>> can comment on his testing separately if he likes) and found  
>>>>>>>> that
>>>>>>>> with
>>>>>>>> pme_order=4, results for free energy calculations involving
>>>>>>>> modifications of electrostatics can be substantially wrong.
>>>>>>>>
>>>>>>>> My particular testing has mostly been done with rcoulomb=1.0
>>>>>>>> or 0.9,
>>>>>>>> although Michael has looked more at varying this, as well.
>>>>>>>> Regardless,
>>>>>>>> I find that results for computed free energies involving
>>>>>>>> changes in
>>>>>>>> electrostatics can be wrong by 1-4 kcal/mol with the default  
>>>>>>>> pme
>>>>>>>> parameters as compared to their values with higher pme_order  
>>>>>>>> and
>>>>>>>> lower
>>>>>>>> fourierspacing. For my own work I think I've ended up  
>>>>>>>> settling on
>>>>>>>> pme_order=6 and fourier_spacing=0.10 with ewald_rtol=1e-6,  
>>>>>>>> as my
>>>>>>>> results don't seem to change substantially if I further reduce
>>>>>>>> fourier_spacing or ewald_rtol, or even increase rcoulomb.  
>>>>>>>> Again,
>>>>>>>> though, I want to emphasize that the difference in free  
>>>>>>>> energies
>>>>>>>> going
>>>>>>>> from the default PME parameters to the ones I suggest here  
>>>>>>>> is 1-4
>>>>>>>> kcal/mol (I tested this on some polar molecules including  
>>>>>>>> phenol;
>>>>>>>> the
>>>>>>>> difference is smaller for less polar molecules), which is  
>>>>>>>> roughly
>>>>>>>> comparable to the total hydration free energy for such  
>>>>>>>> molecules.
>>>>>>>>
>>>>>>>> If for some reason you decide it's not a good idea to change  
>>>>>>>> the
>>>>>>>> default PME parameters (at least for free_energy=yes), I would
>>>>>>>> strongly encourage adding some suitable warnings if
>>>>>>>> free_energy=yes
>>>>>>>> with these parameters, and probably also in the documentation.
>>>>>>>>
>>>>>>>> I can provide exact values for the computed free energies,
>>>>>>>> etc. as
>>>>>>>> needed. It is possible of course that the *best* parameters for
>>>>>>>> these
>>>>>>>> things are system-dependent, but it seems pretty clear that
>>>>>>>> pme_order=4 with fourier_spacing=0.12 is NOT adequate for free
>>>>>>>> energy
>>>>>>>> calculations involving modification of charges, at least  
>>>>>>>> when the
>>>>>>>> charges are reasonably polar.
>>>>>>>>
>>>>>>>> Thanks,
>>>>>>>> David Mobley
>>>>>>>> UCSF
>>>>>>>> _______________________________________________
>>>>>>>> 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.
>>>
>>
>> _______________________________________________
>> 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