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

need to explicitly clamp the forces and energy in single function to zero at (implicit) pairwise cutoff

Browse Source
This commit is contained in:
Axel Kohlmeyer
2020-08-01 14:05:04 -04:00
parent 891be9313a
commit 4c3dc9566c
1 changed files with 12 additions and 7 deletions
Download Patch File Download Diff File Expand all files Collapse all files

19
src/pair_lj_cubic.cpp
View File

@ -129,8 +129,7 @@ void PairLJCubic::compute(int eflag, int vflag)
if (rsq <= cut_inner_sq[itype][jtype])
evdwl = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
else
evdwl = epsilon[itype][jtype]*
(PHIS + DPHIDS*t - A3*t*t*t/6.0);
evdwl = epsilon[itype][jtype]*(PHIS + DPHIDS*t - A3*t*t*t/6.0);
evdwl *= factor_lj;
if (evflag) ev_tally(i,j,nlocal,newton_pair,
@ -322,14 +321,21 @@ void PairLJCubic::read_restart_settings(FILE *fp)
/* ---------------------------------------------------------------------- */
double PairLJCubic::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 r2inv,r6inv,forcelj,philj;
double r,t;
double rmin;
// this is a truncated potential with an implicit cutoff
if (rsq >= cutsq[itype][jtype]) {
fforce=0.0;
return 0.0;
}
r2inv = 1.0/rsq;
if (rsq <= cut_inner_sq[itype][jtype]) {
r6inv = r2inv*r2inv*r2inv;
@ -345,8 +351,7 @@ double PairLJCubic::single(int /*i*/, int /*j*/, int itype, int jtype,
if (rsq <= cut_inner_sq[itype][jtype])
philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
else
philj = epsilon[itype][jtype]*
(PHIS + DPHIDS*t - A3*t*t*t/6.0);
philj = epsilon[itype][jtype]*(PHIS + DPHIDS*t - A3*t*t*t/6.0);
return factor_lj*philj;
}