diff --git a/src/USER-OMP/pair_lj_cut_tip4p_long_omp.cpp b/src/USER-OMP/pair_lj_cut_tip4p_long_omp.cpp
index ef2b5120f676e0c1afdd5ea3b0e87fc40fc10e6e..a0dfc0923b6d16ec58a23feb02627c9e103f2c71 100644
--- a/src/USER-OMP/pair_lj_cut_tip4p_long_omp.cpp
+++ b/src/USER-OMP/pair_lj_cut_tip4p_long_omp.cpp
@@ -120,17 +120,16 @@ void PairLJCutTIP4PLongOMP::compute(int eflag, int vflag)
template <int CTABLE, int EVFLAG, int EFLAG, int VFLAG>
void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
{
- int i,j,ii,jj,jnum,itype,jtype,itable;
+ int i,j,ii,jj,jnum,itype,jtype,itable, key;
int n,vlist[6];
int iH1,iH2,jH1,jH2;
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul;
double fraction,table;
- double delxOM,delyOM,delzOM;
double r,rsq,r2inv,r6inv,forcecoul,forcelj,cforce;
double factor_coul,factor_lj;
- double grij,expm2,prefactor,t,erfc,ddotf;
+ double grij,expm2,prefactor,t,erfc;
double v[6],xH1[3],xH2[3];
- double fdx,fdy,fdz,f1x,f1y,f1z,fOx,fOy,fOz,fHx,fHy,fHz;
+ double fdx,fdy,fdz,fOx,fOy,fOz,fHx,fHy,fHz;
const double *x1,*x2;
int *ilist,*jlist,*numneigh,**firstneigh;
@@ -317,6 +316,7 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
// vlist stores 2,4,6 atoms whose forces contribute to virial
n = 0;
+ key = 0;
if (itype != typeO) {
fxtmp += delx * cforce;
@@ -330,33 +330,23 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
v[3] = x[i][0] * dely * cforce;
v[4] = x[i][0] * delz * cforce;
v[5] = x[i][1] * delz * cforce;
- vlist[n++] = i;
}
+ vlist[n++] = i;
} else {
+ key++;
fdx = delx*cforce;
fdy = dely*cforce;
fdz = delz*cforce;
- delxOM = x[i][0] - x1[0];
- delyOM = x[i][1] - x1[1];
- delzOM = x[i][2] - x1[2];
+ fOx = fdx*(1 - alpha);
+ fOy = fdy*(1 - alpha);
+ fOz = fdz*(1 - alpha);
- ddotf = (delxOM * fdx + delyOM * fdy + delzOM * fdz) /
- (qdist*qdist);
-
- f1x = alpha * (fdx - ddotf * delxOM);
- f1y = alpha * (fdy - ddotf * delyOM);
- f1z = alpha * (fdz - ddotf * delzOM);
-
- fOx = fdx - f1x;
- fOy = fdy - f1y;
- fOz = fdz - f1z;
-
- fHx = 0.5 * f1x;
- fHy = 0.5 * f1y;
- fHz = 0.5 * f1z;
+ fHx = 0.5*alpha * fdx;
+ fHy = 0.5*alpha * fdy;
+ fHz = 0.5*alpha * fdz;
fxtmp += fOx;
fytmp += fOy;
@@ -380,11 +370,10 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
v[3] = x[i][0]*fOy + xH1[0]*fHy + xH2[0]*fHy;
v[4] = x[i][0]*fOz + xH1[0]*fHz + xH2[0]*fHz;
v[5] = x[i][1]*fOz + xH1[1]*fHz + xH2[1]*fHz;
-
- vlist[n++] = i;
- vlist[n++] = iH1;
- vlist[n++] = iH2;
}
+ vlist[n++] = i;
+ vlist[n++] = iH1;
+ vlist[n++] = iH2;
}
if (jtype != typeO) {
@@ -399,33 +388,23 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
v[3] -= x[j][0] * dely * cforce;
v[4] -= x[j][0] * delz * cforce;
v[5] -= x[j][1] * delz * cforce;
- vlist[n++] = j;
}
+ vlist[n++] = j;
} else {
+ key += 2;
fdx = -delx*cforce;
fdy = -dely*cforce;
fdz = -delz*cforce;
- delxOM = x[j][0] - x2[0];
- delyOM = x[j][1] - x2[1];
- delzOM = x[j][2] - x2[2];
+ fOx = fdx*(1 - alpha);
+ fOy = fdy*(1 - alpha);
+ fOz = fdz*(1 - alpha);
- ddotf = (delxOM * fdx + delyOM * fdy + delzOM * fdz) /
- (qdist*qdist);
-
- f1x = alpha * (fdx - ddotf * delxOM);
- f1y = alpha * (fdy - ddotf * delyOM);
- f1z = alpha * (fdz - ddotf * delzOM);
-
- fOx = fdx - f1x;
- fOy = fdy - f1y;
- fOz = fdz - f1z;
-
- fHx = 0.5 * f1x;
- fHy = 0.5 * f1y;
- fHz = 0.5 * f1z;
+ fHx = 0.5*alpha * fdx;
+ fHy = 0.5*alpha * fdy;
+ fHz = 0.5*alpha * fdz;
f[j][0] += fOx;
f[j][1] += fOy;
@@ -449,11 +428,10 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
v[3] += x[j][0]*fOy + xH1[0]*fHy + xH2[0]*fHy;
v[4] += x[j][0]*fOz + xH1[0]*fHz + xH2[0]*fHz;
v[5] += x[j][1]*fOz + xH1[1]*fHz + xH2[1]*fHz;
-
- vlist[n++] = j;
- vlist[n++] = jH1;
- vlist[n++] = jH2;
}
+ vlist[n++] = j;
+ vlist[n++] = jH1;
+ vlist[n++] = jH2;
}
if (EFLAG) {
@@ -466,7 +444,7 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
} else ecoul = 0.0;
- if (EVFLAG) ev_tally_list_thr(this,n,vlist,ecoul,v,thr);
+ if (EVFLAG) ev_tally_list_thr(this,key,vlist,v,ecoul,alpha,thr);
}
}
}
diff --git a/src/USER-OMP/pppm_tip4p_omp.cpp b/src/USER-OMP/pppm_tip4p_omp.cpp
index e59275239b71cc865d538c30ce07dcf24fa15f69..e7c6cbe106e61abdcc93f8ff732e6d0b62a4492d 100644
--- a/src/USER-OMP/pppm_tip4p_omp.cpp
+++ b/src/USER-OMP/pppm_tip4p_omp.cpp
@@ -45,6 +45,7 @@ using namespace MathConst;
PPPMTIP4POMP::PPPMTIP4POMP(LAMMPS *lmp, int narg, char **arg) :
PPPMOMP(lmp, narg, arg)
{
+ tip4pflag = 1;
suffix_flag |= Suffix::OMP;
}
@@ -288,7 +289,6 @@ void PPPMTIP4POMP::fieldforce()
FFT_SCALAR ekx,eky,ekz;
int iH1,iH2;
double xM[3], fx,fy,fz;
- double ddotf, rOMx, rOMy, rOMz, f1x, f1y, f1z;
// this if protects against having more threads than local atoms
if (ifrom < nlocal) {
@@ -342,27 +342,17 @@ void PPPMTIP4POMP::fieldforce()
fz = qfactor * ekz;
find_M(i,iH1,iH2,xM);
- rOMx = xM[0] - x[i][0];
- rOMy = xM[1] - x[i][1];
- rOMz = xM[2] - x[i][2];
+ f[i][0] += fx*(1 - alpha);
+ f[i][1] += fy*(1 - alpha);
+ f[i][2] += fz*(1 - alpha);
- ddotf = (rOMx * fx + rOMy * fy + rOMz * fz) / (qdist * qdist);
+ f[iH1][0] += 0.5*alpha*fx;
+ f[iH1][1] += 0.5*alpha*fy;
+ f[iH1][2] += 0.5*alpha*fz;
- f1x = ddotf * rOMx;
- f1y = ddotf * rOMy;
- f1z = ddotf * rOMz;
-
- f[i][0] += fx - alpha * (fx - f1x);
- f[i][1] += fy - alpha * (fy - f1y);
- f[i][2] += fz - alpha * (fz - f1z);
-
- f[iH1][0] += 0.5*alpha*(fx - f1x);
- f[iH1][1] += 0.5*alpha*(fy - f1y);
- f[iH1][2] += 0.5*alpha*(fz - f1z);
-
- f[iH2][0] += 0.5*alpha*(fx - f1x);
- f[iH2][1] += 0.5*alpha*(fy - f1y);
- f[iH2][2] += 0.5*alpha*(fz - f1z);
+ f[iH2][0] += 0.5*alpha*fx;
+ f[iH2][1] += 0.5*alpha*fy;
+ f[iH2][2] += 0.5*alpha*fz;
}
}
}
diff --git a/src/USER-OMP/thr_omp.cpp b/src/USER-OMP/thr_omp.cpp
index 9deb562a867dd9602035f6ffa930f5d2a5ca9f8e..5bfc496389e29fbbec4cced37274a692a4889e3c 100644
--- a/src/USER-OMP/thr_omp.cpp
+++ b/src/USER-OMP/thr_omp.cpp
@@ -626,15 +626,36 @@ void ThrOMP::ev_tally4_thr(Pair * const pair, const int i, const int j,
changes v values by dividing by n
------------------------------------------------------------------------- */
-void ThrOMP::ev_tally_list_thr(Pair * const pair, const int n,
- const int * const list, const double ecoul,
- const double * const v, ThrData * const thr)
+void ThrOMP::ev_tally_list_thr(Pair * const pair, const int key,
+ const int * const list, const double * const v,
+ const double ecoul, const double alpha,
+ ThrData * const thr)
{
+ int i;
if (pair->eflag_either) {
if (pair->eflag_global) thr->eng_coul += ecoul;
if (pair->eflag_atom) {
- double epairatom = ecoul/static_cast<double>(n);
- for (int i = 0; i < n; i++) thr->eatom_pair[list[i]] += epairatom;
+ if (key == 0) {
+ thr->eatom_pair[list[0]] += 0.5*ecoul;
+ thr->eatom_pair[list[1]] += 0.5*ecoul;
+ } else if (key == 1) {
+ thr->eatom_pair[list[0]] += 0.5*ecoul*(1-alpha);
+ thr->eatom_pair[list[1]] += 0.25*ecoul*alpha;
+ thr->eatom_pair[list[2]] += 0.25*ecoul*alpha;
+ thr->eatom_pair[list[3]] += 0.5*ecoul;
+ } else if (key == 2) {
+ thr->eatom_pair[list[0]] += 0.5*ecoul;
+ thr->eatom_pair[list[1]] += 0.5*ecoul*(1-alpha);
+ thr->eatom_pair[list[2]] += 0.25*ecoul*alpha;
+ thr->eatom_pair[list[3]] += 0.25*ecoul*alpha;
+ } else {
+ thr->eatom_pair[list[0]] += 0.5*ecoul*(1-alpha);
+ thr->eatom_pair[list[1]] += 0.25*ecoul*alpha;
+ thr->eatom_pair[list[2]] += 0.25*ecoul*alpha;
+ thr->eatom_pair[list[3]] += 0.5*ecoul*(1-alpha);
+ thr->eatom_pair[list[4]] += 0.25*ecoul*alpha;
+ thr->eatom_pair[list[5]] += 0.25*ecoul*alpha;
+ }
}
}
@@ -643,19 +664,34 @@ void ThrOMP::ev_tally_list_thr(Pair * const pair, const int n,
v_tally(thr->virial_pair,v);
if (pair->vflag_atom) {
- const double s = 1.0/static_cast<double>(n);
- double vtmp[6];
-
- vtmp[0] = s * v[0];
- vtmp[1] = s * v[1];
- vtmp[2] = s * v[2];
- vtmp[3] = s * v[3];
- vtmp[4] = s * v[4];
- vtmp[5] = s * v[5];
-
- for (int i = 0; i < n; i++) {
- const int j = list[i];
- v_tally(thr->vatom_pair[j],vtmp);
+ if (key == 0) {
+ for (i = 0; i <= 5; i++) {
+ thr->vatom_pair[list[0]][i] += 0.5*v[i];
+ thr->vatom_pair[list[1]][i] += 0.5*v[i];
+ }
+ } else if (key == 1) {
+ for (i = 0; i <= 5; i++) {
+ thr->vatom_pair[list[0]][i] += 0.5*v[i]*(1-alpha);
+ thr->vatom_pair[list[1]][i] += 0.25*v[i]*alpha;
+ thr->vatom_pair[list[2]][i] += 0.25*v[i]*alpha;
+ thr->vatom_pair[list[3]][i] += 0.5*v[i];
+ }
+ } else if (key == 2) {
+ for (i = 0; i <= 5; i++) {
+ thr->vatom_pair[list[0]][i] += 0.5*v[i];
+ thr->vatom_pair[list[1]][i] += 0.5*v[i]*(1-alpha);
+ thr->vatom_pair[list[2]][i] += 0.25*v[i]*alpha;
+ thr->vatom_pair[list[3]][i] += 0.25*v[i]*alpha;
+ }
+ } else {
+ for (i = 0; i <= 5; i++) {
+ thr->vatom_pair[list[0]][i] += 0.5*v[i]*(1-alpha);
+ thr->vatom_pair[list[1]][i] += 0.25*v[i]*alpha;
+ thr->vatom_pair[list[2]][i] += 0.25*v[i]*alpha;
+ thr->vatom_pair[list[3]][i] += 0.5*v[i]*(1-alpha);
+ thr->vatom_pair[list[4]][i] += 0.25*v[i]*alpha;
+ thr->vatom_pair[list[5]][i] += 0.25*v[i]*alpha;
+ }
}
}
}
diff --git a/src/USER-OMP/thr_omp.h b/src/USER-OMP/thr_omp.h
index a9f402af97bd670e95b17557e96a4c1df44bc7ee..3c5bbdec3e2fe1440b9729249ea93e88339f15d7 100644
--- a/src/USER-OMP/thr_omp.h
+++ b/src/USER-OMP/thr_omp.h
@@ -162,7 +162,8 @@ class ThrOMP {
const double * const, const double * const, const double * const,
const double * const, const double * const, ThrData * const);
void ev_tally_list_thr(Pair * const, const int, const int * const,
- const double , const double * const , ThrData * const);
+ const double * const, const double, const double,
+ ThrData * const);
};