Dear Lammps users and developers,
I combined the “Tersoff.cpp” and “lj/cut/coul/dsf.cpp” potentials by modifying the pair_Tersoff.cpp file and adding lj/coul to the original Tersoff file. I was wondering
1- What is “neighbor->request(this)” in init_style() function of PairLJCutCoulDSF? and What is its difference with the following lines in PairTersoff::init_style() ?
int irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
2- I combined these two potentials but I am not sure which of the above methods I should use for combined potential. Which one I have to use? neighbor->request(this) or the lines in Tersoff file?
3- The init_one() function in Tersoff potential returns “cutmax” parameter while the one in lj/cut/coul/dsf returns “cut” parameter which is the maximum of cutoffs for Coulomb and lj interactions. For combination of these two potentials, Should I use the maximum of all three cutoffs? i.e., cut_lj, cut_coul, and cutoff for Tersoff (R+D).
/////////////////////// Tersoff
double PairTersoff::init_one(int i, int j)
{
if (setflag[i][j] == 0) error->all(FLERR,“All pair coeffs are not set”);
return cutmax;
}
/////////////////////// Lj/Coul/dsf
double PairLJCutCoulDSF::init_one(int i, int j)
{
double cut = MAX(cut_lj[i][j],cut_coul);
cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j];
…
…
…
return cut;
}
Thank you very much
Hossein
Dear Lammps users and developers,
I combined the "Tersoff.cpp" and "lj/cut/coul/dsf.cpp" potentials by
modifying the pair_Tersoff.cpp file and adding lj/coul to the original
Tersoff file. I was wondering
why on earth would you need to do such a thing, when you could just
use pair style hybrid/overlay and be done with it?
1- What is "neighbor->request(this)" in init_style() function of
PairLJCutCoulDSF? and What is its difference with the following lines in
PairTersoff::init_style() ?
int irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
both are requests to compute a neighbor list. the one in
PairLJCutCoulDSF asks for a regular "half" neighbor list (every pair
of local atoms is listed only once, with newton pair on, also for
pairs between subdomains)., PairTersoff instead changes the request to
ask for a "full" neighbor list (every pair of atoms is listed).
2- I combined these two potentials but I am not sure which of the above
methods I should use for combined potential. Which one I have to use?
neighbor->request(this) or the lines in Tersoff file?
you need to use the full neighbor list, but you have to adapt tjhe
llj/cut/coul/dsf code accordingly
3- The init_one() function in Tersoff potential returns "cutmax" parameter
while the one in lj/cut/coul/dsf returns "cut" parameter which is the
maximum of cutoffs for Coulomb and lj interactions. For combination of these
two potentials, Should I use the maximum of all three cutoffs? i.e., cut_lj,
cut_coul, and cutoff for Tersoff (R+D).
yes, because that influences what cutoff is used when constructing the
requested neighbor lists and for communication between subdomains.
Thank you very much Axel. You are right, I can use pair hybrid/overlay but for my purpose I need to combine these potentials in order to do some heavy calculations after force calculation.
As you said, I used full neighbor list and removed the highlighted lines from lj/cut/coul/dsf code
fpair = (forcecoul + factor_ljforcelj) * r2inv;
f[i][0] += delxfpair;
f[i][1] += delyfpair;
f[i][2] += delzfpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delxfpair;
f[j][1] -= delyfpair;
f[j][2] -= delz*fpair;
}
By adapting lj/cut/coul/dsf code, did you mean removing the highlighted lines? or I need to do some other modifications?
After applying above modifications, the calculated force on each atom doesn’t match with the one using pair_hybrid/overlay. I am sure that the code works well for each potential(if I active the force calculation for that potential only) but there is something wrong with neighbor list or cut off probably. Do you have more suggestion?
I also have some questions about pair_hybrid/overlay
1- I think combined potential should work faster compared to pair_hybrid/overlay because in pair_hybrid/overlay the forces are calculated for each individual potential and then they are added at each time step. therefore, it should take longer time to do force calculation compared to combined tersoff + lj/coul potential. Am I right?
2- If I use pair_hybrid/overlay, full neighbor list is requested for Tersoff while for lj/coul half neighbor list will be requested. Therefore, do you think that is the reason that forces from my combined potential don’t match with the ones with pair_hybrid/overlay?
thanks
Hossein
Thank you very much Axel. You are right, I can use pair hybrid/overlay but
for my purpose I need to combine these potentials in order to do some heavy
calculations after force calculation.
As you said, I used full neighbor list and removed the highlighted lines
from lj/cut/coul/dsf code
fpair = (forcecoul + factor_lj*forcelj) * r2inv;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
}
By adapting lj/cut/coul/dsf code, did you mean removing the highlighted
lines? or I need to do some other modifications?
it is impossible to say without knowing how exactly you merged the two
compute routines.
the tersoff potential already has code that does compute its pairwise
additive force component without double counting from a full neighbor
list.
After applying above modifications, the calculated force on each atom
doesn't match with the one using pair_hybrid/overlay. I am sure that the
code works well for each potential(if I active the force calculation for
that potential only) but there is something wrong with neighbor list or cut
off probably. Do you have more suggestion?
read and understand more of the code. if you are doing modifications
to the LAMMPS code that go this deep, you better fully understand what
LAMMPS is doing rather than guessing or depending on the say-so of
others. since i cannot look over your shoulder i would have to guess
as well and thus cannot give much more advice beyond that.
I also have some questions about pair_hybrid/overlay
1- I think combined potential should work faster compared to
pair_hybrid/overlay because in pair_hybrid/overlay the forces are calculated
for each individual potential and then they are added at each time step.
therefore, it should take longer time to do force calculation compared to
combined tersoff + lj/coul potential. Am I right?
it depends. not automatically. if you had merged the two pair styles
correctly and adjust the force tally for lj/cut/coul/dsf as above, you
would compute the force for each pair twice and then tally it only
into the force of one atom. that encompasses more computations. also,
lj/cut/coul/dsf would require a much longer pairwise cutoff and thus
your neighborlist would have many more neighbors per atom and thus
when computing the tersoff force terms that go beyond the pairwise
additive component, you would have to compute more distances. the
latter can become quite expensive since you have a 3-fold nested loop
and not only a 2-fold nested loop. however, without making some tests
and checking the output, i cannot say for certain, that with
hybrid/overlay the neighbor lists are suitably trimmed. it may require
to use the "multi" style neighbor list generation flag, rather than
the default "bin".
2- If I use pair_hybrid/overlay, full neighbor list is requested for Tersoff
while for lj/coul half neighbor list will be requested. Therefore, do you
think that is the reason that forces from my combined potential don't match
with the ones with pair_hybrid/overlay?
no.
thanks
Hossein
p.s.: a more general question. are you certain that what you are
computing here is a meaningful model? i would be extremely concerned
about double counting short range forces. that is why conventional
force fields use exclusions, so that for pairs where short range
interactions are computed via bond/angle/dihedral potentials, some or
all of the non-bonded contribution is skipped. AIREBO (which includes
a LJ style long(er) range dispersion component) uses switching
functions to smoothly transition from its (bonded-like) short(er)
range interactions to the long(er) range dispersion interactions from
the LJ potential.
Thank you very much Axel.