The units of "rinv" and "sh_ewald" are inconsistent.
The following expression in "gromacs/gmxlib/nonbonded/nb_generic.cpp":
rinvcorr = (fr->coulomb_modifier == eintmodPOTSHIFT) ? rinv-fr->ic->sh_ewald : rinv;
Unit of "rinv" is inverse of nanometer;
"sh_ewald" is just a number without unit in that
"ic->sh_ewald = gmx_erfc(ic->ewaldcoeff_q*ic->rcoulomb);"
is defined in
Fix value of Ewald shift
In all the short-ranged kernel flavours, sh_ewald is subtracted from
rinv, which have inconsistent dimensions. Fortunately, rcutoff is
often close to 1, and the inconsistent shifts often cancel in
practice, and energy differences computed on neighbour lists of the
same size will have the error cancel. The difference doesn't even
show up in the regressiontests, but would if we had a unit test
of a single interaction.
#1 Updated by Mark Abraham about 3 years ago
There seems to be a minor, but embarrassing, issue here. sh_ewald is used as if it was computed as erfc(beta * cutoff)/ cutoff but computed only as erfc(beta * cutoff). Fortunately, cutoff is approximately one, and the errors from opposite charges approximately cancel, and for a given cutoff energy differences will have even more fortuitous cancelling.
I think we should subtract the cutoff distance from sh_ewald in forcerec.cpp. And plan to use a name that is more descriptive of the quantity, ie that it is a potential shift (albeit in 1/r form).
This issue was introduced along with the Verlet scheme, and now affects the group scheme also, in every Ewald kernel that I have seen.
#4 Updated by Berk Hess about 3 years ago
- Category set to documentation
- Status changed from New to Fix uploaded
Note that this bug only affects the reported energies, not the forces or the sampling.
And, as Mark noted, the effects on the energy are very small, so this is relatively harmless (but still embarrassing).