Might be a typo in pair_coul_debye

LAMMPS LAMMPS Mailing List Mirror
1

Hi,

I am trying to modify based on pair coul/debye, its potential is

_images/pair_debye.jpg

But maybe the following red r2inv should be rinv? It might be a mistake or I understand any place wrongly?

double PairCoulDebye::single(int i, int j, int itype, int jtype,
double rsq, double factor_coul, double factor_lj,
double &fforce)
{
double r2inv,r,rinv,forcecoul,phicoul,screening;

r2inv = 1.0/rsq;
r = sqrt(rsq);
rinv = 1.0/r;
screening = exp(-kappa*r);
forcecoul = force->qqrd2e * atom->q[i]atom->q[j] *
screening * (kappa + rinv);
fforce = factor_coulforcecoul * r2inv;

phicoul = force->qqrd2e * atom->q[i]atom->q[j] * rinv * screening;
return factor_coulphicoul;
}

Best Regards,

NIE PIN(Ms)

Nanyang Technology University, School of Physical and Mathematical Science

Tel: +65 84522947

email: PNIE001@…2775…

2

Hi,

I am trying to modify based on pair coul/debye, its potential is

[image: _images/pair_debye.jpg]
But maybe the following red r2inv should be rinv? It might be a mistake
or I understand any place wrongly?

​no, it may be confusing, but it is correct for the way how Pair::single()
is used. to get the actual force vector, the value of fforce needs to be
multiplied by deltax, deltay, and deltaz, respectively. please see e.g. the
implementation of compute group/group for example, that uses Pair::single()

it is consistent with other pair styles, e.g. you can easily compare to,
e.g. pair coul/cut by setting kappa to 0.0 and similarly to lj/cut/coul/cut
by also setting epsilon to 0.0

axel.