[gmx-developers] suggestion for PME for 3.3.1 for free_energy=yes
David Mobley
dmobley at gmail.com
Sat Apr 1 20:39:59 CEST 2006
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.
>
More information about the gromacs.org_gmx-developers
mailing list