use a half neighbor list for skip_threebody case for further speedup

This commit is contained in:
Axel Kohlmeyer
2022-07-02 17:39:01 -04:00
parent f38a417a32
commit 2fa5f7c97e
4 changed files with 42 additions and 24 deletions

View File

@ -1101,7 +1101,11 @@ void PairSWIntel::allocate()
void PairSWIntel::init_style()
{
// there is no support for skipping threebody loops (yet)
bool tmp_threebody = skip_threebody_flag;
skip_threebody_flag = false;
PairSW::init_style();
skip_threebody_flag = tmp_threebody;
map[0] = map[1];

View File

@ -392,7 +392,11 @@ void PairSWKokkos<DeviceType>::coeff(int narg, char **arg)
template<class DeviceType>
void PairSWKokkos<DeviceType>::init_style()
{
// there is no support for skipping threebody loops (yet)
bool tmp_threebody = skip_threebody_flag;
skip_threebody_flag = false;
PairSW::init_style();
skip_threebody_flag = tmp_threebody;
// adjust neighbor list request for KOKKOS

View File

@ -14,6 +14,7 @@
/* ----------------------------------------------------------------------
Contributing author: Aidan Thompson (SNL)
Optimizations for two-body only: Jackson Elowitt (Univ. of Utah)
------------------------------------------------------------------------- */
#include "pair_sw.h"
@ -138,14 +139,18 @@ void PairSW::compute(int eflag, int vflag)
}
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
// only need to skip if we have a full neighbor list
if (!skip_threebody_flag) {
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
}
}
twobody(&params[ijparam],rsq,fpair,eflag,evdwl);
@ -273,9 +278,12 @@ void PairSW::init_style()
if (force->newton_pair == 0)
error->all(FLERR,"Pair style Stillinger-Weber requires newton pair on");
// need a full neighbor list
// need a full neighbor list for full threebody calculation
neighbor->add_request(this, NeighConst::REQ_FULL);
if (skip_threebody_flag)
neighbor->add_request(this);
else
neighbor->add_request(this, NeighConst::REQ_FULL);
}
/* ----------------------------------------------------------------------

View File

@ -134,14 +134,16 @@ void PairSWOMP::eval(int iifrom, int iito, ThrData * const thr)
}
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (x[j].z < ztmp) continue;
if (x[j].z == ztmp && x[j].y < ytmp) continue;
if (x[j].z == ztmp && x[j].y == ytmp && x[j].x < xtmp) continue;
if (!skip_threebody_flag) {
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (x[j].z < ztmp) continue;
if (x[j].z == ztmp && x[j].y < ytmp) continue;
if (x[j].z == ztmp && x[j].y == ytmp && x[j].x < xtmp) continue;
}
}
twobody(&params[ijparam],rsq,fpair,EFLAG,evdwl);
@ -169,24 +171,24 @@ void PairSWOMP::eval(int iifrom, int iito, ThrData * const thr)
delr1[1] = x[j].y - ytmp;
delr1[2] = x[j].z - ztmp;
rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
double fjxtmp,fjytmp,fjztmp;
fjxtmp = fjytmp = fjztmp = 0.0;
for (kk = jj+1; kk < numshort; kk++) {
k = neighshort_thr[kk];
ktype = map[type[k]];
ikparam = elem3param[itype][ktype][ktype];
ijkparam = elem3param[itype][jtype][ktype];
delr2[0] = x[k].x - xtmp;
delr2[1] = x[k].y - ytmp;
delr2[2] = x[k].z - ztmp;
rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
threebody(&params[ijparam],&params[ikparam],&params[ijkparam],
rsq1,rsq2,delr1,delr2,fj,fk,EFLAG,evdwl);
fxtmp -= fj[0] + fk[0];
fytmp -= fj[1] + fk[1];
fztmp -= fj[2] + fk[2];
@ -196,7 +198,7 @@ void PairSWOMP::eval(int iifrom, int iito, ThrData * const thr)
f[k].x += fk[0];
f[k].y += fk[1];
f[k].z += fk[2];
if (EVFLAG) ev_tally3_thr(this,i,j,k,evdwl,0.0,fj,fk,delr1,delr2,thr);
}
f[j].x += fjxtmp;