ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

wrong qi*qj cuda code

Browse Source 
correction of cutsq[mtype].z instead of extra[j].x !
This commit is contained in:
Eddy Barraud
2025-03-24 15:35:34 +01:00
parent adaa313990
commit bbd057891a
2 changed files with 12 additions and 12 deletions
Download Patch File Download Diff File Expand all files Collapse all files

22
lib/gpu/lal_dpd_coul_slater_long.cu
View File

@ -287,16 +287,16 @@ __kernel void k_dpd_coul_slater_long(const __global numtyp4 *restrict x_,
// apply Slater electrostatic force if distance below Slater cutoff // apply Slater electrostatic force if distance below Slater cutoff
// and the two species have a slater coeff // and the two species have a slater coeff
// cutsq[mtype].z -> Coulombic squared cutoff // cutsq[mtype].z -> Slater cutoff
if ( cutsq[mtype].z != 0.0 && rsq < cutsq[mtype].z){ // extra[j].x -> q[j] ; particle j charge
if ( rsq < cutsq[mtype].z ){
numtyp r2inv=ucl_recip(rsq); numtyp r2inv=ucl_recip(rsq);
numtyp _erfc; numtyp _erfc;
numtyp grij = g_ewald * r; numtyp grij = g_ewald * r;
numtyp expm2 = ucl_exp(-grij*grij); numtyp expm2 = ucl_exp(-grij*grij);
numtyp t = ucl_recip((numtyp)1.0 + EWALD_P*grij); numtyp t = ucl_recip((numtyp)1.0 + EWALD_P*grij);
_erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; _erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
numtyp prefactor = extra[j].x; numtyp prefactor = qqrd2e * extra[j].x * qtmp / r;
prefactor *= qqrd2e * cutsq[mtype].z * qtmp/r;
numtyp rlamdainv = r * lamdainv; numtyp rlamdainv = r * lamdainv;
numtyp exprlmdainv = ucl_exp((numtyp)-2.0*rlamdainv); numtyp exprlmdainv = ucl_exp((numtyp)-2.0*rlamdainv);
numtyp slater_term = exprlmdainv*((numtyp)1.0 + ((numtyp)2.0*rlamdainv*((numtyp)1.0+rlamdainv))); numtyp slater_term = exprlmdainv*((numtyp)1.0 + ((numtyp)2.0*rlamdainv*((numtyp)1.0+rlamdainv)));
@ -306,9 +306,9 @@ __kernel void k_dpd_coul_slater_long(const __global numtyp4 *restrict x_,
if (EVFLAG && eflag) { if (EVFLAG && eflag) {
numtyp e_slater = ((numtyp)1.0 + rlamdainv)*exprlmdainv; numtyp e_slater = ((numtyp)1.0 + rlamdainv)*exprlmdainv;
numtyp e = prefactor*(_erfc-e_slater); numtyp e_sf = prefactor*(_erfc-e_slater);
if (factor_coul > (numtyp)0) e -= factor_coul*prefactor*((numtyp)1.0 - e_slater); if (factor_coul > (numtyp)0) e_sf -= factor_coul*prefactor*((numtyp)1.0 - e_slater);
e_coul += e; e_coul += e_sf;
} }
} // if cut_coulsq } // if cut_coulsq
@ -471,16 +471,16 @@ __kernel void k_dpd_coul_slater_long_fast(const __global numtyp4 *restrict x_,
// apply Slater electrostatic force if distance below Slater cutoff // apply Slater electrostatic force if distance below Slater cutoff
// and the two species have a slater coeff // and the two species have a slater coeff
// cutsq[mtype].z -> Coulombic squared cutoff // cutsq[mtype].z -> Slater cutoff
if ( cutsq[mtype].z != 0.0 && rsq < cutsq[mtype].z){ // extra[j].x -> q[j] ; particle j charge
if ( rsq < cutsq[mtype].z ){
numtyp r2inv=ucl_recip(rsq); numtyp r2inv=ucl_recip(rsq);
numtyp _erfc; numtyp _erfc;
numtyp grij = g_ewald * r; numtyp grij = g_ewald * r;
numtyp expm2 = ucl_exp(-grij*grij); numtyp expm2 = ucl_exp(-grij*grij);
numtyp t = ucl_recip((numtyp)1.0 + EWALD_P*grij); numtyp t = ucl_recip((numtyp)1.0 + EWALD_P*grij);
_erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; _erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
numtyp prefactor = extra[j].x; numtyp prefactor = qqrd2e * extra[j].x * qtmp / r;
prefactor *= qqrd2e * cutsq[mtype].z * qtmp/r;
numtyp rlamdainv = r * lamdainv; numtyp rlamdainv = r * lamdainv;
numtyp exprlmdainv = ucl_exp((numtyp)-2.0*rlamdainv); numtyp exprlmdainv = ucl_exp((numtyp)-2.0*rlamdainv);
numtyp slater_term = exprlmdainv*((numtyp)1.0 + ((numtyp)2.0*rlamdainv*((numtyp)1.0+rlamdainv))); numtyp slater_term = exprlmdainv*((numtyp)1.0 + ((numtyp)2.0*rlamdainv*((numtyp)1.0+rlamdainv)));

2
lib/gpu/lal_dpd_coul_slater_long.h
View File

@ -65,7 +65,7 @@ class DPDCoulSlaterLong : public BaseDPD<numtyp, acctyp> {
/// coeff.x = a0, coeff.y = gamma, coeff.z = sigma, coeff.w = cut_dpd /// coeff.x = a0, coeff.y = gamma, coeff.z = sigma, coeff.w = cut_dpd
UCL_D_Vec<numtyp4> coeff; UCL_D_Vec<numtyp4> coeff;
/// cutsq.x = cutsq, cutsq.y = cut_dpdsq, cutsq.w = cut_slatersq /// cutsq.x = cutsq, cutsq.y = cut_dpdsq, cutsq.z = cut_slatersq
UCL_D_Vec<numtyp4> cutsq; UCL_D_Vec<numtyp4> cutsq;
/// Special LJ values /// Special LJ values