From a8cb638e765b64c06bfa0636f28185d915024d17 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Mon, 11 Feb 2013 15:35:40 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9443 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/USER-OMP/pair_lj_cut_tip4p_long_omp.cpp | 76 ++++++++------------- src/USER-OMP/pppm_tip4p_omp.cpp | 30 +++----- src/USER-OMP/thr_omp.cpp | 72 ++++++++++++++----- src/USER-OMP/thr_omp.h | 3 +- 4 files changed, 93 insertions(+), 88 deletions(-) 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 ef2b5120f6..a0dfc0923b 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 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 e59275239b..e7c6cbe106 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 9deb562a86..5bfc496389 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(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(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 a9f402af97..3c5bbdec3e 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); };