diff --git a/src/KSPACE/pair_born_coul_long.cpp b/src/KSPACE/pair_born_coul_long.cpp index ad45200000..721dc287de 100644 --- a/src/KSPACE/pair_born_coul_long.cpp +++ b/src/KSPACE/pair_born_coul_long.cpp @@ -291,9 +291,9 @@ double PairBornCoulLong::init_one(int i, int j) born3[i][j] = 8.0*d[i][j]; if (offset_flag) { - double rexp = exp(-cut_lj[i][j]/rho[i][j]); - offset[i][j] = a[i][j]*rexp - c[i][j]/pow(cut_lj[i][j],6.0) - + d[i][j]/pow(cut_lj[i][j],8.0); + double rexp = exp(-cut_lj[i][j]*rhoinv[i][j]); + offset[i][j] = a[i][j]*rexp - c[i][j]/pow(cut_lj[i][j],6.0) + + d[i][j]/pow(cut_lj[i][j],8.0); } else offset[i][j] = 0.0; cut_ljsq[j][i] = cut_ljsq[i][j]; @@ -307,6 +307,39 @@ double PairBornCoulLong::init_one(int i, int j) born3[j][i] = born3[i][j]; offset[j][i] = offset[i][j]; + // compute I,J contribution to long-range tail correction + // count total # of atoms of type I and J via Allreduce + + if (tail_flag) { + int *type = atom->type; + int nlocal = atom->nlocal; + + double count[2],all[2]; + count[0] = count[1] = 0.0; + for (int k = 0; k < nlocal; k++) { + if (type[k] == i) count[0] += 1.0; + if (type[k] == j) count[1] += 1.0; + } + MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world); + + double PI = 4.0*atan(1.0); + double rho1 = rho[i][j]; + double rho2 = rho1*rho1; + double rho3 = rho2*rho1; + double rc = cut_lj[i][j]; + double rc2 = rc*rc; + double rc3 = rc2*rc; + double rc5 = rc3*rc2; + etail_ij = 2.0*PI*all[0]*all[1] * + (a[i][j]*exp((sigma[i][j]-rc)/rho1)*rho1* + (rc2 + 2.0*rho1*rc + 2.0*rho2) - + c[i][j]/(3.0*rc3) + d[i][j]/(5.0*rc5)); + ptail_ij = (-1/3.0)*2.0*PI*all[0]*all[1] * + (-a[i][j]*exp((sigma[i][j]-rc)/rho1) * + (rc3 + 3.0*rho1*rc2 + 6.0*rho2*rc + 6.0*rho3) + + 2.0*c[i][j]/rc3 - 8.0*d[i][j]/(5.0*rc5)); + } + return cut; } @@ -442,9 +475,9 @@ double PairBornCoulLong::single(int i, int j, int itype, int jtype, if (rsq < cut_ljsq[itype][jtype]) { r6inv = r2inv*r2inv*r2inv; r = sqrt(rsq); - rexp = exp(-r*rhoinv[itype][jtype]); - forceborn = born1[itype][jtype]*r*rexp - born2[itype][jtype]*r6inv - + born3[itype][jtype]*r2inv*r6inv; + rexp = exp((sigma[itype][jtype]-r)*rhoinv[itype][jtype]); + forceborn = born1[itype][jtype]*r*rexp - born2[itype][jtype]*r6inv + + born3[itype][jtype]*r2inv*r6inv; } else forceborn = 0.0; fforce = (forcecoul + factor_lj*forceborn) * r2inv; @@ -455,8 +488,8 @@ double PairBornCoulLong::single(int i, int j, int itype, int jtype, eng += phicoul; } if (rsq < cut_ljsq[itype][jtype]) { - phiborn = a[itype][jtype]*rexp - c[itype][jtype]*r6inv - + d[itype][jtype]*r2inv*r6inv - offset[itype][jtype]; + phiborn = a[itype][jtype]*rexp - c[itype][jtype]*r6inv + + d[itype][jtype]*r2inv*r6inv - offset[itype][jtype]; eng += factor_lj*phiborn; } return eng; @@ -469,4 +502,3 @@ void *PairBornCoulLong::extract(char *str) if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul; return NULL; } - diff --git a/src/pair_born.cpp b/src/pair_born.cpp index 5c8ab293cb..9fa68e743b 100644 --- a/src/pair_born.cpp +++ b/src/pair_born.cpp @@ -64,7 +64,7 @@ void PairBorn::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; - double rsq,r2inv,r6inv,force_tf,factor_lj; + double rsq,r2inv,r6inv,forceborn,factor_lj; double r,rexp; int *ilist,*jlist,*numneigh,**firstneigh; @@ -247,9 +247,9 @@ double PairBorn::init_one(int i, int j) born3[i][j] = 8.0*d[i][j]; if (offset_flag) { - double rexp = exp((sigma[i][j]-cut_lj[i][j])*rhoinv[i][j]); - offset[i][j] = a[i][j]*rexp - c[i][j]/pow(cut_lj[i][j],6.0) + - d[i][j]/pow(cut_lj[i][j],8.0); + double rexp = exp((sigma[i][j]-cut[i][j])*rhoinv[i][j]); + offset[i][j] = a[i][j]*rexp - c[i][j]/pow(cut[i][j],6.0) + + d[i][j]/pow(cut[i][j],8.0); } else offset[i][j] = 0.0; a[j][i] = a[i][j]; @@ -281,15 +281,16 @@ double PairBorn::init_one(int i, int j) double rho1 = rho[i][j]; double rho2 = rho1*rho1; double rho3 = rho2*rho1; - double rc = sigma[i][j] - cut[i][j]; + double rc = cut[i][j]; double rc2 = rc*rc; double rc3 = rc2*rc; double rc5 = rc3*rc2; etail_ij = 2.0*PI*all[0]*all[1] * - (a[i][j]*exp(-rc/rho1)*rho1* (rc2 + 2.0*rho1*rc + 2.0*rho2) - + (a[i][j]*exp((sigma[i][j]-rc)/rho1)*rho1* + (rc2 + 2.0*rho1*rc + 2.0*rho2) - c[i][j]/(3.0*rc3) + d[i][j]/(5.0*rc5)); ptail_ij = (-1/3.0)*2.0*PI*all[0]*all[1] * - (-a[i][j]*exp(-rc/rho1) * + (-a[i][j]*exp((sigma[i][j]-rc)/rho1) * (rc3 + 3.0*rho1*rc2 + 6.0*rho2*rc + 6.0*rho3) + 2.0*c[i][j]/rc3 - 8.0*d[i][j]/(5.0*rc5)); } @@ -396,7 +397,7 @@ double PairBorn::single(int i, int j, int itype, int jtype, rexp = exp((sigma[itype][jtype]-r)*rhoinv[itype][jtype]); forceborn = born1[itype][jtype]*r*rexp - born2[itype][jtype]*r6inv + born3[itype][jtype]*r2inv*r6inv; - fforce = factor_lj*force_born*r2inv; + fforce = factor_lj*forceborn*r2inv; phiborn = a[itype][jtype]*rexp - c[itype][jtype]*r6inv + d[itype][jtype]*r2inv*r6inv - offset[itype][jtype]; diff --git a/src/pair_born.h b/src/pair_born.h index a51cf0fba1..8b62e47f12 100644 --- a/src/pair_born.h +++ b/src/pair_born.h @@ -42,7 +42,6 @@ class PairBorn : public Pair { double cut_global; double **cut; double **a,**rho,**sigma,**c, **d; - double **rhoinv,**tf1,**tf2,**tf3,**offset; double **rhoinv,**born1,**born2,**born3,**offset; void allocate(); diff --git a/src/pair_gauss.cpp b/src/pair_gauss.cpp index c9b2b18897..dd02501056 100644 --- a/src/pair_gauss.cpp +++ b/src/pair_gauss.cpp @@ -12,7 +12,7 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing author: Sai Jayaraman (University of Notre Dame) + Contributing author: Sai Jayaraman (U Notre Dame) ------------------------------------------------------------------------- */ #include "math.h" @@ -104,7 +104,6 @@ void PairGauss::compute(int eflag, int vflag) r = sqrt(rsq); forcelj = - 2.0*a[itype][jtype]*b[itype][jtype] * rsq * exp(-b[itype][jtype]*rsq); - if (r < 0.00001) forcelj = 0.0; fpair = forcelj*r2inv; f[i][0] += delx*fpair; @@ -116,11 +115,9 @@ void PairGauss::compute(int eflag, int vflag) f[j][2] -= delz*fpair; } - if (eflag) { + if (eflag) evdwl = -(a[itype][jtype]*exp(-b[itype][jtype]*rsq) - offset[itype][jtype]); - evdwl *= MAX(EPSILON,latscale); - } if (evflag) ev_tally(i,j,nlocal,newton_pair, evdwl,0.0,fpair,delx,dely,delz); @@ -312,7 +309,6 @@ double PairGauss::single(int i, int j, int itype, int jtype, double rsq, r2inv = 1.0/rsq; philj = -(a[itype][jtype]*exp(-b[itype][jtype]*rsq) - offset[itype][jtype]); - philj *= MAX(eps,latscale); forcelj = -2.0*a[itype][jtype]*b[itype][jtype]*rsq*exp(-b[itype][jtype]*rsq); fforce = forcelj*r2inv;