Project

General

Profile

Bug #196

Setting rlist > rcoulomb with coulombtype = pme

Added by John D. over 11 years ago. Updated over 11 years ago.

Status:
Closed
Priority:
Normal
Assignee:
Erik Lindahl
Category:
mdrun
Target version:
Affected version - extra info:
Affected version:
Difficulty:
uncategorized
Close

Description

Peter Kasson noticed that the CVS branch version of grompp complains when we try to set rlist > rcoulomb for coulombtype = pme. 'pme-switch' and 'pme-user' both allow this, but 'pme' doesn't seem to:

ERROR: With coulombtype = PME, rcoulomb must be equal to rlist

In order to ensure forces and gradients (as defined by the forcefield) are correctly evaluated at each timestep, we must be able to ensure particles do not drift into the rcoulomb direct-space interaction radius between pairlist updates. Ideally, this is done with rlist > rcoulomb and automatic pairlist updates, but in the absence of automatic pairlist updates, we can fix the update interfal and select rlist > rcoulomb appropriately.

It seems like it should be a trivial change to enable this, provided the nonbonded kernels don't choke on this.

Thanks!

- John

History

#1 Updated by Berk Hess over 11 years ago

I did this on purpose.
I would say that if you want exact integration
you would use pme-switch.
Having this check in for plain pme avoids that
users choose stupid settings.

Berk.

(In reply to comment #0)

Peter Kasson noticed that the CVS branch version of grompp complains when we
try to set rlist > rcoulomb for coulombtype = pme. 'pme-switch' and 'pme-user'
both allow this, but 'pme' doesn't seem to:

ERROR: With coulombtype = PME, rcoulomb must be equal to rlist

In order to ensure forces and gradients (as defined by the forcefield) are
correctly evaluated at each timestep, we must be able to ensure particles do
not drift into the rcoulomb direct-space interaction radius between pairlist
updates. Ideally, this is done with rlist > rcoulomb and automatic pairlist
updates, but in the absence of automatic pairlist updates, we can fix the
update interfal and select rlist > rcoulomb appropriately.

It seems like it should be a trivial change to enable this, provided the
nonbonded kernels don't choke on this.

Thanks!

- John

#2 Updated by John D. over 11 years ago

Berk,

I am very curious as to why you think setting rlist > rcoulomb is a "stupid setting". Could you elaborate?

Only by setting rlist > rcoulomb is it possible to ensure that the forces are equal to the negative gradient of the potential for more than one timestep after each pairlist build if nstlist > 1. How is this "stupid" exactly?

Furthermore, all other molecular dynamics codes I am aware of use a neighborlist build cutoff to be 1-2 A larger than the direct-space cutoff, for precisely this reason -- you need to build a neighborlist that incorporates a larger region than your direct-space cutoff if you expect the forces to be correctly evaluated for more than one step between each neighborlist build.

In fact, I would suggest you REQUIRE rlist > rcoulomb for colombtype = PME if nstlist > 1. When combined with the new automatic neighborlist update triggering, this would prevent the undesirable situation where the gradient is not correctly evaluated in between neighborlist builds.

Or perhaps this is what is going on already, and I'm simply misunderstanding what these options mean?

Cheers,

John

#3 Updated by Berk Hess over 11 years ago

I think this confusion is caused by the different coding
in the forces of Gromacs compared to other packages.

In Gromacs we have SSE and other assembly loops in which
there are no if statements checking for cut-offs.
Therefore potentials are never zero beyond the cut-off,
unless you supply a potential that is really zero at
the cut-off.
Because we use cubic spline interpolation in Gromacs,
potentials can not be truncated exactly at the cut-off.
Also you will never have exact integration when there
is a cut-off. Cutting of the PME real space part
is still a cut-off, even though the effect is quantitatively
much smaller than cutting of a pure Coulomb interaction.
From your reply I infer this is what other packages do.

So for the way you want to use PME with a ns buffer
(which is surely not a stupid setup), you can use pme-switch
with a very short switching region, since the PME direct
space potential is neary zero at the cut-off, but not exactly.

Berk.

#4 Updated by John D. over 11 years ago

Hi Berk,

Thanks for the clarification. We will switch to pme-switch in the lab.

Some points of response:

1. I think your comment about "exact integration" suggests that I was not clear in my argument. I did not mean "energy conserving dynamics" -- I meant "integrating a known, well-defined function that defines the forces". If nstlist > 1 with pme-switch, in between pairlist builds, some molecules outside the cutoff will have their direct coulomb interactions evaluated and some inside the cutoff won't. This is contrary to how the direct-space interaction is defined. The choice of which evaluations will be incorrect depends on which molecules drifted in and out during integration, and will, for all purposes, be random.

2. Almost nobody is going to use pme-switch unless you make it the default.

3. I very much worry that your implementation of PME differs from every other implementation of PME in every other code as a result. I would highly suggest you find some way to fix this in your SSE loops or whatnot -- I am quite sure it is not impossible, and it may not carry a huge performance hit. Even if it does, you could just use pme-switch instead. But it's probably important that what you mean by PME is the same thing that everyone else means by PME. Or maybe you could rename the current method 'pme-broken' and create a new coulombtype = pme method that implements a zeroing of the forces outside of the cutoff so that you can use rlist > rcoulomb.

4. There should probably be a large warming in the manual that the current method differs from standard practice in all other codes.

Cheers,

John

#5 Updated by Erik Lindahl over 11 years ago

1. I think your comment about "exact integration" suggests that I was not clear
in my argument. I did not mean "energy conserving dynamics" -- I meant
"integrating a known, well-defined function that defines the forces". If
nstlist > 1 with pme-switch, in between pairlist builds, some molecules outside
the cutoff will have their direct coulomb interactions evaluated and some
inside the cutoff won't. This is contrary to how the direct-space interaction
is defined. The choice of which evaluations will be incorrect depends on which
molecules drifted in and out during integration, and will, for all purposes, be
random.

Yes, that's a random error, in contrast to switching of PME which introduces a
systematic error, although the latter does conserve energy.

2. Almost nobody is going to use pme-switch unless you make it the default.

PME-switch is NOT the correct direct space interaction prescribed by Tom Darden and coworkers, and in fact introduces an error in the potential. Now, this might be an "error" you like since it helps you conserve energy better, but it nevertheless increases the average error in the forces.

3. I very much worry that your implementation of PME differs from every other
implementation of PME in every other code as a result. I would highly suggest
you find some way to fix this in your SSE loops or whatnot -- I am quite sure
it is not impossible, and it may not carry a huge performance hit. Even if it
does, you could just use pme-switch instead. But it's probably important that
what you mean by PME is the same thing that everyone else means by PME. Or
maybe you could rename the current method 'pme-broken' and create a new
coulombtype = pme method that implements a zeroing of the forces outside of the
cutoff so that you can use rlist > rcoulomb.

The meaning of PME is defined by the equations in the Darden papers, period.

Frankly, calling an implementation is "broken" when it can conserve energy better in single precision than either Desmond, NAMD, Amber or Charmm do in double is plain stupid.

However, you're more than welcome to contribute code that does the above, and we'll include it as an option.

#6 Updated by John D. over 11 years ago

Hi Erik,

Yes, that's a random error, in contrast to switching of PME which introduces a
systematic error, although the latter does conserve energy.

This is my point. It's a random error that can't be overcome unless we set nstlist = 1. With nstlist > 1, it differs from the implementation described by T. Darden and coworkers.

PME-switch is NOT the correct direct space interaction prescribed by Tom Darden
and coworkers, and in fact introduces an error in the potential. Now, this
might be an "error" you like since it helps you conserve energy better, but it
nevertheless increases the average error in the forces.

I completely agree that it's not the same as the PME described by Darden, but it was Berk who advocated we use pme-switch instead of we wanted gromacs to actually respect rlist > rcoulomb so that nstlist could be set > 1!

" I would say that if you want exact integration you would use pme-switch." - Berk Hess

"So for the way you want to use PME with a ns buffer (which is surely not a stupid setup), you can use pme-switch with a very short switching region, since the PME direct space potential is neary zero at the cut-off, but not exactly." - Berk Hess

Don't yell at me, yell at Berk! :)

The meaning of PME is defined by the equations in the Darden papers, period.

Agreed.

Frankly, calling an implementation is "broken" when it can conserve energy better in single precision
than either Desmond, NAMD, Amber or Charmm do in double is plain stupid.

If you have read the thread here, you'll see that nowhere in my complaint do I even bring up energy conservation as an issue. I'm only concerned with whether the direct-space interaction is "correctly" evaluated, as defined by the Darden papers, each timestep. Because we cannot set rlist > rcoulomb with coulombtype = pme, molecules will drift in and out of the direct-space cutoff radius in between pairlist updates. This means that the direct-space interaction isn't being properly evaluated for all of the atoms within the direct-space cutoff, contrary to what is specified in the Darden papers. That's all I'm worried about here, Erik.

I was hoping that either (1) someone could explain my misconception to me or (2) this would be easy to rectify to bring the implementation in line with the Darden papers. What I'm hearing so far is that "the only way this implementation corresponds to the Darden papers is if you set coulombtype = pme, rlist = rcoulomb, and nstlist = 1". Is that accurate?

However, you're more than welcome to contribute code that does the above, and
we'll include it as an option.

You and I both know there's no way on earth I would be able to figure out how to do anything of the sort on a reasonable timescale, both due to a lack of assembly knowledge on my end and a lack of useful documentation (FLOP counts don't qualify) on the gromacs end.

- John

#7 Updated by Erik Lindahl over 11 years ago

Hi,

Sorry if I sounded a bit grumpy yesterday; long day and I was a bit tired - nothing personal ;-)

My main reason for closing the bug is that it's not "unintended behavior", although it might be useful to have the archived discussion.

(In reply to comment #6)

This is my point. It's a random error that can't be overcome unless we set
nstlist = 1. With nstlist > 1, it differs from the implementation described by
T. Darden and coworkers.

I think that depends a bit on your interpretation of Tom's equations.
First, as I'm sure you're aware, the effect of the switching term is to make the direct space interaction negligible (or rather, decay to less than e.g. 1e-5 relative magnitude) beyond a certain distance.

In practice this means it is possible to use a cutoff. The quality of the solution is judged by the average error magnitude in the forces, and in the reference implementation Tom also uses an "exact PME" alternative where he still includes neighbors significantly beyond the default point where the potential has decayed to 1e-5 to improve the accuracy until it reaches numerical precision. The actual "cutoff" never enters the equation (apart from weighing direct/reciprocal space, of course).

Now, there are other reasons for some of the choices in Gromacs, but when it comes to PME the best solution is clearly the one that produces the smallest error in the forces.

There are a couple of separate issues here. Some can be fixed with minor changes to the tables, but for others there are good technical reasons behind our choices:

1. It's of course technically possible to modify neighborlists such that all atoms that might be within the cutoff are in the list - that's only neighborsearch mods. However, the net result would be that we end up spending a whole lot of time calculating zeros. That might still be smart if it provides is with a more accurate solution, but it doesn't! If we anyway calculate the distances we will reduce the average force error if we evaluate the correct functional form even beyond the cutoff (compare Darden's "exact PME") rather than artificially setting it to zero.

Now, I fully respect that you prefer a different solution with a "well-defined" function to integrate, but compared to the "exact PME" that error is at least of the same magnitude, if not larger.

Of course we could consider including a "sharply truncated PME" table that still sets the potential to zero outside the cutoff (but does not use a switch), but the current implementation is not broken - it makes good use of every cycle available to improve the accuracy!

I completely agree that it's not the same as the PME described by Darden, but
it was Berk who advocated we use pme-switch instead of we wanted gromacs to
actually respect rlist > rcoulomb so that nstlist could be set > 1!

" I would say that if you want exact integration you would use pme-switch." -
Berk Hess

"So for the way you want to use PME with a ns buffer (which is surely not a
stupid setup), you can use pme-switch with a very short switching region, since
the PME direct space potential is neary zero at the cut-off, but not exactly."
- Berk Hess

Don't yell at me, yell at Berk! :)

Berk, consider yourself yelled at ;-)

If you have read the thread here, you'll see that nowhere in my complaint do I
even bring up energy conservation as an issue. I'm only concerned with whether
the direct-space interaction is "correctly" evaluated, as defined by the Darden
papers, each timestep. Because we cannot set rlist > rcoulomb with coulombtype = pme, molecules will drift in and out of the direct-space cutoff radius in
between pairlist updates. This means that the direct-space interaction isn't
being properly evaluated for all of the atoms within the direct-space cutoff,
contrary to what is specified in the Darden papers. That's all I'm worried
about here, Erik.

The actual equations don't allow for ANY atom pairs (inside/outside cutoff)
to be excluded, but the beta factor makes it possible to neglect the part of it
outside e.g. 9 Å. With the Ewald solution idefined by Eqs. 2.3-2.5 in the Essman95
paper, we're simply saying that the best numerical PME solution is the one that
minimizes the error in these equations, whether it's random or systematic.

Of course, it would be quite straightforward is to enable an alternative table
where the table data is truncated at the cutoff, and filled with zero outside. I'll
see what we can do about that for 4.0!

However, there is also a deeper technical reason for our choice of cutoff treatment.
Classically it's a good idea to check if an interaction is outside a cutoff, since you can skip the rest of the code on that case (done e.g. in Amber/Charmm).
However, already with SSE we use 4-way SIMD, and although I can do the check and set one of the four registers to zero I still have to do the other interactions and don't save anything, unless all 4 are outside the cutoff. This is going to get even worse (well, better performance ;-) in 2009/2010 when Intel's Sandy Bridge CPUs arrive with 256-bit wide AVX registers, and we're facing similar issues on GPUs when using neighborsearching.

So, the big question in that case is if we should start wasting more and more cycles, or try to use them the best way we can (at least by default)? It might be 10-20% extra cost to extend the neighborlist with SSE loops, but when we have PME working on GPUs it might be 100%!

On a sidenote, I think this is a valuable discussion, but what I'd really like is a set of higher-level tests to assess the accuracy of free energies, deviation from "exact" trajectories, cRMS of small proteins, etc., since that would make it much easier to compare pros and cons of the different approaches, in particular if we're interested e.g. in "getting a free energy as accurate as possible in X CPU hours".

Cheers,

Erik

#8 Updated by Berk Hess over 11 years ago

There are two different issues here.
One is the actual accuracy of the energy/force evaluation,
the other if things are defined in a unique way so one can
compare packages.

I think both discussions are largely academical.
ewald_rtol is by default 1e-5 at the cut-off, which means that
the error in the energy or force there is negligibly small
(unless you want near-perfect energy conservation over
long times). For two particles of charge e the energy
at 1 nm is 0.0014 kJ/mol.
In most cases other terms would probably cause larger
differences between different packages,
for instance the way LJ switching/shifting
and the dispersion correction is treated.

I understand that you want exact comparable energies.
If there was a simple way to do this, we would implement
it, but the cubic spline interpolation does not allow
for this. In principle we could go through the effort
to implement real abrupt cut-off's in assembly, but this
seems a waste of time to me, as this would make integration
less accurate.

The only way to currently get exactly what you want
is to use charge groups of single atoms, and set
rcoulomb=rlist with nstlist=1.

But again, we are talking about energies in the range
of 1e-3 kJ/mol.

Even we would have the option of exact cut-off, I would
only use if for comparison purposes. The electrostatics
gets more accurate when you include the interactions
between rcoulomb and rlist. The real electrostatics
is 1/r, not Darden's equations which also depend on the choice
of cut-off radius and ewald_rtol.

We should aim for the most accurate force evaluation,
not for the most accurate comparison with other packages.
However, it is useful to be able to do exact comparisons
with other packages. For PME this is possible
(as described above), but at a relatively high computational
cost.

Berk.

#9 Updated by Berk Hess over 11 years ago

I have now added a hint for PME-Switch in the grompp error message
when PME is set with rcoulomb<rlist.

Berk.

#10 Updated by John D. over 11 years ago

Berk, Erik:

Apologies for my absence from this thread -- I was in Banff for a conference and then traveling, and wasn't able to attend to my mail.

Berk: You're right, comparison with other codes is possible with nstlist = 1. That should suffice for purposes of comparison. You're also right that the errors will be small for a single interaction at the cutoff, though I haven't yet managed to convince myself that the errors will remain small after nine or so steps after a pairlist build -- this is something that gave us a good deal of trouble in the reaction-field case by inducing heating as waters drifted into the cutoff and their interactions were suddenly included when the next pairlist build occurred. It is something that can certainly be tested, however, which we should do. The modified warning in grompp is much appreciated!

Erik: Your Comment #7 is an excellent discussion of the real issues at hand here -- accuracy in terms of both error differences for single force evaluations (compared to including all atoms in the direct-space Ewald sum), but also in terms of computed properties of interest, like free energies of binding or hydration. Stephen Bond (UIUC) gave an excellent talk at the Banff meeting about how simulation parameters such as the choice of integrator and timestep can perturb the probability distribution being sampled, and how this has real consequences for computed expectations. I'm trying to put together a "regression suite" for free energy calculations (just hydration free energies right now) to examine their sensitivity to parameter choices, mostly to inform our own work and ensure our free energy protocols are robust with respect to parameters. But it would be nice to see this applied to later versions of gromacs as well.

I'll inform you when we have some interesting results here.

Thanks again to you both!

- John

Also available in: Atom PDF