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

[David Mobley]

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