piece-wise pair coul long potential

LAMMPS LAMMPS Mailing List Mirror
1

Dear Lammps user and developer,

I am currently working on modifying the pair_coul_long potential to describe the Coulomb interaction between 100 electrons and a large perovskite system. The goal is to add another cutoff distant to the pair_coul_long potential so that when the distance between an electron and perovskite atoms is smaller than cutoff1 (r<cutoff1 ), the E=q1q2/cutoff1 (constant) and F=0; while r>=cutoff1 and < cutoff2, E and F just follow the normal Coulomb interaction.

I am trying to achieve this in two ways:

The first one is modifying the pair_coul_long.cpp, by changing the PairCoulLong2::compute and PairCoulLong2::single subroutine, I can set the real space F=0 and E=E_cutoff1 when r<cutoff1. But I did not find an efficient way to modify the kspace energy and force.

I am also trying to let the real space F= - kspace F by defining a vector equals to “f2group” in the Ewlad.cpp and call this vector in the pair_coul_long.cpp:

Ewald::compute_group_group:

…
f2group[0] *= qscale;

f2group[1] *= qscale;

f2group[2] *= qscale;

yf2group[0] += f2group[0];

yf2group[1] += f2group[1];

yf2group[2] += f2group[2];

…
I also defined the yf2group vector in ksapce.h.
Then in PairCoulLong2::compute:
…
if (rsq < 9.61) {

r2inv = 1.0/9.61;

r = sqrt(9.61);

grij = g_ewald * r;

expm2 = exp(-grij*grij);

t = 1.0 / (1.0 + EWALD_P*grij);

erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;

prefactor = qqrd2e * scale[itype][jtype] * qtmp*q[j]/r;
fpair = 0;

f[i][0] -= yvector_kspace[0];

f[i][1] -= yvector_kspace[1];

f[i][2] -= yvector_kspace[2];

if (newton_pair || j < nlocal) {

f[j][0] += yvector_kspace[0];

f[j][1] += yvector_kspace[1];

f[j][2] += yvector_kspace[2];

}
…
else if (rsq >= 9.61 && rsq < cut_coulsq) {

r2inv = 1.0/rsq;

r = sqrt(rsq);

grij = g_ewald * r;

expm2 = exp(-grij*grij);

t = 1.0 / (1.0 + EWALD_P*grij);

erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;

prefactor = qqrd2e * scale[itype][jtype] * qtmp*q[j]/r;

forcecoul = prefactor * (erfc + EWALD_Fgrijexpm2);

…

But when I run the simulation, the kspace contribution is still there and when I echo the “yf2group”, it is 0. I think it might because the compute_group_group.cpp calculates the real space and kspace energy and force simultaneously, so I cannot let the real space force equals the negative force calculated by the kspace.cpp.

Is there any way to exclude the kspace force and energy when r smaller than cuoff1?

The other way is following the previous answer https://lammps.sandia.gov/threads/msg83130.html

you can use logical expressions in atom style and equal style expressions. here you need to use, of course, an atom style variable. those logicals would evaluate to either 1.0 (if true) or 0.0 (if false) and thus if you multiply those with the force contribution for a given section and sum each segment, you can define a piecewise force.

axel.

But I am not sure how to write the true or false command, like how to determine if the distance between each electron and each perovskite atoms is smaller than the cutoff1?

I am using lammps-7Aug19.

Does anyone have other suggestions than these two ways?

Thank you very much for your help!

Best,

Yang Hong

2

you cannot make a modification like this as this will - as you already noticed - create inconsistencies with kspace.
please note that kspace_styles will always consider all charged particles, since its very definition is based on all-to-all interactions across a lattice.

also an atom style variable based approach cannot work, since this is about pairs of atoms.

please note that whatever workaround you will find, it will implement an inconsistent coulomb behavior.
with that in mind, you could not set charges for the “electron particles” and instead write a pair style similar to coul/cut where you assign the charges used by the atom type (or only for specific atoms) and then include those interactions using pair style hybrid without affecting the (long-range) coulomb treatment for the rest of the system.

axel.

3

Hi Axel,

Thank you for your suggestions. I have made a pair_coul_cut potential and successfully applied it to the system.

But I am still curious about the coul/long potential. Is compute_group_group.cpp the right place to combine the real space and kspace energy and force, and send the net force to do the velocity verlet? If so, can I add a statement there such as if rsq < cutoff1sq, then net force = net force* 0? Something works like the “factor_coul”, it will not affect the kspace solver and will not be performed until both component forces are calculated.

Also, I did not find the “factor_coul” term in ewald.cpp or kspace.cpp, so how does the special_bond affect the kspace energy calculations?

Best,
Yang Hong

4

Hi Axel,

Thank you for your suggestions. I have made a pair_coul_cut potential and successfully applied it to the system.

But I am still curious about the coul/long potential. Is compute_group_group.cpp the right place to combine the real space and kspace energy and force, and send the net force to do the velocity verlet?

no. compute group/group is very inefficient and not used during the regular for computation. please read the information provided in the manual about modifying and extending LAMMPS. also, you may want to have a look at section 3 of https://lammps.sandia.gov/doc/Developer.pdf for additional background information.

If so, can I add a statement there such as if rsq < cutoff1sq, then net force = net force* 0? Something works like the “factor_coul”, it will not affect the kspace solver and will not be performed until both component forces are calculated.

no. compute group/group is a diagnostic tool that is used in addition to the regular (and much more efficient) force computation for analysis purposes only. please see the warning in the documentation for compute group/group how it can be costly (i.e. will slow your calculation down) when used too often.

Also, I did not find the “factor_coul” term in ewald.cpp or kspace.cpp, so how does the special_bond affect the kspace energy calculations?

as already explained, kspace can be rationalized as an interaction of each particle with an infinite lattice of the replica of all other charges on all other atoms (and itself). this is is not a pairwise interaction. exclusion factors, however, apply only to the nearest pair of two atoms, i.e. the pair that is connected through the bond topology as a 1-2, 1-3, or 1-4 pair. if you look closely at the code in the corresponding pair styles, you should recognize that for excluded pairs, the factor is not applied to the damped coulomb (like for the unscaled interactions), but rather to the full coulomb interaction for that specific pair. that way kspace does not need to know about it and all other interactions with the replica of the excluded atom are computed in full, as they need to be.

axel.