pair_lj_cut_coul_long's treatment of damping factor

Dear list,

I am reading the codes of pair style which is to be coupled with ewald summation. From my understanding, the short-range term of electrostatic energy could be found by that of cut style multiplied with an erfc function, but I found out in pair_lj_cut_coul_long.cpp that this part is divided into two parts: one is by a power-to-the-fifth monomial, when r^2 is less than 2, and the other is by multiplying a table which is related to coulombic force calculation. I would like to know why erfc function, which is simple to use, is changed to this. Is it something related to efficiency?

Looking forward to hearing from you.

Andrew

Paul can give you details, but yes it is a faster way to compute erfc().

Steve

Paul can give you details, but yes it is a faster way to compute erfc().

it is a polynomial approximation of erfc()

the tabulation uses the same approximation to fill the tables and then
a very fast way to extract the tabulation.

axel.

Yes, both the polynomial and the tabulation are approximations used for computational efficiency.

Paul

Hi Paul,

I noticed in manual of pair_modify

"For N = 0, forces and energies are computed directly, using a polynomial fit for the needed erfc() function evaluation, which is what earlier versions of LAMMPS did. "

So I thought it was changed recently since erfc() is included in C++11 already. However I just checked 6Aug14 version and the polynomial approximation is still used. What is “earlier versions” mean here?

Sincerely,
Zhenxing

We’re talking about three ways of doing the same damped coulomb calculation for long-range electrostatics:

  1. erfc()

  2. polynomial fit of erfc()

  3. tabulation

LAMMPS uses all three, but in different places in the code. Before the tabulation code was written about a decade ago, LAMMPS did only (2). That’s what that statement is referring to. The ability to do exclusively (2) is still available in LAMMPS, by using “pair_modify table 0”.

The only place erfc() is used, AFAIK, is in building the tabulation. It is the most expensive, but the most accurate. The tabulation is the fastest and least accurate, but the accuracy is usually adequate and is tunable, so that is what LAMMPS mostly uses.

The polynomial fit is still used when R^2 is less than 2, which is typically relatively rare.

Paul

Hi Paul,

I noticed in manual of pair_modify

"For N = 0, forces and energies are computed directly, using a polynomial
fit for the needed erfc() function evaluation, which is what earlier
versions of LAMMPS did. "

So I thought it was changed recently since erfc() is included in C++11

there are two issues here: a) this sentence refers to using tabulation
(which is done in single precision floating point math, btw) vs.
direct evaluation (using the polynomial approximation). b) C++ 11
standard level is not relevant, since LAMMPS aims to be more portable.
C++ 11 is barely available in the latest C++ compiler versions.
however, most environments that LAMMPS has to run on, predate this.
also LAMMPS is mostly C with classes using the C runtime environment
and not so much the C++ runtime. that being said, erfc() is part of
C-99 standard and LAMMPS *does* depend on this for a few years.

already. However I just checked 6Aug14 version and the polynomial
approximation is still used. What is "earlier versions" mean here?

waaaay back, before tabulation for the short-range part of coulomb was
available.

i would also like to point out that the argument of erfc() that is
used in typical MD calculations with LAMMPS will be in a very limited
range (and the argument is always positive) and for that the
approximation is pretty good. please check out and compare to:

http://en.wikipedia.org/wiki/Error_function#Approximation_with_elementary_functions

you will have much larger errors, e.g. from having a cutoff on your LJ
interactions, or not using a very tight convergence parameter on the
long range solver. the often used 1.e-4 convergence for ewald
summation is not very tight and is usually the major cause for errors
on long-range coulomb forces. if you care about accuracy, you want to
set this to 1e-6 or smaller.

axel.