Hi all,

I’m trying out a screened coulombic pair style similar to coul/debye (in my case, using the complementary error function for attenuation), and I noticed a line in pair_coul_debye.cpp that struck me as a bit odd. It’s very possible that I’m mis-reading the code, but in case I’m misunderstanding the implementation style (or in the off chance that it’s a bug), I thought I’d ask about it. In both PairCoulDebye::compute and PairCoulDebye::single, the coulomb’s law force is calculated as F=((k/epsilon)*q1*q2*screening) * (kappa + 1/r) * (1/r**2)*

with screening = exp(-kappar)

where k is the coulomb’s law constant, q1 and q2 are the respective charges on each atom in the pair, epsilon is the dielectric constant, kappa is the debye length, and r is the distance between the interacting atoms.

To be clear and not speak too abstractly, in the former instance, at lines 96-97 in PairCoulDebye::compute, this is coded as

96: forcecoul = qqrd2e * qtmp*q[j] * screening * (kappa + rinv);*

97: fpair = factor_coulforcecoul * r2inv;

Now, since F=-dE/dr here, and E=(k/e)*q1*q2 * (1/r) * exp(-kappa*r), I’d expect to calculate force as*

F = (k/e)*q1*q2 * (1/r**2) * exp(-kappar) + (k/e)*q1*q2*(1/r) * kappa*exp(-kappa*r)

or, reorganized as it is in the code,

F = (k/e)*q1*q1 * exp(-kappa*r) * (kappa + 1/r) * (1/r)*

In other words, I’d expect line 97 to read

fpair = factor_coulforcecoul * rinv;

and similarly, on line 184,

fforce = factor_coul*forcecoul * rinv;

Can anyone who’s looked at and used this code before confirm this? Am I somehow just completely screwing up my chain rule here, or should that r2inv on lines 97 and 184 actually be rinv?

Thanks for having a look at this. Hope all those who celebrated it had a great thanksgiving.

Tyler

Ack! I was indeed misunderstanding the implementation. I’d naively assumed that fpair and fforce were actually the magnitudes of force and didn’t look carefully enough at how the directional components of force were computed. Sorry about bothering everybody with my misunderstanding!

Tyler

hi tyler,

Hi all,

I'm trying out a screened coulombic pair style similar to coul/debye (in my

case, using the complementary error function for attenuation), and I noticed

a line in pair_coul_debye.cpp that struck me as a bit odd. It's very

possible that I'm mis-reading the code, but in case I'm misunderstanding the

implementation style (or in the off chance that it's a bug), I thought I'd

ask about it. In both PairCoulDebye::compute and PairCoulDebye::single, the

coulomb's law force is calculated as F=((k/epsilon)*q1*q2*screening) *

(kappa + 1/r) * (1/r**2)

with screening = exp(-kappa*r)

where k is the coulomb's law constant, q1 and q2 are the respective charges

on each atom in the pair, epsilon is the dielectric constant, kappa is the

debye length, and r is the distance between the interacting atoms.

To be clear and not speak too abstractly, in the former instance, at lines

96-97 in PairCoulDebye::compute, this is coded as

96: forcecoul = qqrd2e * qtmp*q[j] * screening * (kappa + rinv);

97: fpair = factor_coul*forcecoul * r2inv;

Now, since F=-dE/dr here, and E=(k/e)*q1*q2 * (1/r) * exp(-kappa*r), I'd

expect to calculate force as

F = (k/e)*q1*q2 * (1/r**2) * exp(-kappa*r) + (k/e)*q1*q2*(1/r) *

kappa*exp(-kappa*r)

or, reorganized as it is in the code,

F = (k/e)*q1*q1 * exp(-kappa*r) * (kappa + 1/r) * (1/r)

In other words, I'd expect line 97 to read

fpair = factor_coul*forcecoul * rinv;

and similarly, on line 184,

fforce = factor_coul*forcecoul * rinv;

Can anyone who's looked at and used this code before confirm this? Am I

somehow just completely screwing up my chain rule here, or should that r2inv

on lines 97 and 184 actually be rinv?

you are overlooking one little detail.

in line 99 the forces are incremented into the x-direction as

f[i][x] += delx * fpair;

if you change the code as you propose, then this

could would have to be:

f[i][x] += delx*rinv * fpair;

so, the fforce part already contains the *rinv from later down.

this doesn't make so much sense for coulomb, but makes

a lot of sense when doing 12-6 lennard-jones, as this little

change allows to avoid the sqrt(rsq) completely.

hope that helps,

axel.

Excellent point–sorry that I didn’t catch myself quick enough! Good to know the motivation for this form, too.

Tyler