ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

avoid sqrt computation on cutoffs. make single functions consistent with compute

Browse Source
This commit is contained in:
Axel Kohlmeyer
2020-08-07 18:34:10 -04:00
parent 9f469623c0
commit 17b7476217
1 changed files with 28 additions and 21 deletions
Download Patch File Download Diff File Expand all files Collapse all files

49
src/USER-MISC/pair_coul_shield.cpp
View File

@ -115,8 +115,8 @@ void PairCoulShield::compute(int eflag, int vflag)
// turn on/off taper function
if (tap_flag) {
Tap = calc_Tap(r,sqrt(cutsq[itype][jtype]));
dTap = calc_dTap(r,sqrt(cutsq[itype][jtype]));
Tap = calc_Tap(r,cut[itype][jtype]);
dTap = calc_dTap(r,cut[itype][jtype]);
} else {Tap = 1.0; dTap = 0.0;}
forcecoul = qqrd2e*qtmp*q[j]*r*depsdr;
@ -340,34 +340,41 @@ void PairCoulShield::read_restart_settings(FILE *fp)
/* ---------------------------------------------------------------------- */
double PairCoulShield::single(int i, int j, int itype, int jtype,
double rsq, double factor_coul, double /*factor_lj*/,
double &fforce)
double rsq, double factor_coul, double /*factor_lj*/,
double &fforce)
{
double r, rarg,Vc,fvc,forcecoul,phishieldec;
double r3,th,epsr,depsdr,Tap,dTap;
double *q = atom->q;
double qqrd2e = force->qqrd2e;
r = sqrt(rsq);
r3 = rsq*r;
rarg = 1.0/sigmae[itype][jtype];
th = r3 + MathSpecial::cube(rarg);
epsr = 1.0/pow(th,0.333333333333333333);
depsdr = epsr*epsr;
depsdr *= depsdr;
Vc = qqrd2e*q[i]*q[j]*epsr;
// only computed between different layers as indicated by different molecule ids.
// turn on/off taper function
if (tap_flag) {
Tap = calc_Tap(r,sqrt(cutsq[itype][jtype]));
dTap = calc_dTap(r,sqrt(cutsq[itype][jtype]));
} else {Tap = 1.0; dTap = 0.0;}
if (atom->molecule[i] == atom->molecule[j]) {
fforce = 0.0;
return 0.0;
}
forcecoul = qqrd2e*q[i]*q[j]*r*depsdr;
fvc = forcecoul*Tap - Vc*dTap/r;
fforce = factor_coul*fvc;
r = sqrt(rsq);
r3 = rsq*r;
rarg = 1.0/sigmae[itype][jtype];
th = r3 + MathSpecial::cube(rarg);
epsr = 1.0/pow(th,0.333333333333333333);
depsdr = epsr*epsr;
depsdr *= depsdr;
Vc = qqrd2e*q[i]*q[j]*epsr;
if (tap_flag) phishieldec = factor_coul*Vc*Tap;
// turn on/off taper function
if (tap_flag) {
Tap = calc_Tap(r,cut[itype][jtype]);
dTap = calc_dTap(r,cut[itype][jtype]);
} else {Tap = 1.0; dTap = 0.0;}
forcecoul = qqrd2e*q[i]*q[j]*r*depsdr;
fvc = forcecoul*Tap - Vc*dTap/r;
fforce = factor_coul*fvc;
if (tap_flag) phishieldec = Vc*Tap;
else phishieldec = Vc - offset[itype][jtype];
return factor_coul*phishieldec;
}