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

update conditions, comments, better eflag handle

Browse Source
This commit is contained in:
Eddy Barraud
2024-06-10 11:33:58 +02:00
parent 85c345cf2d
commit e0c2009525
3 changed files with 41 additions and 41 deletions
Download Patch File Download Diff File Expand all files Collapse all files

9
lib/gpu/lal_dpd_charged.cu
View File

@ -425,7 +425,7 @@ __kernel void k_dpd_charged_fast(const __global numtyp4 *restrict x_,
int mtype=itype+jx.w;
// cutsq[mtype].x -> global squared cutoff
/// cutsq.x = cutsq, cutsq.y = cut_dpdsq, cutsq.z = scale, cutsq.w = cut_slatersq
if (rsq<cutsq[mtype].x) {
numtyp r=ucl_sqrt(rsq);
numtyp force_dpd = (numtyp)0.0;
@ -453,6 +453,7 @@ __kernel void k_dpd_charged_fast(const __global numtyp4 *restrict x_,
// conservative force = a0 * wd, or 0 if tstat only
// drag force = -gamma * wd^2 * (delx dot delv) / r
// random force = sigma * wd * rnd * dtinvsqrt;
/// coeff.x = a0, coeff.y = gamma, coeff.z = sigma, coeff.w = cut_dpd
if (!tstat_only) force_dpd = coeff[mtype].x*wd;
force_dpd -= coeff[mtype].y*wd*wd*dot*rinv;
@ -491,9 +492,9 @@ __kernel void k_dpd_charged_fast(const __global numtyp4 *restrict x_,
if (EVFLAG && eflag) {
numtyp e_slater = ((numtyp)1.0 + rlamdainv)*exprlmdainv;
numtyp e = prefactor*(_erfc-e_slater);
if (factor_coul > (numtyp)0) e -= factor_coul*prefactor*((numtyp)1.0 - e_slater);
e_coul += e;
numtyp e_sf = prefactor*(_erfc-e_slater);
if (factor_coul > (numtyp)0) e_sf -= factor_coul*prefactor*((numtyp)1.0 - e_slater);
e_coul += e_sf;
}
} // if cut_coulsq

39
src/DPD-BASIC/pair_dpd_charged.cpp
View File

@ -85,11 +85,9 @@ void PairDPDCharged::compute(int eflag, int vflag)
double rsq,r,rinv,dot,wd,randnum,factor_dpd,factor_sqrt;
int *ilist,*jlist,*numneigh,**firstneigh;
double slater_term;
evdwl = 0.0;
evdwl = ecoul = 0.0;
ev_init(eflag,vflag);
ecoul = 0.0;
double **x = atom->x;
double **v = atom->v;
@ -140,9 +138,10 @@ void PairDPDCharged::compute(int eflag, int vflag)
// forces if below maximum cutoff
if (rsq < cutsq[itype][jtype]) {
r = sqrt(rsq);
if (r < EPSILON) continue; // r can be 0.0 in DPD systems
if (evflag) evdwl = ecoul = 0.0;
// apply DPD force if distance below DPD cutoff
if (rsq < cut_dpdsq[itype][jtype]) {
if (rsq < cut_dpdsq[itype][jtype] && r > EPSILON) {
rinv = 1.0/r;
delvx = vxtmp - v[j][0];
delvy = vytmp - v[j][1];
@ -161,11 +160,18 @@ void PairDPDCharged::compute(int eflag, int vflag)
forcedpd *= factor_dpd;
forcedpd += factor_sqrt*sigma[itype][jtype]*wd*randnum*dtinvsqrt;
forcedpd *= rinv;
if (eflag) {
// eng shifted to 0.0 at cutoff
evdwl = 0.5*a0[itype][jtype]*cut_dpd[itype][jtype] * wd*wd;
evdwl *= factor_dpd;
}
} else forcedpd = 0.0;
// apply Slater electrostatic force if distance below Slater cutoff
// and the two species are charged
if (cut_slater[itype][jtype] != 0.0 && rsq < cut_slatersq[itype][jtype]){
if (rsq < cut_slatersq[itype][jtype]){
r2inv = 1.0/rsq;
grij = g_ewald * r;
expm2 = exp(-grij*grij);
@ -175,7 +181,13 @@ void PairDPDCharged::compute(int eflag, int vflag)
prefactor = qqrd2e * scale[itype][jtype] * qtmp*q[j]/r;
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - slater_term);
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor*(1-slater_term);
forcecoul *= r2inv;
forcecoul *= r2inv;
if (eflag) {
ecoul = prefactor*(erfc - (1 + r/lamda)*exp(-2*r/lamda));
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor*(1.0-(1 + r/lamda)*exp(-2*r/lamda));
}
} else forcecoul = 0.0;
fpair = forcedpd + forcecoul;
@ -189,19 +201,6 @@ void PairDPDCharged::compute(int eflag, int vflag)
f[j][2] -= delz*fpair;
}
// tallies global or per-atom energy and virial only if needed
if (eflag) {
if (rsq < cut_dpdsq[itype][jtype]) {
// eng shifted to 0.0 at cutoff
evdwl = 0.5*a0[itype][jtype]*cut_dpd[itype][jtype] * wd*wd;
evdwl *= factor_dpd;
} else evdwl = 0.0;
if (cut_slater[itype][jtype] != 0.0 && rsq < cut_slatersq[itype][jtype]){
ecoul = prefactor*(erfc - (1 + r/lamda)*exp(-2*r/lamda));
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor*(1.0-(1 + r/lamda)*exp(-2*r/lamda));
} else ecoul = 0.0;
}
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,ecoul,fpair,delx,dely,delz);

34
src/GPU/pair_dpd_charged_gpu.cpp
View File

@ -399,9 +399,9 @@ void PairDPDChargedGPU::cpu_compute(int start, int inum, int eflag, int /* vflag
// forces if below maximum cutoff
if (rsq < cutsq[itype][jtype]) {
r = sqrt(rsq);
if (r < EPSILON) continue; // r can be 0.0 in DPD systems
if (evflag) evdwl = ecoul = 0.0;
// apply DPD force if distance below DPD cutoff
if (rsq < cut_dpdsq[itype][jtype]) {
if (rsq < cut_dpdsq[itype][jtype] && r > EPSILON ) {
rinv = 1.0 / r;
delvx = vxtmp - v[j][0];
delvy = vytmp - v[j][1];
@ -426,11 +426,18 @@ void PairDPDChargedGPU::cpu_compute(int start, int inum, int eflag, int /* vflag
forcedpd *= factor_dpd;
forcedpd += factor_sqrt*sigma[itype][jtype]*wd*randnum*dtinvsqrt;
forcedpd *= rinv;
if (eflag) {
// eng shifted to 0.0 at cutoff
evdwl = 0.5*a0[itype][jtype]*cut_dpd[itype][jtype] * wd*wd;
evdwl *= factor_dpd;
}
} else forcedpd = 0.0;
// apply Slater electrostatic force if distance below Slater cutoff
// and the two species are charged
if (cut_slater[itype][jtype] != 0.0 && rsq < cut_slatersq[itype][jtype]){
if (rsq < cut_slatersq[itype][jtype]){
r2inv = 1.0/rsq;
grij = g_ewald * r;
expm2 = exp(-grij*grij);
@ -440,7 +447,13 @@ void PairDPDChargedGPU::cpu_compute(int start, int inum, int eflag, int /* vflag
prefactor = qqrd2e * scale[itype][jtype] * qtmp*q[j]/r;
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - slater_term);
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor*(1-slater_term);
forcecoul *= r2inv;
forcecoul *= r2inv;
if (eflag) {
ecoul = prefactor*(erfc - (1 + r/lamda)*exp(-2*r/lamda));
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor*(1.0-(1 + r/lamda)*exp(-2*r/lamda));
}
} else forcecoul = 0.0;
fpair = forcedpd + forcecoul;
@ -449,19 +462,6 @@ void PairDPDChargedGPU::cpu_compute(int start, int inum, int eflag, int /* vflag
f[i][1] += dely * fpair;
f[i][2] += delz * fpair;
if (eflag) {
if (rsq < cut_dpdsq[itype][jtype]) {
// eng shifted to 0.0 at cutoff
evdwl = 0.5*a0[itype][jtype]*cut_dpd[itype][jtype] * wd*wd;
evdwl *= factor_dpd;
} else evdwl = 0.0;
if (cut_slater[itype][jtype] != 0.0 && rsq < cut_slatersq[itype][jtype]){
ecoul = prefactor*(erfc - (1 + r/lamda)*exp(-2*r/lamda));
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor*(1.0-(1 + r/lamda)*exp(-2*r/lamda));
} else ecoul = 0.0;
}
if (evflag) ev_tally_full(i, evdwl, ecoul, fpair, delx, dely, delz);
}
}