[gmx-developers] Gromacs 3.3.4 PME energies keep changing with longer cutoffs

[Mark Abraham]

On 17/02/10 06:25, Gabriel Rocklin wrote:
> Hi all,
>
> I'm seeing some confusing behavior with Gromacs 3.3.4 and PME
> electrostatics.  My potential energies seem to not to converge a single
> "correct" value as I make the PME parameters more and more exact.

You're actually not converging PME meaningfully. ewald_rtol and rcoulomb 
determine the value of a parameter (beta of manual 4.9.1) that controls 
the relative importance of the direct- and reciprocal-space terms. As 
manual 7.3.13 notes, decreasing ewald_rtol has contrary effects on the 
accuracy of the two components. GROMACS does not have a way of making 
both halves more accurate at constant beta. Converging accuracy requires 
having sufficiently large rcoulomb and small fourierspacing for the 
prevailing ewald_rtol. You can vary all three sensibly, and need to do 
so in order to maximize speed for given accuracy.

> I took 10 snapshots from a single precision, constant volume simulation
> I'm doing of a protein in water and saved them as .gro files with -ndec
> = 8.  I then used mdrun -rerun to re-evaluate their energies using a
> number of different PME parameters.  I found that the energies were not
> sensitive to ewald_rtol or fourierspacing (below certain values), but
> they were very sensitive to rlist and rcoulomb even out to very long
> distances.
>
> For example, using pme_order=6, spacing=0.1, and tolerance=0.00000001,

Your choice of ewald_rtol=1e-8 is a lot smaller than typical values, so 
already your reciprocal-space contribution is much less accurate. 
Increasing pme_order will help, but you should also observe the effect 
of decreasing the spacing. Since the spacing is not used directly, but 
is used to determine the dimensions of a grid subject to various 
constraints, I suggest you try explicit grid sizes fourier_n[xyz] that 
are multiples of your number of PME-only nodes and suitable multiples of 
the Fourier factors (see 4.9.4).

> as I increase rlist and rcoulomb together from 0.7 to 2.5nm, I get these
> potential energies for my 10 snapshots:
>
> 0.7 -656927.375 -657225.4375 -656475.0 -656451.875 -657604.8125
> -656090.6875 -657553.25 -657241.8125 -657747.9375 -656863.5625
> 0.9 -654605.375 -654916.1875 -654163.9375 -654132.125 -655300.875
> -653771.5 -655239.875 -654936.875 -655420.75 -654539.25
> 1.0 -654553.5625 -654862.0625 -654108.9375 -654078.0 -655246.875
> -653718.4375 -655182.6875 -654884.3125 -655365.8125 -654486.375
> 1.1 -654536.875 -654842.5 -654094.3125 -654061.25 -655227.9375
> -653701.4375 -655166.5 -654864.3125 -655349.625 -654468.75
> 1.2 -654533.5625 -654842.625 -654087.1875 -654054.5 -655224.625
> -653695.125 -655164.25 -654863.5625 -655340.5 -654463.125
> 1.3 -654531.9375 -654837.0625 -654087.5625 -654054.5625 -655223.6875
> -653694.125 -655161.6875 -654857.5 -655345.0625 -654462.5
> 1.5 -654532.375 -654839.8125 -654086.1875 -654057.25 -655223.25
> -653693.3125 -655164.4375 -654864.1875 -655340.1875 -654461.5625
> 1.8 -654527.125 -654840.3125 -654086.9375 -654054.0625 -655221.25
> -653696.6875 -655155.75 -654860.0625 -655343.5 -654465.8125
> 2.0 -654532.5625 -654843.125 -654083.75 -654053.6875 -655223.375
> -653695.375 -655158.8125 -654860.1875 -655342.3125 -654464.4375
> 2.5 -654535.875 -654838.75 -654089.875 -654055.0625 -655220.25
> -653692.75 -655163.4375 -654857.375 -655337.125 -654461.375
>
> So not only do the numbers keep changing, but the relative energies
> change as well, so that for example going from 2.0 to 2.5 nm, the
> standard deviation of the change in energy across the 10 snapshots is 4 kJ.
>
> The system is a dodecahedral box with a volume of about 500 nm^3.
>
> Is this normal behavior?  Should these values converge to specific
> energies?  And how far out would you normally have to set rcoulomb/rlist
> to see this?

You ought to see the real-space component of PE converging to machine 
precision as you vary rcoulomb up to about 2 nm for ewald_rtol around 
1e-5 and nk[xyz] a Fourier multiple of npme that gives a spacing around 
0.1nm.

Mark