I have started using CORESHELL package (thanks Hendrik Heenen) in LAMMPS. I have the following questions-
1. In the documentation, following is mentioned under the context of Coulomb interaction between a core and its shell: "Coulomb interaction between a core and its shell should be turned off using the special_bonds command ... ``This induces a long-range correction approximation which fails at small distances (~< 10e-8). Therefore, the Coulomb term which is used to calculate the correction factor is extended by a minimal distance (r_min = 1.0-6) when the interaction between a core/shell pair is treated."
What is "```long-range correction approximation"? I can understand that the approximated complementary error function’s accuracy is around 1.0e-7 and when the distance between a core and its shell ```~< 1.0e-7, short-range computations(force, energy) are erroneous.``
2. Is there a reason why Wolf summation is not implemented (considering that it might not involve lot of effort) alongside Ewald to treat Coulombic interactions? Or are there any complications that are not apparent?
3. Can someone point me to a source where the constants(used in the source code) B0, B1, ..., B5 are mentioned? Or how they have come to be? I was able to identify constants A0, ..., A5 that are used for approximating complementary error function in (normal) Ewald summation.
>What is "```long-range correction approximation"? I can understand that the approximated >complementary error function’s accuracy is around 1.0e-7 and when the distance between a >core and its shell ```~< 1.0e-7, short-range computations(force, energy) are erroneous.``
You are right. LAMMPS calculates the Ewald sum for all point charges per default. That means one has
to remove non-counted interactions (via special_bonds) from the forces afterwards, hence “correction”. The
existing error functions accuracy suffers from oscillations at distances smaller than 1.0e-6 - therefore the internal
cutoff is used.
>3. Can someone point me to a source where the constants(used in the source code) B0, B1, >..., B5 are mentioned? Or how they have come to be? I was able to identify constants A0, >..., A5 that are used for approximating complementary error function in (normal) Ewald >summation.
These constants have been proposed by Alain Dequidt and coworkers (authors of the USER-DRUDE package - CCd).
This parametrization is more accurate for small distances so one can use a smaller internal cutoff - as I remember the
constants A0-A5 already fail at 1.0e-3 which is a quite probable distance between a core and its attached shell.
Note that these precautions with internal cutoffs and new erfc function have to be taken as the core-shell distances
are so unphysically small!
>2. Is there a reason why Wolf summation is not implemented (considering that it might not >involve lot of effort) alongside Ewald to treat Coulombic interactions? Or are there any >complications that are not apparent?
No there is no apparent complication I can think of. We simply have not got around to doing it.
on the topic of computing the error function complement, i have one
comment. please see below.
Hello Vishal,
thanks for your email,
What is "long-range correction approximation"? I can understand that the
approximated >complementary error function's accuracy is around 1.0e-7 and
when the distance between a >core and its shell ~< 1.0e-7, short-range
computations(force, energy) are erroneous.
You are right. LAMMPS calculates the Ewald sum for all point charges per
default. That means one has
to remove non-counted interactions (via special_bonds) from the forces
afterwards, hence "correction". The
existing error functions accuracy suffers from oscillations at distances
smaller than 1.0e-6 - therefore the internal
cutoff is used.
3. Can someone point me to a source where the constants(used in the source
code) B0, B1, >..., B5 are mentioned? Or how they have come to be? I was
able to identify constants A0, >..., A5 that are used for approximating
complementary error function in (normal) Ewald >summation.
These constants have been proposed by Alain Dequidt and coworkers (authors
of the USER-DRUDE package - CCd).
This parametrization is more accurate for small distances so one can use a
smaller internal cutoff - as I remember the
constants A0-A5 already fail at 1.0e-3 which is a quite probable distance
between a core and its attached shell.
please note, that we are about to replace the (rather inaccurate) erfc
(or rather erfcx) approximation throughout the LAMMPS code with a
version that is derived from the one proposed in libcerf ( http://apps.jcns.fz-juelich.de/libcerf ) this has superior precision
and is - after a few tweaks - within less than 10% of the speed of the
5th degree polynomial approximation. with this, the lines:
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
u = 1.0 - t;
erfc = t * (1.+u*(B0+u*(B1+u*(B2+u*(B3+u*(B4+u*B5)))))) * expm2;
in my tests, this implementation has provided close to double
precision accuracy across the entire argument value range, and thus is
significantly superior to either of the two approximations mentioned
above. the functions are already available in the MathSpecial
namespace of the standard distribution; have not yet converted all
long-range coulomb styles to them. i'd like to encourage both the
CORESHELL and USER-DRUDE package developers to make a test with this
new version, as there is hope to consolidate a lot of replicated code
with this.
I’m using also the CORESHELL model and found myself needing a Wolf like approach (for double checking the code I’m using to develop the potentials), and also avoids forcing pbc to use ewald sum.
Actually, rather than the plain Wolf sum the shifted version introduced by Fennell and Gezelter called damped shifted force (DSF) is preferred, which avoids problems with derivative at cutoff.
This model is already implemented in LAMMPS pair styles coul/dsf and lj/cut/coul/dsf
Since the Core shell model generally consists of Coulomb + Buckingham, I created based on coul/dsf the born/coul/long the pair styles
born/coul/dsf and the core-shell version
born/coul/dsf/cs
I used the born one, since it’s more general and includes buckingham.
Taking note of axel’s advice I’m making use of the MathSpecial::my_erfcx function to avoid the polynomials for the error function.
By doing so, the trick of adding the extra EWALD_EPS is avoided, and the /cs version consists basically of just adding the small distance for the potential r=0 cases.
I’ve been testing with the Na-Cl example in the CORESHELL package, comparing the runs with the provided born/coul/long/cs and everything seems to be working fine.
If you want to test yourself this is the commit adding these pair styles to the latest lammps-icms branch as of today.
this is a response only to the part pertaining to submitting a pull request.
with the recent switch to officially support contributions via github
and merging them directly into the "master" branch, the "lammps-icms"
branch has become obsolete.
pull requests should now be submitted relative to "master". in this
particular case, this should not be a big issue, since these are only
new files. you may still submit it relative to lammps-icms, but then i
would "convert" it.
for inclusion into LAMMPS, you should add proper documentation to the
manual. in this case a few lines at the top and a paragraph or two in
the middle of doc/src/pair_cs.txt should do.
including a small/fast input deck (with reference output) for the
examples/coreshell folder would be appreciated. since this is a more
complex to set up model, it is generally helpful for users to see a
fully working input. if you can add some comments about why which
settings/parameters were chose, it would be awesome.
i will send you an invite as an external collaborator of the LAMMPS
project on github. that will allow us to assign issues or pull
requests to you concerning your contribution.
if you have questions about these steps or anything else related to
contributing your code, don't hesitate to respond.