correct nofdotr handling for tersoff styles

This commit is contained in:
Axel Kohlmeyer
2021-07-11 00:55:22 -04:00
parent bfc9df1302
commit f82096c46c
8 changed files with 68 additions and 46 deletions

View File

@ -88,10 +88,10 @@ void PairTersoff::compute(int eflag, int vflag)
if (shift_flag) {
if (evflag) {
if (eflag) {
if (vflag_atom) eval<1,1,1,1>();
if (vflag_either) eval<1,1,1,1>();
else eval<1,1,1,0>();
} else {
if (vflag_atom) eval<1,1,0,1>();
if (vflag_either) eval<1,1,0,1>();
else eval<1,1,0,0>();
}
} else eval<1,0,0,0>();
@ -100,17 +100,17 @@ void PairTersoff::compute(int eflag, int vflag)
if (evflag) {
if (eflag) {
if (vflag_atom) eval<0,1,1,1>();
if (vflag_either) eval<0,1,1,1>();
else eval<0,1,1,0>();
} else {
if (vflag_atom) eval<0,1,0,1>();
if (vflag_either) eval<0,1,0,1>();
else eval<0,1,0,0>();
}
} else eval<0,0,0,0>();
}
}
template <int SHIFT_FLAG, int EVFLAG, int EFLAG, int VFLAG_ATOM>
template <int SHIFT_FLAG, int EVFLAG, int EFLAG, int VFLAG_EITHER>
void PairTersoff::eval()
{
int i,j,k,ii,jj,kk,inum,jnum;
@ -315,7 +315,7 @@ void PairTersoff::eval()
f[k][1] += fk[1];
f[k][2] += fk[2];
if (VFLAG_ATOM) v_tally3(i,j,k,fj,fk,delr1,delr2);
if (VFLAG_EITHER) v_tally3(i,j,k,fj,fk,delr1,delr2);
}
f[j][0] += fjxtmp;
f[j][1] += fjytmp;

View File

@ -60,20 +60,20 @@ void PairTersoffMODCOMP::compute(int eflag, int vflag)
if (shift_flag) {
if (evflag) {
if (eflag) {
if (vflag_atom) eval<1,1,1,1>(ifrom, ito, thr);
if (vflag_either) eval<1,1,1,1>(ifrom, ito, thr);
else eval<1,1,1,0>(ifrom, ito, thr);
} else {
if (vflag_atom) eval<1,1,0,1>(ifrom, ito, thr);
if (vflag_either) eval<1,1,0,1>(ifrom, ito, thr);
else eval<1,1,0,0>(ifrom, ito, thr);
}
} else eval<1,0,0,0>(ifrom, ito, thr);
} else {
if (evflag) {
if (eflag) {
if (vflag_atom) eval<0,1,1,1>(ifrom, ito, thr);
if (vflag_either) eval<0,1,1,1>(ifrom, ito, thr);
else eval<0,1,1,0>(ifrom, ito, thr);
} else {
if (vflag_atom) eval<0,1,0,1>(ifrom, ito, thr);
if (vflag_either) eval<0,1,0,1>(ifrom, ito, thr);
else eval<0,1,0,0>(ifrom, ito, thr);
}
} else eval<0,0,0,0>(ifrom, ito, thr);
@ -84,7 +84,7 @@ void PairTersoffMODCOMP::compute(int eflag, int vflag)
} // end of omp parallel region
}
template <int SHIFT_FLAG, int EVFLAG, int EFLAG, int VFLAG_ATOM>
template <int SHIFT_FLAG, int EVFLAG, int EFLAG, int VFLAG_EITHER>
void PairTersoffMODCOMP::eval(int iifrom, int iito, ThrData * const thr)
{
int i,j,k,ii,jj,kk,jnum;
@ -282,7 +282,7 @@ void PairTersoffMODCOMP::eval(int iifrom, int iito, ThrData * const thr)
f[k].y += fk[1];
f[k].z += fk[2];
if (VFLAG_ATOM) v_tally3_thr(i,j,k,fj,fk,delr1,delr2,thr);
if (VFLAG_EITHER) v_tally3_thr(this,i,j,k,fj,fk,delr1,delr2,thr);
}
f[j].x += fjxtmp;
f[j].y += fjytmp;

View File

@ -60,20 +60,20 @@ void PairTersoffMODOMP::compute(int eflag, int vflag)
if (shift_flag) {
if (evflag) {
if (eflag) {
if (vflag_atom) eval<1,1,1,1>(ifrom, ito, thr);
if (vflag_either) eval<1,1,1,1>(ifrom, ito, thr);
else eval<1,1,1,0>(ifrom, ito, thr);
} else {
if (vflag_atom) eval<1,1,0,1>(ifrom, ito, thr);
if (vflag_either) eval<1,1,0,1>(ifrom, ito, thr);
else eval<1,1,0,0>(ifrom, ito, thr);
}
} else eval<1,0,0,0>(ifrom, ito, thr);
} else {
if (evflag) {
if (eflag) {
if (vflag_atom) eval<0,1,1,1>(ifrom, ito, thr);
if (vflag_either) eval<0,1,1,1>(ifrom, ito, thr);
else eval<0,1,1,0>(ifrom, ito, thr);
} else {
if (vflag_atom) eval<0,1,0,1>(ifrom, ito, thr);
if (vflag_either) eval<0,1,0,1>(ifrom, ito, thr);
else eval<0,1,0,0>(ifrom, ito, thr);
}
} else eval<0,0,0,0>(ifrom, ito, thr);
@ -84,7 +84,7 @@ void PairTersoffMODOMP::compute(int eflag, int vflag)
} // end of omp parallel region
}
template <int SHIFT_FLAG, int EVFLAG, int EFLAG, int VFLAG_ATOM>
template <int SHIFT_FLAG, int EVFLAG, int EFLAG, int VFLAG_EITHER>
void PairTersoffMODOMP::eval(int iifrom, int iito, ThrData * const thr)
{
int i,j,k,ii,jj,kk,jnum;
@ -282,7 +282,7 @@ void PairTersoffMODOMP::eval(int iifrom, int iito, ThrData * const thr)
f[k].y += fk[1];
f[k].z += fk[2];
if (VFLAG_ATOM) v_tally3_thr(i,j,k,fj,fk,delr1,delr2,thr);
if (VFLAG_EITHER) v_tally3_thr(this,i,j,k,fj,fk,delr1,delr2,thr);
}
f[j].x += fjxtmp;
f[j].y += fjytmp;

View File

@ -61,10 +61,10 @@ void PairTersoffOMP::compute(int eflag, int vflag)
if (shift_flag) {
if (evflag) {
if (eflag) {
if (vflag_atom) eval<1,1,1,1>(ifrom, ito, thr);
if (vflag_either) eval<1,1,1,1>(ifrom, ito, thr);
else eval<1,1,1,0>(ifrom, ito, thr);
} else {
if (vflag_atom) eval<1,1,0,1>(ifrom, ito, thr);
if (vflag_either) eval<1,1,0,1>(ifrom, ito, thr);
else eval<1,1,0,0>(ifrom, ito, thr);
}
} else eval<1,0,0,0>(ifrom, ito, thr);
@ -73,10 +73,10 @@ void PairTersoffOMP::compute(int eflag, int vflag)
if (evflag) {
if (eflag) {
if (vflag_atom) eval<0,1,1,1>(ifrom, ito, thr);
if (vflag_either) eval<0,1,1,1>(ifrom, ito, thr);
else eval<0,1,1,0>(ifrom, ito, thr);
} else {
if (vflag_atom) eval<0,1,0,1>(ifrom, ito, thr);
if (vflag_either) eval<0,1,0,1>(ifrom, ito, thr);
else eval<0,1,0,0>(ifrom, ito, thr);
}
} else eval<0,0,0,0>(ifrom, ito, thr);
@ -87,7 +87,7 @@ void PairTersoffOMP::compute(int eflag, int vflag)
} // end of omp parallel region
}
template <int SHIFT_FLAG, int EVFLAG, int EFLAG, int VFLAG_ATOM>
template <int SHIFT_FLAG, int EVFLAG, int EFLAG, int VFLAG_EITHER>
void PairTersoffOMP::eval(int iifrom, int iito, ThrData * const thr)
{
int i,j,k,ii,jj,kk,jnum,maxshort_thr;
@ -293,7 +293,7 @@ void PairTersoffOMP::eval(int iifrom, int iito, ThrData * const thr)
f[k].y += fk[1];
f[k].z += fk[2];
if (VFLAG_ATOM) v_tally3_thr(i,j,k,fj,fk,delr1,delr2,thr);
if (VFLAG_EITHER) v_tally3_thr(this,i,j,k,fj,fk,delr1,delr2,thr);
}
f[j].x += fjxtmp;
f[j].y += fjytmp;

View File

@ -80,7 +80,7 @@ void PairTersoffTableOMP::compute(int eflag, int vflag)
ev_setup_thr(eflag, vflag, nall, eatom, vatom, nullptr, thr);
if (evflag)
if (vflag_atom) eval<1,1>(ifrom, ito, thr);
if (vflag_either) eval<1,1>(ifrom, ito, thr);
else eval<1,0>(ifrom, ito, thr);
else eval<0,0>(ifrom, ito, thr);
@ -89,7 +89,7 @@ void PairTersoffTableOMP::compute(int eflag, int vflag)
} // end of omp parallel region
}
template <int EVFLAG, int VFLAG_ATOM>
template <int EVFLAG, int VFLAG_EITHER>
void PairTersoffTableOMP::eval(int iifrom, int iito, ThrData * const thr)
{
int i,j,k,ii,jnum;
@ -427,7 +427,7 @@ void PairTersoffTableOMP::eval(int iifrom, int iito, ThrData * const thr)
fytmp += f_ij[1] + f_ik[1];
fztmp += f_ij[2] + f_ik[2];
if (VFLAG_ATOM) v_tally3_thr(i,j,k,f_ij,f_ik,dr_ij,dr_ik,thr);
if (VFLAG_EITHER) v_tally3_thr(this,i,j,k,f_ij,f_ik,dr_ij,dr_ik,thr);
}
// second loop over neighbors of atom i except j, forces and virial only - part 2/2
@ -493,7 +493,7 @@ void PairTersoffTableOMP::eval(int iifrom, int iito, ThrData * const thr)
fytmp += f_ij[1] + f_ik[1];
fztmp += f_ij[2] + f_ik[2];
if (VFLAG_ATOM) v_tally3_thr(i,j,k,f_ij,f_ik,dr_ij,dr_ik,thr);
if (VFLAG_EITHER) v_tally3_thr(this,i,j,k,f_ij,f_ik,dr_ij,dr_ik,thr);
}
} // loop on J

View File

@ -1486,23 +1486,34 @@ void ThrOMP::v_tally2_thr(const int i, const int j, const double fpair,
called by AIREBO and Tersoff potential, newton_pair is always on
------------------------------------------------------------------------- */
void ThrOMP::v_tally3_thr(const int i, const int j, const int k,
void ThrOMP::v_tally3_thr(Pair *pair, const int i, const int j, const int k,
const double * const fi, const double * const fj,
const double * const drik, const double * const drjk,
ThrData * const thr)
{
double v[6];
v[0] = THIRD * (drik[0]*fi[0] + drjk[0]*fj[0]);
v[1] = THIRD * (drik[1]*fi[1] + drjk[1]*fj[1]);
v[2] = THIRD * (drik[2]*fi[2] + drjk[2]*fj[2]);
v[3] = THIRD * (drik[0]*fi[1] + drjk[0]*fj[1]);
v[4] = THIRD * (drik[0]*fi[2] + drjk[0]*fj[2]);
v[5] = THIRD * (drik[1]*fi[2] + drjk[1]*fj[2]);
v[0] = (drik[0]*fi[0] + drjk[0]*fj[0]);
v[1] = (drik[1]*fi[1] + drjk[1]*fj[1]);
v[2] = (drik[2]*fi[2] + drjk[2]*fj[2]);
v[3] = (drik[0]*fi[1] + drjk[0]*fj[1]);
v[4] = (drik[0]*fi[2] + drjk[0]*fj[2]);
v[5] = (drik[1]*fi[2] + drjk[1]*fj[2]);
v_tally(thr->vatom_pair[i],v);
v_tally(thr->vatom_pair[j],v);
v_tally(thr->vatom_pair[k],v);
if (pair->vflag_global) v_tally(thr->virial_pair,v);
if (pair->vflag_atom) {
v[0] *= THIRD;
v[1] *= THIRD;
v[2] *= THIRD;
v[3] *= THIRD;
v[4] *= THIRD;
v[5] *= THIRD;
v_tally(thr->vatom_pair[i],v);
v_tally(thr->vatom_pair[j],v);
v_tally(thr->vatom_pair[k],v);
}
}
/* ----------------------------------------------------------------------

View File

@ -138,6 +138,8 @@ class ThrOMP {
void ev_tally3_thr(Pair *const, const int, const int, const int, const double, const double,
const double *const, const double *const, const double *const,
const double *const, ThrData *const);
void v_tally3_thr(Pair *const, const int, const int, const int, const double *const,
const double *const, const double *const, const double *const, ThrData *const);
void ev_tally4_thr(Pair *const, const int, const int, const int, const int, const double,
const double *const, const double *const, const double *const,
const double *const, const double *const, const double *const, ThrData *const);
@ -170,8 +172,6 @@ class ThrOMP {
// style independent versions
void v_tally2_thr(const int, const int, const double, const double *const, ThrData *const);
void v_tally3_thr(const int, const int, const int, const double *const, const double *const,
const double *const, const double *const, ThrData *const);
void v_tally4_thr(const int, const int, const int, const int, const double *const,
const double *const, const double *const, const double *const,
const double *const, const double *const, ThrData *const);

View File

@ -1529,12 +1529,23 @@ void Pair::v_tally3(int i, int j, int k,
v[4] = THIRD * (drik[0]*fi[2] + drjk[0]*fj[2]);
v[5] = THIRD * (drik[1]*fi[2] + drjk[1]*fj[2]);
vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2];
vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5];
vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2];
vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5];
vatom[k][0] += v[0]; vatom[k][1] += v[1]; vatom[k][2] += v[2];
vatom[k][3] += v[3]; vatom[k][4] += v[4]; vatom[k][5] += v[5];
if (vflag_global) {
virial[0] += 3.0*v[0];
virial[1] += 3.0*v[1];
virial[2] += 3.0*v[2];
virial[3] += 3.0*v[3];
virial[4] += 3.0*v[4];
virial[5] += 3.0*v[5];
}
if (vflag_atom) {
vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2];
vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5];
vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2];
vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5];
vatom[k][0] += v[0]; vatom[k][1] += v[1]; vatom[k][2] += v[2];
vatom[k][3] += v[3]; vatom[k][4] += v[4]; vatom[k][5] += v[5];
}
}
/* ----------------------------------------------------------------------