Convergence test | kspace precision

LAMMPS LAMMPS Mailing List Mirror
1

Dear LAMMPS developers,

I am trying to compute Coulomb energies up to a high precision for
testing purposes.

The system I consider is very simple: one SPC/E water molecule in a
cubic box (please find the input script and data file attached). I am
only interested in the sum ecoul+elong at timestep 0.

Here are the results for kspace ewald with a precision ranging from
10^-6 to 10^-14 (version 25Feb16):

-0.17174293258645
-0.17170600041624
-0.1716601235286
-0.17167025873511
-0.17167495005634
-0.17167500295002
-0.17167126200061
-0.1716649212963
-0.17165706647864

These are the values I got with "pair_modify table 0" in addition:

-0.17173894259253
-0.17170520561771
-0.17166206555518
-0.17167364960872
-0.17167846349162
-0.17167762652758
-0.17167238340016
-0.17166428487184
-0.17165468286819

These energy values do not seem change significantly upon increasing the
precision by 8 orders of magnitude (the absolute error is about 10^-5
kcal/mol).

I wonder why this is the case. Is it possible that this behaviour is due
to some pre-calculated constants in pair_lj_cut_coul_long.cpp which are
only accurate up to, say 8-9 digits (e.g. EWALD_F, EWALD_P, A1-A5)?

Peter

spce.dat (723 Bytes)

energy.lmp (1.06 KB)

2

Dear LAMMPS developers,

I am trying to compute Coulomb energies up to a high precision for
testing purposes.

The system I consider is very simple: one SPC/E water molecule in a
cubic box (please find the input script and data file attached). I am
only interested in the sum ecoul+elong at timestep 0.

Here are the results for kspace ewald with a precision ranging from
10^-6 to 10^-14 (version 25Feb16):

-0.17174293258645
-0.17170600041624
-0.1716601235286
-0.17167025873511
-0.17167495005634
-0.17167500295002
-0.17167126200061
-0.1716649212963
-0.17165706647864

These are the values I got with "pair_modify table 0" in addition:

-0.17173894259253
-0.17170520561771
-0.17166206555518
-0.17167364960872
-0.17167846349162
-0.17167762652758
-0.17167238340016
-0.17166428487184
-0.17165468286819

These energy values do not seem change significantly upon increasing the
precision by 8 orders of magnitude (the absolute error is about 10^-5
kcal/mol).

I wonder why this is the case. Is it possible that this behaviour is due
to some pre-calculated constants in pair_lj_cut_coul_long.cpp which are
only accurate up to, say 8-9 digits (e.g. EWALD_F, EWALD_P, A1-A5)?

this is quite likely. those are part of an approximation for erfc(x),
that looks to me good for single precision floating point accuracy
only.
however, after paul crozier's latest refactoring of the coulomb tables
erfc() from the math library is used to fill the tables. but then the
tabulation is only in single precision...

originally, erfc(x) wasn't part of the C/C++ standard for a while and
thus not available on all machines. it was adopted only in ANSI-C99,
and not in ANSI-C89. by now this is probably a non-issue (since nobody
complained).

the problem now is more that using the math library functions from the
math library bundled with glibc often incurs a massive performance
penalty due to floating point mode switching. so it may still be
advantageous to have a custom (and inlining compatible) implementation
of erfc() in LAMMPS (and potentially also exp(), log(), and pow() with
a special case for pow() for integer powers, which is already
available as MathSpecial::powint() in math_special.h) as those offer
3-10x speedup for those functions over glibc's implementation.

BTW: for people using intel compilers, please note that one of reasons
for its improvements over gcc is that intel ships with its own version
of libm, which does not suffer from glibc's paranoia.

axel.

3

These energy values do not seem change significantly upon increasing
the precision by 8 orders of magnitude (the absolute error is about
10^-5 kcal/mol).

What were you expecting? The sum is converging, so the values aren't going to change significantly upon increasing the precision. How are you calculating the absolute error? If you are referring to the output in LAMMPS, that is only an estimate. The only way to get the true absolute error is to compare to a very, very high precision Ewald calculation, and that means increasing *both* the real space cutoff and the kspace precision to ridiculously high values so you know that it is converged.

Remember that the error is composed of two parts, real and reciprocal space. In theory you can hold the real space cutoff constant and still increase the accuracy by increasing the kspace precision, but that is only true if you have the correct g_ewald parameter. The kspace error estimator is based on approximations and it isn't perfect, especially if you are going to extremes, so you may also need to set g_ewald manually.

If you really think the problem is the approximation for erfc, then you could always just try it out the built-in erfc function... (see below) :slight_smile:

Stan

Index: KSPACE/pair_lj_cut_coul_long.cpp