tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

update and synchronize with implementation of the non-OPENMP version

Browse Source
This commit is contained in:
Axel Kohlmeyer
2025-01-28 21:33:03 -05:00
parent 8b85ee22a3
commit 759a37cc75
2 changed files with 41 additions and 34 deletions
Download Patch File Download Diff File Expand all files Collapse all files

44
src/OPENMP/pair_hbond_dreiding_lj_omp.cpp
View File

@ -120,15 +120,16 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr)
int i,j,k,m,ii,jj,kk,jnum,knum,itype,jtype,ktype,iatom,imol;
tagint tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq,rsq1,rsq2,r1,r2;
double factor_hb,force_angle,force_kernel,evdwl,eng_lj;
double c,s,a,b,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2;
double factor_hb,force_angle,force_kernel,evdwl,eng_lj,ehbond,force_switch;
double c,s,a,b,d,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2;
double fi[3],fj[3],delr1[3],delr2[3];
double r2inv,r10inv;
double switch1,switch2;
int *ilist,*jlist,*numneigh,**firstneigh;
const tagint *klist;
evdwl = 0.0;
evdwl = ehbond = 0.0;
int hbcount = 0;
const auto * _noalias const x = (dbl3_t *) atom->x[0];
auto * _noalias const f = (dbl3_t *) thr->get_f()[0];
@ -151,9 +152,6 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr)
// jj = loop over acceptors
// kk = loop over hydrogens bonded to donor
int hbcount = 0;
double hbeng = 0.0;
for (ii = iifrom; ii < iito; ++ii) {
i = ilist[ii];
itype = type[i];
@ -237,30 +235,34 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr)
r2inv = 1.0/rsq;
r10inv = r2inv*r2inv*r2inv*r2inv*r2inv;
force_kernel = r10inv*(pm.lj1*r2inv - pm.lj2)*r2inv *
powint(c,pm.ap);
force_angle = pm.ap * r10inv*(pm.lj3*r2inv - pm.lj4) *
powint(c,pm.ap-1)*s;
force_kernel = r10inv*(pm.lj1*r2inv - pm.lj2)*r2inv * powint(c,pm.ap);
force_angle = pm.ap * r10inv*(pm.lj3*r2inv - pm.lj4) * powint(c,pm.ap-1)*s;
force_switch = 0.0;
eng_lj = r10inv*(pm.lj3*r2inv - pm.lj4);
if (rsq > pm.cut_innersq) {
switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) *
(pm.cut_outersq + 2.0*rsq - 3.0*pm.cut_innersq) /
pm.denom_vdw;
switch2 = 12.0*rsq * (pm.cut_outersq-rsq) *
(rsq-pm.cut_innersq) / pm.denom_vdw;
force_kernel = force_kernel*switch1 + eng_lj*switch2/rsq;
force_angle *= switch1;
eng_lj *= switch1;
force_kernel *= switch1;
force_angle *= switch1;
force_switch = eng_lj*switch2/rsq;
eng_lj *= switch1;
}
if (EFLAG) {
evdwl = eng_lj * powint(c,pm.ap);
evdwl *= factor_hb;
ehbond += evdwl;
}
a = factor_hb*force_angle/s;
b = factor_hb*force_kernel;
d = factor_hb*force_switch;
a11 = a*c / rsq1;
a12 = -a / (r1*r2);
@ -273,12 +275,12 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr)
vz1 = a11*delr1[2] + a12*delr2[2];
vz2 = a22*delr2[2] + a12*delr1[2];
fi[0] = vx1 + b*delx;
fi[1] = vy1 + b*dely;
fi[2] = vz1 + b*delz;
fj[0] = vx2 - b*delx;
fj[1] = vy2 - b*dely;
fj[2] = vz2 - b*delz;
fi[0] = vx1 + b*delx + d*delx;
fi[1] = vy1 + b*dely + d*dely;
fi[2] = vz1 + b*delz + d*delz;
fj[0] = vx2 - b*delx - d*delx;
fj[1] = vy2 - b*dely - d*dely;
fj[2] = vz2 - b*delz - d*delz;
fxtmp += fi[0];
fytmp += fi[1];
@ -297,7 +299,7 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr)
if (EVFLAG) ev_tally3_thr(this,k,i,j,evdwl,0.0,fi,fj,delr1,delr2,thr);
if (EFLAG) {
hbcount++;
hbeng += evdwl;
ehbond += evdwl;
}
}
}
@ -309,7 +311,7 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr)
}
const int tid = thr->get_tid();
hbcount_thr[tid] = static_cast<double>(hbcount);
hbeng_thr[tid] = hbeng;
hbeng_thr[tid] = ehbond;
}
/* ---------------------------------------------------------------------- */

31
src/OPENMP/pair_hbond_dreiding_morse_omp.cpp
View File

@ -120,14 +120,14 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr)
int i,j,k,m,ii,jj,kk,jnum,knum,itype,jtype,ktype,imol,iatom;
tagint tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq,rsq1,rsq2,r1,r2;
double factor_hb,force_angle,force_kernel,evdwl;
double c,s,a,b,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2;
double factor_hb,force_angle,force_kernel,force_switch,evdwl,ehbond;
double c,s,a,b,d,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2;
double fi[3],fj[3],delr1[3],delr2[3];
double r,dr,dexp,eng_morse,switch1,switch2;
int *ilist,*jlist,*numneigh,**firstneigh;
const tagint *klist;
evdwl = 0.0;
evdwl = ehbond = 0.0;
const auto * _noalias const x = (dbl3_t *) atom->x[0];
auto * _noalias const f = (dbl3_t *) thr->get_f()[0];
@ -151,7 +151,6 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr)
// kk = loop over hydrogens bonded to donor
int hbcount = 0;
double hbeng = 0.0;
for (ii = iifrom; ii < iito; ++ii) {
@ -239,8 +238,10 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr)
dr = r - pm.r0;
dexp = exp(-pm.alpha * dr);
eng_morse = pm.d0 * (dexp*dexp - 2.0*dexp);
force_kernel = pm.morse1*(dexp*dexp - dexp)/r * powint(c,pm.ap);
force_angle = pm.ap * eng_morse * powint(c,pm.ap-1)*s;
force_switch = 0.0;
if (rsq > pm.cut_innersq) {
switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) *
@ -248,18 +249,22 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr)
pm.denom_vdw;
switch2 = 12.0*rsq * (pm.cut_outersq-rsq) *
(rsq-pm.cut_innersq) / pm.denom_vdw;
force_kernel = force_kernel*switch1 + eng_morse*switch2/rsq;
force_kernel *= switch1;
force_angle *= switch1;
force_switch = eng_morse*switch2/rsq;
eng_morse *= switch1;
}
if (EFLAG) {
evdwl = eng_morse * powint(c,pm.ap);
evdwl *= factor_hb;
ehbond += evdwl;
}
a = factor_hb*force_angle/s;
b = factor_hb*force_kernel;
d = factor_hb*force_switch;
a11 = a*c / rsq1;
a12 = -a / (r1*r2);
@ -272,12 +277,12 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr)
vz1 = a11*delr1[2] + a12*delr2[2];
vz2 = a22*delr2[2] + a12*delr1[2];
fi[0] = vx1 + b*delx;
fi[1] = vy1 + b*dely;
fi[2] = vz1 + b*delz;
fj[0] = vx2 - b*delx;
fj[1] = vy2 - b*dely;
fj[2] = vz2 - b*delz;
fi[0] = vx1 + b*delx + d*delx;
fi[1] = vy1 + b*dely + d*dely;
fi[2] = vz1 + b*delz + d*delz;
fj[0] = vx2 - b*delx - d*delx;
fj[1] = vy2 - b*dely - d*dely;
fj[2] = vz2 - b*delz - d*delz;
fxtmp += fi[0];
fytmp += fi[1];
@ -296,7 +301,7 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr)
if (EVFLAG) ev_tally3_thr(this,k,i,j,evdwl,0.0,fi,fj,delr1,delr2,thr);
if (EFLAG) {
hbcount++;
hbeng += evdwl;
ehbond += evdwl;
}
}
}
@ -308,7 +313,7 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr)
}
const int tid = thr->get_tid();
hbcount_thr[tid] = static_cast<double>(hbcount);
hbeng_thr[tid] = hbeng;
hbeng_thr[tid] = ehbond;
}
/* ---------------------------------------------------------------------- */