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

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9443 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2013-02-11 15:35:40 +00:00
parent 59b8846dc5
commit a8cb638e76
4 changed files with 93 additions and 88 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

76
src/USER-OMP/pair_lj_cut_tip4p_long_omp.cpp
View File

@ -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);
}
}
}

30
src/USER-OMP/pppm_tip4p_omp.cpp
View File

@ -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;
}
}
}

72
src/USER-OMP/thr_omp.cpp
View File

@ -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;
}
}
}
}

3
src/USER-OMP/thr_omp.h
View File

@ -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);
};