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

avoid division by zero when using cutoff 0.0 with pair_modify shift yes

Browse Source
This commit is contained in:
Axel Kohlmeyer
2017-07-14 23:33:00 -04:00
parent 3f297382ac
commit 4ec07422f0
42 changed files with 97 additions and 88 deletions
Download Patch File Download Diff File Expand all files Collapse all files

94
src/USER-MISC/pair_list.cpp
View File

@ -81,7 +81,7 @@ void PairList::compute(int eflag, int vflag)
{
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = eflag_global =
vflag_global = eflag_atom = vflag_atom = 0;
vflag_global = eflag_atom = vflag_atom = 0;
const int nlocal = atom->nlocal;
const int newton_pair = force->newton_pair;
@ -126,46 +126,46 @@ void PairList::compute(int eflag, int vflag)
const double r2inv = 1.0/rsq;
if (style[n] == HARM) {
const double r = sqrt(rsq);
const double dr = par.parm.harm.r0 - r;
fpair = 2.0*par.parm.harm.k*dr/r;
const double r = sqrt(rsq);
const double dr = par.parm.harm.r0 - r;
fpair = 2.0*par.parm.harm.k*dr/r;
if (eflag_either)
epair = par.parm.harm.k*dr*dr - par.offset;
if (eflag_either)
epair = par.parm.harm.k*dr*dr - par.offset;
} else if (style[n] == MORSE) {
const double r = sqrt(rsq);
const double dr = par.parm.morse.r0 - r;
const double dexp = exp(par.parm.morse.alpha * dr);
fpair = 2.0*par.parm.morse.d0*par.parm.morse.alpha
* (dexp*dexp - dexp) / r;
const double r = sqrt(rsq);
const double dr = par.parm.morse.r0 - r;
const double dexp = exp(par.parm.morse.alpha * dr);
fpair = 2.0*par.parm.morse.d0*par.parm.morse.alpha
* (dexp*dexp - dexp) / r;
if (eflag_either)
epair = par.parm.morse.d0 * (dexp*dexp - 2.0*dexp) - par.offset;
if (eflag_either)
epair = par.parm.morse.d0 * (dexp*dexp - 2.0*dexp) - par.offset;
} else if (style[n] == LJ126) {
const double r6inv = r2inv*r2inv*r2inv;
const double sig6 = mypow(par.parm.lj126.sigma,6);
fpair = 24.0*par.parm.lj126.epsilon*r6inv
* (2.0*sig6*sig6*r6inv - sig6) * r2inv;
const double r6inv = r2inv*r2inv*r2inv;
const double sig6 = mypow(par.parm.lj126.sigma,6);
fpair = 24.0*par.parm.lj126.epsilon*r6inv
* (2.0*sig6*sig6*r6inv - sig6) * r2inv;
if (eflag_either)
epair = 4.0*par.parm.lj126.epsilon*r6inv
* (sig6*sig6*r6inv - sig6) - par.offset;
if (eflag_either)
epair = 4.0*par.parm.lj126.epsilon*r6inv
* (sig6*sig6*r6inv - sig6) - par.offset;
}
if (newton_pair || i < nlocal) {
f[i].x += dx*fpair;
f[i].y += dy*fpair;
f[i].z += dz*fpair;
f[i].x += dx*fpair;
f[i].y += dy*fpair;
f[i].z += dz*fpair;
}
if (newton_pair || j < nlocal) {
f[j].x -= dx*fpair;
f[j].y -= dy*fpair;
f[j].z -= dz*fpair;
f[j].x -= dx*fpair;
f[j].y -= dy*fpair;
f[j].z -= dz*fpair;
}
if (evflag) ev_tally(i,j,nlocal,newton_pair,epair,0.0,fpair,dx,dy,dz);
@ -174,10 +174,10 @@ void PairList::compute(int eflag, int vflag)
if (vflag_fdotr) virial_fdotr_compute();
if (check_flag) {
int tmp;
MPI_Allreduce(&pc,&tmp,1,MPI_INT,MPI_SUM,world);
if (tmp != 2*npairs)
error->all(FLERR,"Not all pairs processed in pair_style list");
int tmp;
MPI_Allreduce(&pc,&tmp,1,MPI_INT,MPI_SUM,world);
if (tmp != 2*npairs)
error->all(FLERR,"Not all pairs processed in pair_style list");
}
}
@ -263,12 +263,12 @@ void PairList::settings(int narg, char **arg)
ptr = strtok(NULL," \t\n\r\f");
if ((ptr == NULL) || (*ptr == '#'))
error->all(FLERR,"Incorrectly formatted harmonic pair parameters");
error->all(FLERR,"Incorrectly formatted harmonic pair parameters");
par.parm.harm.k = force->numeric(FLERR,ptr);
ptr = strtok(NULL," \t\n\r\f");
if ((ptr == NULL) || (*ptr == '#'))
error->all(FLERR,"Incorrectly formatted harmonic pair parameters");
error->all(FLERR,"Incorrectly formatted harmonic pair parameters");
par.parm.harm.r0 = force->numeric(FLERR,ptr);
++nharm;
@ -279,17 +279,17 @@ void PairList::settings(int narg, char **arg)
ptr = strtok(NULL," \t\n\r\f");
if (!ptr)
error->all(FLERR,"Incorrectly formatted morse pair parameters");
error->all(FLERR,"Incorrectly formatted morse pair parameters");
par.parm.morse.d0 = force->numeric(FLERR,ptr);
ptr = strtok(NULL," \t\n\r\f");
if (!ptr)
error->all(FLERR,"Incorrectly formatted morse pair parameters");
error->all(FLERR,"Incorrectly formatted morse pair parameters");
par.parm.morse.alpha = force->numeric(FLERR,ptr);
ptr = strtok(NULL," \t\n\r\f");
if (!ptr)
error->all(FLERR,"Incorrectly formatted morse pair parameters");
error->all(FLERR,"Incorrectly formatted morse pair parameters");
par.parm.morse.r0 = force->numeric(FLERR,ptr);
++nmorse;
@ -300,12 +300,12 @@ void PairList::settings(int narg, char **arg)
ptr = strtok(NULL," \t\n\r\f");
if (!ptr)
error->all(FLERR,"Incorrectly formatted 12-6 LJ pair parameters");
error->all(FLERR,"Incorrectly formatted 12-6 LJ pair parameters");
par.parm.lj126.epsilon = force->numeric(FLERR,ptr);
ptr = strtok(NULL," \t\n\r\f");
if (!ptr)
error->all(FLERR,"Incorrectly formatted 12-6 LJ pair parameters");
error->all(FLERR,"Incorrectly formatted 12-6 LJ pair parameters");
par.parm.lj126.sigma = force->numeric(FLERR,ptr);
++nlj126;
@ -332,10 +332,10 @@ void PairList::settings(int narg, char **arg)
if (comm->me == 0) {
if (screen)
fprintf(screen,"Read %d (%d/%d/%d) interacting pairs from %s\n",
npairs, nharm, nmorse, nlj126, arg[0]);
npairs, nharm, nmorse, nlj126, arg[0]);
if (logfile)
fprintf(logfile,"Read %d (%d/%d/%d) interacting pairs from %s\n",
npairs, nharm, nmorse, nlj126, arg[0]);
npairs, nharm, nmorse, nlj126, arg[0]);
}
}
@ -380,18 +380,18 @@ void PairList::init_style()
list_parm_t &par = params[n];
if (style[n] == HARM) {
const double dr = sqrt(par.cutsq) - par.parm.harm.r0;
par.offset = par.parm.harm.k*dr*dr;
const double dr = sqrt(par.cutsq) - par.parm.harm.r0;
par.offset = par.parm.harm.k*dr*dr;
} else if (style[n] == MORSE) {
const double dr = par.parm.morse.r0 - sqrt(par.cutsq);
const double dexp = exp(par.parm.morse.alpha * dr);
par.offset = par.parm.morse.d0 * (dexp*dexp - 2.0*dexp);
const double dr = par.parm.morse.r0 - sqrt(par.cutsq);
const double dexp = exp(par.parm.morse.alpha * dr);
par.offset = par.parm.morse.d0 * (dexp*dexp - 2.0*dexp);
} else if (style[n] == LJ126) {
const double r6inv = par.cutsq*par.cutsq*par.cutsq;
const double sig6 = mypow(par.parm.lj126.sigma,6);
par.offset = 4.0*par.parm.lj126.epsilon*r6inv * (sig6*sig6*r6inv - sig6);
const double r6inv = par.cutsq*par.cutsq*par.cutsq;
const double sig6 = mypow(par.parm.lj126.sigma,6);
par.offset = 4.0*par.parm.lj126.epsilon*r6inv * (sig6*sig6*r6inv - sig6);
}
}
}