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

Erik Lindahl lindahl at sbc.su.se
Sat Apr 1 20:33:35 CEST 2006 

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

More information about the gromacs.org_gmx-developers mailing list