jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Propoaged SW cutoff fix into similar code elsewhere 

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14118 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
athomps 2015-10-15 18:41:59 +00:00
parent c7170a4296
commit 834731ddb4
4 changed files with 16 additions and 80 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
src/KOKKOS/pair_sw_kokkos.cpp
View File

@ -229,7 +229,7 @@ void PairSWKokkos<DeviceType>::operator()(TagPairSWComputeHalf<NEIGHFLAG,EVFLAG>
const F_FLOAT rsq = delx*delx + dely*dely + delz*delz; const F_FLOAT rsq = delx*delx + dely*dely + delz*delz;
const int ijparam = d_elem2param(itype,jtype,jtype); const int ijparam = d_elem2param(itype,jtype,jtype);
if (rsq > d_params[ijparam].cutsq) continue; if (rsq >= d_params[ijparam].cutsq) continue;
twobody(d_params[ijparam],rsq,fpair,eflag,evdwl); twobody(d_params[ijparam],rsq,fpair,eflag,evdwl);
@ -257,7 +257,7 @@ void PairSWKokkos<DeviceType>::operator()(TagPairSWComputeHalf<NEIGHFLAG,EVFLAG>
delr1[1] = x(j,1) - ytmp; delr1[1] = x(j,1) - ytmp;
delr1[2] = x(j,2) - ztmp; delr1[2] = x(j,2) - ztmp;
const F_FLOAT rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2]; const F_FLOAT rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
if (rsq1 > d_params[ijparam].cutsq) continue; if (rsq1 >= d_params[ijparam].cutsq) continue;
F_FLOAT fxtmpj = 0.0; F_FLOAT fxtmpj = 0.0;
F_FLOAT fytmpj = 0.0; F_FLOAT fytmpj = 0.0;
@ -275,7 +275,7 @@ void PairSWKokkos<DeviceType>::operator()(TagPairSWComputeHalf<NEIGHFLAG,EVFLAG>
delr2[2] = x(k,2) - ztmp; delr2[2] = x(k,2) - ztmp;
const F_FLOAT rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; const F_FLOAT rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
if (rsq2 > d_params[ikparam].cutsq) continue; if (rsq2 >= d_params[ikparam].cutsq) continue;
threebody(d_params[ijparam],d_params[ikparam],d_params[ijkparam], threebody(d_params[ijparam],d_params[ikparam],d_params[ijkparam],
rsq1,rsq2,delr1,delr2,fj,fk,eflag,evdwl); rsq1,rsq2,delr1,delr2,fj,fk,eflag,evdwl);
@ -355,7 +355,7 @@ void PairSWKokkos<DeviceType>::operator()(TagPairSWComputeFullA<NEIGHFLAG,EVFLAG
const int ijparam = d_elem2param(itype,jtype,jtype); const int ijparam = d_elem2param(itype,jtype,jtype);
if (rsq > d_params[ijparam].cutsq) continue; if (rsq >= d_params[ijparam].cutsq) continue;
twobody(d_params[ijparam],rsq,fpair,eflag,evdwl); twobody(d_params[ijparam],rsq,fpair,eflag,evdwl);
@ -381,7 +381,7 @@ void PairSWKokkos<DeviceType>::operator()(TagPairSWComputeFullA<NEIGHFLAG,EVFLAG
delr1[2] = x(j,2) - ztmp; delr1[2] = x(j,2) - ztmp;
const F_FLOAT rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2]; const F_FLOAT rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
if (rsq1 > d_params[ijparam].cutsq) continue; if (rsq1 >= d_params[ijparam].cutsq) continue;
for (int kk = jj+1; kk < jnum; kk++) { for (int kk = jj+1; kk < jnum; kk++) {
int k = d_neighbors(i,kk); int k = d_neighbors(i,kk);
@ -395,7 +395,7 @@ void PairSWKokkos<DeviceType>::operator()(TagPairSWComputeFullA<NEIGHFLAG,EVFLAG
delr2[2] = x(k,2) - ztmp; delr2[2] = x(k,2) - ztmp;
const F_FLOAT rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; const F_FLOAT rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
if (rsq2 > d_params[ikparam].cutsq) continue; if (rsq2 >= d_params[ikparam].cutsq) continue;
threebody(d_params[ijparam],d_params[ikparam],d_params[ijkparam], threebody(d_params[ijparam],d_params[ikparam],d_params[ijkparam],
rsq1,rsq2,delr1,delr2,fj,fk,eflag,evdwl); rsq1,rsq2,delr1,delr2,fj,fk,eflag,evdwl);
@ -463,7 +463,7 @@ void PairSWKokkos<DeviceType>::operator()(TagPairSWComputeFullB<NEIGHFLAG,EVFLAG
delr1[2] = ztmpi - ztmpj; delr1[2] = ztmpi - ztmpj;
const F_FLOAT rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2]; const F_FLOAT rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
if (rsq1 > d_params[jiparam].cutsq) continue; if (rsq1 >= d_params[jiparam].cutsq) continue;
const int j_jnum = d_numneigh[j]; const int j_jnum = d_numneigh[j];
@ -480,7 +480,7 @@ void PairSWKokkos<DeviceType>::operator()(TagPairSWComputeFullB<NEIGHFLAG,EVFLAG
delr2[2] = x(k,2) - ztmpj; delr2[2] = x(k,2) - ztmpj;
const F_FLOAT rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; const F_FLOAT rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
if (rsq2 > d_params[jkparam].cutsq) continue; if (rsq2 >= d_params[jkparam].cutsq) continue;
if (vflag_atom) if (vflag_atom)
threebody(d_params[jiparam],d_params[jkparam],d_params[jikparam], threebody(d_params[jiparam],d_params[jkparam],d_params[jikparam],
@ -905,4 +905,4 @@ void PairSWKokkos<DeviceType>::ev_tally3_atom(EV_FLOAT &ev, const int &i,
template class PairSWKokkos<LMPDeviceType>; template class PairSWKokkos<LMPDeviceType>;
#ifdef KOKKOS_HAVE_CUDA #ifdef KOKKOS_HAVE_CUDA
template class PairSWKokkos<LMPHostType>; template class PairSWKokkos<LMPHostType>;
#endif #endif

37
src/MANYBODY/pair_nb3b_harmonic.cpp
View File

@ -79,9 +79,8 @@ void PairNb3bHarmonic::compute(int eflag, int vflag)
{ {
int i,j,k,ii,jj,kk,inum,jnum,jnumm1; int i,j,k,ii,jj,kk,inum,jnum,jnumm1;
int itype,jtype,ktype,ijparam,ikparam,ijkparam; int itype,jtype,ktype,ijparam,ikparam,ijkparam;
tagint itag,jtag; double xtmp,ytmp,ztmp,evdwl;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl; double rsq1,rsq2;
double rsq,rsq1,rsq2;
double delr1[3],delr2[3],fj[3],fk[3]; double delr1[3],delr2[3],fj[3],fk[3];
int *ilist,*jlist,*numneigh,**firstneigh; int *ilist,*jlist,*numneigh,**firstneigh;
@ -91,7 +90,6 @@ void PairNb3bHarmonic::compute(int eflag, int vflag)
double **x = atom->x; double **x = atom->x;
double **f = atom->f; double **f = atom->f;
tagint *tag = atom->tag;
int *type = atom->type; int *type = atom->type;
inum = list->inum; inum = list->inum;
@ -103,44 +101,13 @@ void PairNb3bHarmonic::compute(int eflag, int vflag)
for (ii = 0; ii < inum; ii++) { for (ii = 0; ii < inum; ii++) {
i = ilist[ii]; i = ilist[ii];
itag = tag[i];
itype = map[type[i]]; itype = map[type[i]];
xtmp = x[i][0]; xtmp = x[i][0];
ytmp = x[i][1]; ytmp = x[i][1];
ztmp = x[i][2]; ztmp = x[i][2];
// two-body interactions, skip half of them
jlist = firstneigh[i]; jlist = firstneigh[i];
jnum = numneigh[i]; jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
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;
}
jtype = map[type[j]];
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
ijparam = elem2param[itype][jtype][jtype];
if (rsq > params[ijparam].cutsq) continue;
}
jnumm1 = jnum - 1; jnumm1 = jnum - 1;
for (jj = 0; jj < jnumm1; jj++) { for (jj = 0; jj < jnumm1; jj++) {

35
src/USER-OMP/pair_nb3b_harmonic_omp.cpp
View File

@ -72,10 +72,9 @@ template <int EVFLAG, int EFLAG>
void PairNb3bHarmonicOMP::eval(int iifrom, int iito, ThrData * const thr) void PairNb3bHarmonicOMP::eval(int iifrom, int iito, ThrData * const thr)
{ {
int i,j,k,ii,jj,kk,jnum,jnumm1; int i,j,k,ii,jj,kk,jnum,jnumm1;
tagint itag,jtag;
int itype,jtype,ktype,ijparam,ikparam,ijkparam; int itype,jtype,ktype,ijparam,ikparam,ijkparam;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl; double xtmp,ytmp,ztmp,evdwl;
double rsq,rsq1,rsq2; double rsq1,rsq2;
double delr1[3],delr2[3],fj[3],fk[3]; double delr1[3],delr2[3],fj[3],fk[3];
int *ilist,*jlist,*numneigh,**firstneigh; int *ilist,*jlist,*numneigh,**firstneigh;
@ -83,7 +82,6 @@ void PairNb3bHarmonicOMP::eval(int iifrom, int iito, ThrData * const thr)
const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0]; const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0]; dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0];
const tagint * _noalias const tag = atom->tag;
const int * _noalias const type = atom->type; const int * _noalias const type = atom->type;
ilist = list->ilist; ilist = list->ilist;
@ -97,7 +95,6 @@ void PairNb3bHarmonicOMP::eval(int iifrom, int iito, ThrData * const thr)
for (ii = iifrom; ii < iito; ++ii) { for (ii = iifrom; ii < iito; ++ii) {
i = ilist[ii]; i = ilist[ii];
itag = tag[i];
itype = map[type[i]]; itype = map[type[i]];
xtmp = x[i].x; xtmp = x[i].x;
ytmp = x[i].y; ytmp = x[i].y;
@ -108,34 +105,6 @@ void PairNb3bHarmonicOMP::eval(int iifrom, int iito, ThrData * const thr)
jlist = firstneigh[i]; jlist = firstneigh[i];
jnum = numneigh[i]; jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
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;
}
jtype = map[type[j]];
delx = xtmp - x[j].x;
dely = ytmp - x[j].y;
delz = ztmp - x[j].z;
rsq = delx*delx + dely*dely + delz*delz;
ijparam = elem2param[itype][jtype][jtype];
if (rsq > params[ijparam].cutsq) continue;
}
jnumm1 = jnum - 1; jnumm1 = jnum - 1;
for (jj = 0; jj < jnumm1; jj++) { for (jj = 0; jj < jnumm1; jj++) {

6
src/USER-OMP/pair_sw_omp.cpp
View File

@ -133,7 +133,7 @@ void PairSWOMP::eval(int iifrom, int iito, ThrData * const thr)
rsq = delx*delx + dely*dely + delz*delz; rsq = delx*delx + dely*dely + delz*delz;
ijparam = elem2param[itype][jtype][jtype]; ijparam = elem2param[itype][jtype][jtype];
if (rsq > params[ijparam].cutsq) continue; if (rsq >= params[ijparam].cutsq) continue;
twobody(&params[ijparam],rsq,fpair,EFLAG,evdwl); twobody(&params[ijparam],rsq,fpair,EFLAG,evdwl);
@ -159,7 +159,7 @@ void PairSWOMP::eval(int iifrom, int iito, ThrData * const thr)
delr1[1] = x[j].y - ytmp; delr1[1] = x[j].y - ytmp;
delr1[2] = x[j].z - ztmp; delr1[2] = x[j].z - ztmp;
rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2]; rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
if (rsq1 > params[ijparam].cutsq) continue; if (rsq1 >= params[ijparam].cutsq) continue;
double fjxtmp,fjytmp,fjztmp; double fjxtmp,fjytmp,fjztmp;
fjxtmp = fjytmp = fjztmp = 0.0; fjxtmp = fjytmp = fjztmp = 0.0;
@ -175,7 +175,7 @@ void PairSWOMP::eval(int iifrom, int iito, ThrData * const thr)
delr2[1] = x[k].y - ytmp; delr2[1] = x[k].y - ytmp;
delr2[2] = x[k].z - ztmp; delr2[2] = x[k].z - ztmp;
rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
if (rsq2 > params[ikparam].cutsq) continue; if (rsq2 >= params[ikparam].cutsq) continue;
threebody(&params[ijparam],&params[ikparam],&params[ijkparam], threebody(&params[ijparam],&params[ikparam],&params[ijkparam],
rsq1,rsq2,delr1,delr2,fj,fk,EFLAG,evdwl); rsq1,rsq2,delr1,delr2,fj,fk,EFLAG,evdwl);