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

bugfix for exclusion problem in coul/dsf reported by nathan houtz and peter wirnsberger

Browse Source
This commit is contained in:
Axel Kohlmeyer
2015-03-14 08:54:44 -04:00
parent 96b60c7aea
commit 2215d34587
3 changed files with 14 additions and 11 deletions
Download Patch File Download Diff File Expand all files Collapse all files

10
src/pair_coul_dsf.cpp
View File

@ -65,7 +65,7 @@ void PairCoulDSF::compute(int eflag, int vflag)
int i,j,ii,jj,inum,jnum;
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;
double r,rsq,r2inv,forcecoul,factor_coul;
double prefactor,erfcc,erfcd,e_self,t;
double prefactor,erfcc,erfcd,t;
int *ilist,*jlist,*numneigh,**firstneigh;
ecoul = 0.0;
@ -97,7 +97,7 @@ void PairCoulDSF::compute(int eflag, int vflag)
jnum = numneigh[i];
if (eflag) {
e_self = -(e_shift/2.0 + alpha/MY_PIS) * qtmp*qtmp*qqrd2e;
double e_self = -(e_shift/2.0 + alpha/MY_PIS) * qtmp*qtmp*qqrd2e;
ev_tally(i,i,nlocal,0,0.0,e_self,0.0,0.0,0.0,0.0);
}
@ -115,14 +115,15 @@ void PairCoulDSF::compute(int eflag, int vflag)
r2inv = 1.0/rsq;
r = sqrt(rsq);
prefactor = factor_coul * qqrd2e*qtmp*q[j]/r;
prefactor = qqrd2e*qtmp*q[j]/r;
erfcd = exp(-alpha*alpha*rsq);
t = 1.0 / (1.0 + EWALD_P*alpha*r);
erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd;
forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS * erfcd +
r*f_shift) * r;
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
fpair = forcecoul * r2inv;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
@ -134,6 +135,7 @@ void PairCoulDSF::compute(int eflag, int vflag)
if (eflag) {
ecoul = prefactor * (erfcc - r*e_shift - rsq*f_shift);
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
} else ecoul = 0.0;
if (evflag) ev_tally(i,j,nlocal,newton_pair,

8
src/pair_coul_wolf.cpp
View File

@ -131,11 +131,9 @@ void PairCoulWolf::compute(int eflag, int vflag)
}
if (eflag) {
if (rsq < cut_coulsq) {
ecoul = v_sh;
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
} else ecoul = 0.0;
}
ecoul = v_sh;
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
} else ecoul = 0.0;
if (evflag) ev_tally(i,j,nlocal,newton_pair,
0.0,ecoul,fpair,delx,dely,delz);

7
src/pair_lj_cut_coul_dsf.cpp
View File

@ -136,12 +136,14 @@ void PairLJCutCoulDSF::compute(int eflag, int vflag)
if (rsq < cut_coulsq) {
r = sqrt(rsq);
prefactor = factor_coul * qqrd2e*qtmp*q[j]/r;
prefactor = qqrd2e*qtmp*q[j]/r;
erfcd = exp(-alpha*alpha*r*r);
t = 1.0 / (1.0 + EWALD_P*alpha*r);
erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd;
forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS * erfcd +
r*f_shift) * r;
r*f_shift) * r;
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
fpair = forcecoul * r2inv;
} else forcecoul = 0.0;
fpair = (forcecoul + factor_lj*forcelj) * r2inv;
@ -163,6 +165,7 @@ void PairLJCutCoulDSF::compute(int eflag, int vflag)
if (rsq < cut_coulsq) {
ecoul = prefactor * (erfcc - r*e_shift - rsq*f_shift);
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
} else ecoul = 0.0;
}