diff --git a/src/INTEL/pair_dpd_intel.cpp b/src/INTEL/pair_dpd_intel.cpp
index 5ed50e54ba..011979223e 100644
--- a/src/INTEL/pair_dpd_intel.cpp
+++ b/src/INTEL/pair_dpd_intel.cpp
@@ -325,8 +325,11 @@ void PairDPDIntel::eval(const int offload, const int vflag,
           const flt_t rinv = (flt_t)1.0/sqrt(rsq);
 
           if (rinv > icut) {
-            flt_t factor_dpd;
-            if (!ONETYPE) factor_dpd = special_lj[sbindex];
+            flt_t factor_dpd, factor_sqrt;
+            if (!ONETYPE) {
+              factor_dpd = special_lj[sbindex];
+              factor_sqrt = special_lj[sbindex];
+            }
 
             flt_t delvx = vxtmp - v[j].x;
             flt_t delvy = vytmp - v[j].y;
@@ -342,8 +345,11 @@ void PairDPDIntel::eval(const int offload, const int vflag,
               gamma = parami[jtype].gamma;
               sigma = parami[jtype].sigma;
             }
-            flt_t fpair = a0 - iwd * gamma * dot + sigma * randnum * dtinvsqrt;
-            if (!ONETYPE) fpair *= factor_dpd;
+            flt_t fpair = a0 - iwd * gamma * dot;
+            if (!ONETYPE) {
+              fpair *= factor_dpd;
+              fpair += factor_sqrt * sigma * randnum * dtinvsqrt;
+            } else fpair += sigma * randnum * dtinvsqrt;
             fpair *= iwd;
 
             const flt_t fpx = fpair * delx;
@@ -493,8 +499,7 @@ void PairDPDIntel::init_style()
   fix->pair_init_check();
   #ifdef _LMP_INTEL_OFFLOAD
   if (fix->offload_balance() != 0.0)
-    error->all(FLERR,
-          "Offload for dpd/intel is not yet available. Set balance to 0.");
+    error->all(FLERR, "Offload for dpd/intel is not yet available. Set balance to 0.");
   #endif
 
   if (fix->precision() == FixIntel::PREC_MODE_MIXED)
diff --git a/src/KOKKOS/pair_dpd_ext_kokkos.cpp b/src/KOKKOS/pair_dpd_ext_kokkos.cpp
index 42fcff0596..a75497450a 100644
--- a/src/KOKKOS/pair_dpd_ext_kokkos.cpp
+++ b/src/KOKKOS/pair_dpd_ext_kokkos.cpp
@@ -41,8 +41,8 @@ using namespace LAMMPS_NS;
 
 
 template<class DeviceType>
-PairDPDExtKokkos<DeviceType>::PairDPDExtKokkos(class LAMMPS *lmp) :
-  PairDPDExt(lmp) ,
+PairDPDExtKokkos<DeviceType>::PairDPDExtKokkos(class LAMMPS *_lmp) :
+  PairDPDExt(_lmp) ,
 #ifdef DPD_USE_RAN_MARS
   rand_pool(0 /* unused */, lmp)
 #else
@@ -134,6 +134,10 @@ void PairDPDExtKokkos<DeviceType>::compute(int eflagin, int vflagin)
   special_lj[1] = force->special_lj[1];
   special_lj[2] = force->special_lj[2];
   special_lj[3] = force->special_lj[3];
+  special_rf[0] = sqrt(force->special_lj[0]);
+  special_rf[1] = sqrt(force->special_lj[1]);
+  special_rf[2] = sqrt(force->special_lj[2]);
+  special_rf[3] = sqrt(force->special_lj[3]);
 
   nlocal = atom->nlocal;
   dtinvsqrt = 1.0/sqrt(update->dt);
@@ -232,7 +236,8 @@ void PairDPDExtKokkos<DeviceType>::operator() (TagDPDExtKokkos<NEIGHFLAG,EVFLAG>
   int i,j,jj,jnum,itype,jtype;
   double xtmp,ytmp,ztmp,delx,dely,delz,fpairx,fpairy,fpairz,fpair;
   double vxtmp,vytmp,vztmp,delvx,delvy,delvz;
-  double rsq,r,rinv,dot,wd,wdPar,wdPerp,randnum,randnumx,randnumy,randnumz,factor_dpd;
+  double rsq,r,rinv,dot,wd,wdPar,wdPerp,randnum,randnumx,randnumy,randnumz;
+  double prefactor_g,prefactor_s,factor_dpd,factor_sqrt;
   double fx = 0,fy = 0,fz = 0;
   double evdwl = 0;
   i = d_ilist[ii];
@@ -249,6 +254,7 @@ void PairDPDExtKokkos<DeviceType>::operator() (TagDPDExtKokkos<NEIGHFLAG,EVFLAG>
     double P[3][3];
     j = d_neighbors(i,jj);
     factor_dpd = special_lj[sbmask(j)];
+    factor_sqrt = special_rf[sbmask(j)];
     j &= NEIGHMASK;
 
     delx = xtmp - x(j,0);
@@ -291,33 +297,26 @@ void PairDPDExtKokkos<DeviceType>::operator() (TagDPDExtKokkos<NEIGHFLAG,EVFLAG>
 
       // drag force - parallel
       fpair -= params(itype,jtype).gamma*wdPar*wdPar*dot*rinv;
+      fpair *= factor_dpd;
 
       // random force - parallel
-      fpair += params(itype,jtype).sigma*wdPar*randnum*dtinvsqrt;
+      fpair += factor_sqrt*params(itype,jtype).sigma*wdPar*randnum*dtinvsqrt;
 
       fpairx = fpair*rinv*delx;
       fpairy = fpair*rinv*dely;
       fpairz = fpair*rinv*delz;
 
       // drag force - perpendicular
-      fpairx -= params(itype,jtype).gammaT*wdPerp*wdPerp*
-                (P[0][0]*delvx + P[0][1]*delvy + P[0][2]*delvz);
-      fpairy -= params(itype,jtype).gammaT*wdPerp*wdPerp*
-                (P[1][0]*delvx + P[1][1]*delvy + P[1][2]*delvz);
-      fpairz -= params(itype,jtype).gammaT*wdPerp*wdPerp*
-                (P[2][0]*delvx + P[2][1]*delvy + P[2][2]*delvz);
+      prefactor_g = factor_dpd*params(itype,jtype).gammaT*wdPerp*wdPerp;
+      fpairx -= prefactor_g * (P[0][0]*delvx + P[0][1]*delvy + P[0][2]*delvz);
+      fpairy -= prefactor_g * (P[1][0]*delvx + P[1][1]*delvy + P[1][2]*delvz);
+      fpairz -= prefactor_g * (P[2][0]*delvx + P[2][1]*delvy + P[2][2]*delvz);
 
       // random force - perpendicular
-      fpairx += params(itype,jtype).sigmaT*wdPerp*
-                (P[0][0]*randnumx + P[0][1]*randnumy + P[0][2]*randnumz)*dtinvsqrt;
-      fpairy += params(itype,jtype).sigmaT*wdPerp*
-                (P[1][0]*randnumx + P[1][1]*randnumy + P[1][2]*randnumz)*dtinvsqrt;
-      fpairz += params(itype,jtype).sigmaT*wdPerp*
-                (P[2][0]*randnumx + P[2][1]*randnumy + P[2][2]*randnumz)*dtinvsqrt;
-
-      fpairx *= factor_dpd;
-      fpairy *= factor_dpd;
-      fpairz *= factor_dpd;
+      prefactor_s = factor_sqrt*params(itype,jtype).sigmaT*wdPerp;
+      fpairx += prefactor_s * (P[0][0]*randnumx + P[0][1]*randnumy + P[0][2]*randnumz)*dtinvsqrt;
+      fpairy += prefactor_s * (P[1][0]*randnumx + P[1][1]*randnumy + P[1][2]*randnumz)*dtinvsqrt;
+      fpairz += prefactor_s * (P[2][0]*randnumx + P[2][1]*randnumy + P[2][2]*randnumz)*dtinvsqrt;
 
       fx += fpairx;
       fy += fpairy;
diff --git a/src/KOKKOS/pair_dpd_ext_kokkos.h b/src/KOKKOS/pair_dpd_ext_kokkos.h
index 8df218f82c..28965fcb17 100644
--- a/src/KOKKOS/pair_dpd_ext_kokkos.h
+++ b/src/KOKKOS/pair_dpd_ext_kokkos.h
@@ -80,7 +80,7 @@ class PairDPDExtKokkos : public PairDPDExt {
       const F_FLOAT &fx, const F_FLOAT &fy, const F_FLOAT &fz,
       const F_FLOAT &delx, const F_FLOAT &dely, const F_FLOAT &delz) const;
  private:
-  double special_lj[4];
+  double special_lj[4], special_rf[4];
   int eflag,vflag;
   int neighflag,nlocal;
   double dtinvsqrt;
diff --git a/src/KOKKOS/pair_dpd_ext_tstat_kokkos.cpp b/src/KOKKOS/pair_dpd_ext_tstat_kokkos.cpp
index a6c5da37fa..75e9c869f7 100644
--- a/src/KOKKOS/pair_dpd_ext_tstat_kokkos.cpp
+++ b/src/KOKKOS/pair_dpd_ext_tstat_kokkos.cpp
@@ -41,8 +41,8 @@ using namespace LAMMPS_NS;
 
 
 template<class DeviceType>
-PairDPDExtTstatKokkos<DeviceType>::PairDPDExtTstatKokkos(class LAMMPS *lmp) :
-  PairDPDExtTstat(lmp) ,
+PairDPDExtTstatKokkos<DeviceType>::PairDPDExtTstatKokkos(class LAMMPS *_lmp) :
+  PairDPDExtTstat(_lmp) ,
 #ifdef DPD_USE_RAN_MARS
   rand_pool(0 /* unused */, lmp)
 #else
@@ -149,6 +149,10 @@ void PairDPDExtTstatKokkos<DeviceType>::compute(int eflagin, int vflagin)
   special_lj[1] = force->special_lj[1];
   special_lj[2] = force->special_lj[2];
   special_lj[3] = force->special_lj[3];
+  special_rf[0] = sqrt(force->special_lj[0]);
+  special_rf[1] = sqrt(force->special_lj[1]);
+  special_rf[2] = sqrt(force->special_lj[2]);
+  special_rf[3] = sqrt(force->special_lj[3]);
 
   nlocal = atom->nlocal;
   dtinvsqrt = 1.0/sqrt(update->dt);
@@ -233,11 +237,11 @@ void PairDPDExtTstatKokkos<DeviceType>::operator() (TagDPDExtTstatKokkos<NEIGHFL
   auto v_f = ScatterViewHelper<NeedDup_v<NEIGHFLAG,DeviceType>,decltype(dup_f),decltype(ndup_f)>::get(dup_f,ndup_f);
   auto a_f = v_f.template access<AtomicDup_v<NEIGHFLAG,DeviceType>>();
 
-
   int i,j,jj,jnum,itype,jtype;
   double xtmp,ytmp,ztmp,delx,dely,delz,fpairx,fpairy,fpairz,fpair;
   double vxtmp,vytmp,vztmp,delvx,delvy,delvz;
-  double rsq,r,rinv,dot,wd,wdPar,wdPerp,randnum,randnumx,randnumy,randnumz,factor_dpd;
+  double rsq,r,rinv,dot,wd,wdPar,wdPerp,randnum,randnumx,randnumy,randnumz;
+  double prefactor_g,prefactor_s,factor_dpd,factor_sqrt;
   double fx = 0,fy = 0,fz = 0;
 
   i = d_ilist[ii];
@@ -254,6 +258,7 @@ void PairDPDExtTstatKokkos<DeviceType>::operator() (TagDPDExtTstatKokkos<NEIGHFL
     double P[3][3];
     j = d_neighbors(i,jj);
     factor_dpd = special_lj[sbmask(j)];
+    factor_sqrt = special_rf[sbmask(j)];
     j &= NEIGHMASK;
 
     delx = xtmp - x(j,0);
@@ -292,34 +297,27 @@ void PairDPDExtTstatKokkos<DeviceType>::operator() (TagDPDExtTstatKokkos<NEIGHFL
       randnumz = rand_gen.normal();
 
       // drag force - parallel
-      fpair = -params(itype,jtype).gamma*wdPar*wdPar*dot*rinv;
+      fpair = params(itype,jtype).gamma*wdPar*wdPar*dot*rinv;
+      fpair *= factor_dpd;
 
       // random force - parallel
-      fpair += params(itype,jtype).sigma*wdPar*randnum*dtinvsqrt;
+      fpair += factor_sqrt*params(itype,jtype).sigma*wdPar*randnum*dtinvsqrt;
 
       fpairx = fpair*rinv*delx;
       fpairy = fpair*rinv*dely;
       fpairz = fpair*rinv*delz;
 
       // drag force - perpendicular
-      fpairx -= params(itype,jtype).gammaT*wdPerp*wdPerp*
-                (P[0][0]*delvx + P[0][1]*delvy + P[0][2]*delvz);
-      fpairy -= params(itype,jtype).gammaT*wdPerp*wdPerp*
-                (P[1][0]*delvx + P[1][1]*delvy + P[1][2]*delvz);
-      fpairz -= params(itype,jtype).gammaT*wdPerp*wdPerp*
-                (P[2][0]*delvx + P[2][1]*delvy + P[2][2]*delvz);
+      prefactor_g = factor_dpd*params(itype,jtype).gammaT*wdPerp*wdPerp;
+      fpairx -= prefactor_g * (P[0][0]*delvx + P[0][1]*delvy + P[0][2]*delvz);
+      fpairy -= prefactor_g * (P[1][0]*delvx + P[1][1]*delvy + P[1][2]*delvz);
+      fpairz -= prefactor_g * (P[2][0]*delvx + P[2][1]*delvy + P[2][2]*delvz);
 
       // random force - perpendicular
-      fpairx += params(itype,jtype).sigmaT*wdPerp*
-                (P[0][0]*randnumx + P[0][1]*randnumy + P[0][2]*randnumz)*dtinvsqrt;
-      fpairy += params(itype,jtype).sigmaT*wdPerp*
-                (P[1][0]*randnumx + P[1][1]*randnumy + P[1][2]*randnumz)*dtinvsqrt;
-      fpairz += params(itype,jtype).sigmaT*wdPerp*
-                (P[2][0]*randnumx + P[2][1]*randnumy + P[2][2]*randnumz)*dtinvsqrt;
-
-      fpairx *= factor_dpd;
-      fpairy *= factor_dpd;
-      fpairz *= factor_dpd;
+      prefactor_s = factor_sqrt*params(itype,jtype).sigmaT*wdPerp;
+      fpairx += prefactor_s * (P[0][0]*randnumx + P[0][1]*randnumy + P[0][2]*randnumz)*dtinvsqrt;
+      fpairy += prefactor_s * (P[1][0]*randnumx + P[1][1]*randnumy + P[1][2]*randnumz)*dtinvsqrt;
+      fpairz += prefactor_s * (P[2][0]*randnumx + P[2][1]*randnumy + P[2][2]*randnumz)*dtinvsqrt;
 
       fx += fpairx;
       fy += fpairy;
diff --git a/src/KOKKOS/pair_dpd_ext_tstat_kokkos.h b/src/KOKKOS/pair_dpd_ext_tstat_kokkos.h
index 99490894d3..b519204078 100644
--- a/src/KOKKOS/pair_dpd_ext_tstat_kokkos.h
+++ b/src/KOKKOS/pair_dpd_ext_tstat_kokkos.h
@@ -79,7 +79,7 @@ class PairDPDExtTstatKokkos : public PairDPDExtTstat {
       const F_FLOAT &fx,const F_FLOAT &fy, const F_FLOAT &fz,
       const F_FLOAT &delx,const F_FLOAT &dely, const F_FLOAT &delz) const;
  private:
-  double special_lj[4];
+  double special_lj[4], special_rf[4];
   int eflag,vflag;
   int neighflag,nlocal;
   double dtinvsqrt;
diff --git a/src/KOKKOS/pair_dpd_kokkos.cpp b/src/KOKKOS/pair_dpd_kokkos.cpp
index 9ce29528c0..ed9efa572b 100644
--- a/src/KOKKOS/pair_dpd_kokkos.cpp
+++ b/src/KOKKOS/pair_dpd_kokkos.cpp
@@ -41,8 +41,8 @@ using namespace LAMMPS_NS;
 
 
 template<class DeviceType>
-PairDPDKokkos<DeviceType>::PairDPDKokkos(class LAMMPS *lmp) :
-  PairDPD(lmp) ,
+PairDPDKokkos<DeviceType>::PairDPDKokkos(class LAMMPS *_lmp) :
+  PairDPD(_lmp) ,
 #ifdef DPD_USE_RAN_MARS
   rand_pool(0 /* unused */, lmp)
 #else
@@ -134,6 +134,10 @@ void PairDPDKokkos<DeviceType>::compute(int eflagin, int vflagin)
   special_lj[1] = force->special_lj[1];
   special_lj[2] = force->special_lj[2];
   special_lj[3] = force->special_lj[3];
+  special_rf[0] = sqrt(force->special_lj[0]);
+  special_rf[1] = sqrt(force->special_lj[1]);
+  special_rf[2] = sqrt(force->special_lj[2]);
+  special_rf[3] = sqrt(force->special_lj[3]);
 
   nlocal = atom->nlocal;
   dtinvsqrt = 1.0/sqrt(update->dt);
@@ -248,7 +252,7 @@ void PairDPDKokkos<DeviceType>::operator() (TagDPDKokkos<NEIGHFLAG,EVFLAG>, cons
   for (jj = 0; jj < jnum; jj++) {
     j = d_neighbors(i,jj);
     factor_dpd = special_lj[sbmask(j)];
-    factor_sqrt = special_sqrt[sbmask(i)];
+    factor_sqrt = special_rf[sbmask(j)];
     j &= NEIGHMASK;
 
     delx = xtmp - x(j,0);
diff --git a/src/KOKKOS/pair_dpd_kokkos.h b/src/KOKKOS/pair_dpd_kokkos.h
index 4fce75db80..b9f24e3e5f 100644
--- a/src/KOKKOS/pair_dpd_kokkos.h
+++ b/src/KOKKOS/pair_dpd_kokkos.h
@@ -76,10 +76,10 @@ class PairDPDKokkos : public PairDPD {
   template<int NEIGHFLAG>
   KOKKOS_INLINE_FUNCTION
   void ev_tally(EV_FLOAT &ev, const int &i, const int &j,
-      const F_FLOAT &epair, const F_FLOAT &fpair, const F_FLOAT &delx,
-                  const F_FLOAT &dely, const F_FLOAT &delz) const;
+                const F_FLOAT &epair, const F_FLOAT &fpair,
+                const F_FLOAT &delx, const F_FLOAT &dely, const F_FLOAT &delz) const;
  private:
-  double special_lj[4];
+  double special_lj[4], special_rf[4];
   int eflag,vflag;
   int neighflag,nlocal;
   double dtinvsqrt;
diff --git a/src/KOKKOS/pair_dpd_tstat_kokkos.cpp b/src/KOKKOS/pair_dpd_tstat_kokkos.cpp
index 9f835fde56..4c92ac8770 100644
--- a/src/KOKKOS/pair_dpd_tstat_kokkos.cpp
+++ b/src/KOKKOS/pair_dpd_tstat_kokkos.cpp
@@ -41,8 +41,8 @@ using namespace LAMMPS_NS;
 
 
 template<class DeviceType>
-PairDPDTstatKokkos<DeviceType>::PairDPDTstatKokkos(class LAMMPS *lmp) :
-  PairDPDTstat(lmp) ,
+PairDPDTstatKokkos<DeviceType>::PairDPDTstatKokkos(class LAMMPS *_lmp) :
+  PairDPDTstat(_lmp) ,
 #ifdef DPD_USE_RAN_MARS
   rand_pool(0 /* unused */, lmp)
 #else
@@ -128,7 +128,6 @@ void PairDPDTstatKokkos<DeviceType>::compute(int eflagin, int vflagin)
     memory->create(eatom,maxeatom,"pair:eatom");
     memset(&eatom[0], 0, maxeatom * sizeof(double));
   }
-
   if (vflag_atom) {
     memoryKK->destroy_kokkos(k_vatom,vatom);
     memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"pair:vatom");
@@ -149,6 +148,10 @@ void PairDPDTstatKokkos<DeviceType>::compute(int eflagin, int vflagin)
   special_lj[1] = force->special_lj[1];
   special_lj[2] = force->special_lj[2];
   special_lj[3] = force->special_lj[3];
+  special_rf[0] = sqrt(force->special_lj[0]);
+  special_rf[1] = sqrt(force->special_lj[1]);
+  special_rf[2] = sqrt(force->special_lj[2]);
+  special_rf[3] = sqrt(force->special_lj[3]);
 
   nlocal = atom->nlocal;
   dtinvsqrt = 1.0/sqrt(update->dt);
@@ -236,7 +239,7 @@ void PairDPDTstatKokkos<DeviceType>::operator() (TagDPDTstatKokkos<NEIGHFLAG,VFL
   int i,j,jj,jnum,itype,jtype;
   double xtmp,ytmp,ztmp,delx,dely,delz,fpair;
   double vxtmp,vytmp,vztmp,delvx,delvy,delvz;
-  double rsq,r,rinv,dot,wd,randnum,factor_dpd;
+  double rsq,r,rinv,dot,wd,randnum,factor_dpd,factor_sqrt;
   double fx = 0,fy = 0,fz = 0;
 
   i = d_ilist[ii];
@@ -252,6 +255,7 @@ void PairDPDTstatKokkos<DeviceType>::operator() (TagDPDTstatKokkos<NEIGHFLAG,VFL
   for (jj = 0; jj < jnum; jj++) {
     j = d_neighbors(i,jj);
     factor_dpd = special_lj[sbmask(j)];
+    factor_sqrt = special_rf[sbmask(j)];
     j &= NEIGHMASK;
 
     delx = xtmp - x(j,0);
@@ -274,10 +278,11 @@ void PairDPDTstatKokkos<DeviceType>::operator() (TagDPDTstatKokkos<NEIGHFLAG,VFL
 
       // drag force - parallel
       fpair = -params(itype,jtype).gamma*wd*wd*dot*rinv;
+      fpair *= factor_dpd;
 
       // random force - parallel
-      fpair += params(itype,jtype).sigma*wd*randnum*dtinvsqrt;
-      fpair *= factor_dpd*rinv;
+      fpair += factor_sqrt*params(itype,jtype).sigma*wd*randnum*dtinvsqrt;
+      fpair *= rinv;
 
       fx += fpair*delx;
       fy += fpair*dely;
diff --git a/src/KOKKOS/pair_dpd_tstat_kokkos.h b/src/KOKKOS/pair_dpd_tstat_kokkos.h
index 6a655e6892..c277094dba 100644
--- a/src/KOKKOS/pair_dpd_tstat_kokkos.h
+++ b/src/KOKKOS/pair_dpd_tstat_kokkos.h
@@ -79,7 +79,7 @@ class PairDPDTstatKokkos : public PairDPDTstat {
       const F_FLOAT &fpair, const F_FLOAT &delx,
                   const F_FLOAT &dely, const F_FLOAT &delz) const;
  private:
-  double special_lj[4];
+  double special_lj[4], special_rf[4];
   int eflag,vflag;
   int neighflag,nlocal;
   double dtinvsqrt;