diff --git a/lib/gpu/lal_sph_heatconduction.cu b/lib/gpu/lal_sph_heatconduction.cu
index e2ba40db0c..8e4ec6ff19 100644
--- a/lib/gpu/lal_sph_heatconduction.cu
+++ b/lib/gpu/lal_sph_heatconduction.cu
@@ -29,23 +29,23 @@ _texture_2d( vel_tex,int4);
 
 #if (SHUFFLE_AVAIL == 0)
 
-#define store_dE(dEacc, ii, inum, tid, t_per_atom, offset, dE)              \
+#define store_dE(dEacc, ii, inum, tid, t_per_atom, offset, i, dE)           \
   if (t_per_atom>1) {                                                       \
     simdsync();                                                             \
     simd_reduce_add1(t_per_atom, red_acc, offset, tid, dEacc);              \
   }                                                                         \
   if (offset==0 && ii<inum) {                                               \
-    dE[ii]=dEacc;                                                           \
+    dE[i]=dEacc;                                                            \
   }
 #else
-#define store_drhoE(dEacc, ii, inum, tid, t_per_atom, offset, dE)           \
+#define store_drhoE(dEacc, ii, inum, tid, t_per_atom, offset, i, dE)        \
   if (t_per_atom>1) {                                                       \
     for (unsigned int s=t_per_atom/2; s>0; s>>=1) {                         \
       dEacc += shfl_down(dEacc, s, t_per_atom);                             \
     }                                                                       \
   }                                                                         \
   if (offset==0 && ii<inum) {                                               \
-    dE[ii]=dEacc;                                                           \
+    dE[i]=dEacc;                                                           \
   }
 #endif
 
@@ -66,7 +66,7 @@ __kernel void k_sph_heatconduction(const __global numtyp4 *restrict x_,
                        const int inum, const int nbor_pitch,
                        const __global numtyp4 *restrict v_,
                        const int dimension, const int t_per_atom) {
-  int tid, ii, offset;
+  int tid, ii, offset, i;
   atom_info(t_per_atom,ii,tid,offset);
 
   int n_stride;
@@ -77,7 +77,7 @@ __kernel void k_sph_heatconduction(const __global numtyp4 *restrict x_,
   acctyp dEacc = (acctyp)0;
 
   if (ii<inum) {
-    int i, numj, nbor, nbor_end;
+    int numj, nbor, nbor_end;
     nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,i,numj,
               n_stride,nbor_end,nbor);
 
@@ -140,7 +140,7 @@ __kernel void k_sph_heatconduction(const __global numtyp4 *restrict x_,
     } // for nbor
   } // if ii
 
-  store_drhoE(dEacc,ii,inum,tid,t_per_atom,offset,dE);
+  store_drhoE(dEacc,ii,inum,tid,t_per_atom,offset,i,dE);
 }
 
 __kernel void k_sph_heatconduction_fast(const __global numtyp4 *restrict x_,
@@ -157,7 +157,7 @@ __kernel void k_sph_heatconduction_fast(const __global numtyp4 *restrict x_,
                             const int inum, const int nbor_pitch,
                             const __global numtyp4 *restrict v_,
                             const int dimension, const int t_per_atom) {
-  int tid, ii, offset;
+  int tid, ii, offset, i;
   atom_info(t_per_atom,ii,tid,offset);
 
   #ifndef ONETYPE
@@ -180,7 +180,7 @@ __kernel void k_sph_heatconduction_fast(const __global numtyp4 *restrict x_,
   acctyp dEacc = (acctyp)0;
 
   if (ii<inum) {
-    int i, numj, nbor, nbor_end;
+    int numj, nbor, nbor_end;
     nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,i,numj,
               n_stride,nbor_end,nbor);
 
@@ -204,7 +204,7 @@ __kernel void k_sph_heatconduction_fast(const __global numtyp4 *restrict x_,
       #endif
 
       numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j];
-      int jtype = jx.w;
+      int jtype=jx.w;
       #ifndef ONETYPE
       int mtype=itype+jx.w;
       const numtyp cutsq_p=coeff[mtype].z;
@@ -252,6 +252,6 @@ __kernel void k_sph_heatconduction_fast(const __global numtyp4 *restrict x_,
     } // for nbor
   } // if ii
 
-  store_drhoE(dEacc,ii,inum,tid,t_per_atom,offset,dE);
+  store_drhoE(dEacc,ii,inum,tid,t_per_atom,offset,i,dE);
 }
 
diff --git a/lib/gpu/lal_sph_lj.cu b/lib/gpu/lal_sph_lj.cu
index f376dfc533..ae8e5aecb5 100644
--- a/lib/gpu/lal_sph_lj.cu
+++ b/lib/gpu/lal_sph_lj.cu
@@ -29,17 +29,18 @@ _texture_2d( vel_tex,int4);
 
 #if (SHUFFLE_AVAIL == 0)
 
-#define store_drhoE(drhoEacc, ii, inum, tid, t_per_atom, offset, drhoE)      \
-  if (t_per_atom>1) {                                                        \
-    simdsync();                                                              \
-    simd_reduce_add2(t_per_atom, red_acc, offset, tid,                       \
-                     drhoEacc.x, drhoEacc.y);                                \
-  }                                                                          \
-  if (offset==0 && ii<inum) {                                                \
-    drhoE[ii]=drhoEacc;                                                      \
+#define store_drhoE(drhoEacc, ii, inum, tid, t_per_atom, offset, i, drhoE)  \
+  if (t_per_atom>1) {                                                       \
+    simdsync();                                                             \
+    simd_reduce_add2(t_per_atom, red_acc, offset, tid,                      \
+                     drhoEacc.x, drhoEacc.y);                               \
+  }                                                                         \
+  if (offset==0 && ii<inum) {                                               \
+    drhoE[i]=drhoEacc.x;                                                    \
+    drhoE[i+inum]=drhoEacc.y;                                               \
   }
 #else
-#define store_drhoE(drhoEacc, ii, inum, tid, t_per_atom, offset, drhoE)     \
+#define store_drhoE(drhoEacc, ii, inum, tid, t_per_atom, offset, i, drhoE)  \
   if (t_per_atom>1) {                                                       \
     for (unsigned int s=t_per_atom/2; s>0; s>>=1) {                         \
       drhoEacc.x += shfl_down(drhoEacc.x, s, t_per_atom);                   \
@@ -47,7 +48,8 @@ _texture_2d( vel_tex,int4);
     }                                                                       \
   }                                                                         \
   if (offset==0 && ii<inum) {                                               \
-    drhoE[ii]=drhoEacc;                                                     \
+    drhoE[i]=drhoEacc.x;                                                    \
+    drhoE[i+inum]=drhoEacc.y;                                               \
   }
 #endif
 
@@ -105,12 +107,12 @@ __kernel void k_sph_lj(const __global numtyp4 *restrict x_,
                        const __global int * dev_packed,
                        __global acctyp3 *restrict ans,
                        __global acctyp *restrict engv,
-                       __global acctyp2 *restrict drhoE,
+                       __global acctyp *restrict drhoE,
                        const int eflag, const int vflag,
                        const int inum, const int nbor_pitch,
                        const __global numtyp4 *restrict v_,
                        const int dimension, const int t_per_atom) {
-  int tid, ii, offset;
+  int tid, ii, offset, i;
   atom_info(t_per_atom,ii,tid,offset);
 
   int n_stride;
@@ -124,10 +126,10 @@ __kernel void k_sph_lj(const __global numtyp4 *restrict x_,
     for (int i=0; i<6; i++) virial[i]=(acctyp)0;
   }
   acctyp2 drhoEacc;
-  drhoEacc.x = drhoEacc.x = (acctyp)0;
+  drhoEacc.x = drhoEacc.y = (acctyp)0;
 
   if (ii<inum) {
-    int i, numj, nbor, nbor_end;
+    int numj, nbor, nbor_end;
     nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,i,numj,
               n_stride,nbor_end,nbor);
 
@@ -246,7 +248,7 @@ __kernel void k_sph_lj(const __global numtyp4 *restrict x_,
   } // if ii
   store_answers(f,energy,virial,ii,inum,tid,t_per_atom,offset,eflag,vflag,
                 ans,engv);
-  store_drhoE(drhoEacc,ii,inum,tid,t_per_atom,offset,drhoE);
+  store_drhoE(drhoEacc,ii,inum,tid,t_per_atom,offset,i,drhoE);
 }
 
 __kernel void k_sph_lj_fast(const __global numtyp4 *restrict x_,
@@ -258,12 +260,12 @@ __kernel void k_sph_lj_fast(const __global numtyp4 *restrict x_,
                             const __global int * dev_packed,
                             __global acctyp3 *restrict ans,
                             __global acctyp *restrict engv,
-                            __global acctyp2 *restrict drhoE,
+                            __global acctyp *restrict drhoE,
                             const int eflag, const int vflag,
                             const int inum, const int nbor_pitch,
                             const __global numtyp4 *restrict v_,
                             const int dimension, const int t_per_atom) {
-  int tid, ii, offset;
+  int tid, ii, offset, i;
   atom_info(t_per_atom,ii,tid,offset);
 
   #ifndef ONETYPE
@@ -289,10 +291,10 @@ __kernel void k_sph_lj_fast(const __global numtyp4 *restrict x_,
     for (int i=0; i<6; i++) virial[i]=(acctyp)0;
   }
   acctyp2 drhoEacc;
-  drhoEacc.x = drhoEacc.x = (acctyp)0;
+  drhoEacc.x = drhoEacc.y = (acctyp)0;
 
   if (ii<inum) {
-    int i, numj, nbor, nbor_end;
+    int numj, nbor, nbor_end;
     nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,i,numj,
               n_stride,nbor_end,nbor);
 
@@ -325,7 +327,7 @@ __kernel void k_sph_lj_fast(const __global numtyp4 *restrict x_,
       #endif
 
       numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j];
-      int jtype = jx.w;
+      int jtype=jx.w;
       #ifndef ONETYPE
       int mtype=itype+jx.w;
       const numtyp cutsq_p=coeff[mtype].z; // cutsq[itype][jtype];
@@ -415,12 +417,11 @@ __kernel void k_sph_lj_fast(const __global numtyp4 *restrict x_,
           virial[4] += delx*delz*force;
           virial[5] += dely*delz*force;
         }
-
       }
     } // for nbor
   } // if ii
 
   store_answers(f,energy,virial,ii,inum,tid,t_per_atom,offset,eflag,vflag, ans,engv);
-  store_drhoE(drhoEacc,ii,inum,tid,t_per_atom,offset,drhoE);
+  store_drhoE(drhoEacc,ii,inum,tid,t_per_atom,offset,i,drhoE);
 }
 
diff --git a/lib/gpu/lal_sph_taitwater.cu b/lib/gpu/lal_sph_taitwater.cu
index 9424d58996..969ce548de 100644
--- a/lib/gpu/lal_sph_taitwater.cu
+++ b/lib/gpu/lal_sph_taitwater.cu
@@ -29,17 +29,18 @@ _texture_2d( vel_tex,int4);
 
 #if (SHUFFLE_AVAIL == 0)
 
-#define store_drhoE(drhoEacc, ii, inum, tid, t_per_atom, offset, drhoE)      \
-  if (t_per_atom>1) {                                                        \
-    simdsync();                                                              \
-    simd_reduce_add2(t_per_atom, red_acc, offset, tid,                       \
-                     drhoEacc.x, drhoEacc.y);                                \
-  }                                                                          \
-  if (offset==0 && ii<inum) {                                                \
-    drhoE[ii]=drhoEacc;                                                      \
+#define store_drhoE(drhoEacc, ii, inum, tid, t_per_atom, offset, i, drhoE)  \
+  if (t_per_atom>1) {                                                       \
+    simdsync();                                                             \
+    simd_reduce_add2(t_per_atom, red_acc, offset, tid,                      \
+                     drhoEacc.x, drhoEacc.y);                               \
+  }                                                                         \
+  if (offset==0 && ii<inum) {                                               \
+    drhoE[i]=drhoEacc.x;                                                    \
+    drhoE[i+inum]=drhoEacc.y;                                               \
   }
 #else
-#define store_drhoE(drhoEacc, ii, inum, tid, t_per_atom, offset, drhoE)     \
+#define store_drhoE(drhoEacc, ii, inum, tid, t_per_atom, offset, i, drhoE)  \
   if (t_per_atom>1) {                                                       \
     for (unsigned int s=t_per_atom/2; s>0; s>>=1) {                         \
       drhoEacc.x += shfl_down(drhoEacc.x, s, t_per_atom);                   \
@@ -47,7 +48,8 @@ _texture_2d( vel_tex,int4);
     }                                                                       \
   }                                                                         \
   if (offset==0 && ii<inum) {                                               \
-    drhoE[ii]=drhoEacc;                                                     \
+    drhoE[i]=drhoEacc.x;                                                    \
+    drhoE[i+inum]=drhoEacc.y;                                               \
   }
 #endif
 
@@ -61,12 +63,12 @@ __kernel void k_sph_taitwater(const __global numtyp4 *restrict x_,
                               const __global int * dev_packed,
                               __global acctyp3 *restrict ans,
                               __global acctyp *restrict engv,
-                              __global acctyp2 *restrict drhoE,
+                              __global acctyp *restrict drhoE,
                               const int eflag, const int vflag,
                               const int inum, const int nbor_pitch,
                               const __global numtyp4 *restrict v_,
                               const int dimension, const int t_per_atom) {
-  int tid, ii, offset;
+  int tid, ii, offset, i;
   atom_info(t_per_atom,ii,tid,offset);
 
   int n_stride;
@@ -80,10 +82,10 @@ __kernel void k_sph_taitwater(const __global numtyp4 *restrict x_,
     for (int i=0; i<6; i++) virial[i]=(acctyp)0;
   }
   acctyp2 drhoEacc;
-  drhoEacc.x = drhoEacc.x = (acctyp)0;
+  drhoEacc.x = drhoEacc.y = (acctyp)0;
 
   if (ii<inum) {
-    int i, numj, nbor, nbor_end;
+    int numj, nbor, nbor_end;
     nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,i,numj,
               n_stride,nbor_end,nbor);
 
@@ -163,9 +165,8 @@ __kernel void k_sph_taitwater(const __global numtyp4 *restrict x_,
         numtyp fvisc = (numtyp)0;
         if (delVdotDelR < (numtyp)0) {
           numtyp mu = h * delVdotDelR / (rsq + (numtyp)0.01 * h * h);
-          fvisc = -coeffx * (soundspeed_itype
-              + soundspeed_jtype) * mu / (rhoi + rhoj);
-        }
+          fvisc = -coeffx * (soundspeed_itype + soundspeed_jtype) * mu / (rhoi + rhoj);
+        } 
 
         // total pair force & thermal energy increment
         numtyp force = -mass_itype * mass_jtype * (fi + fj + fvisc) * wfd;
@@ -176,15 +177,11 @@ __kernel void k_sph_taitwater(const __global numtyp4 *restrict x_,
         f.z+=delz*force;
 
         // and change in density, drho[i]
-        drhoEacc.x += mass_jtype* delVdotDelR * wfd;
+        drhoEacc.x += mass_jtype * delVdotDelR * wfd;
 
         // change in thermal energy, desph[i]
         drhoEacc.y += deltaE;
 
-        if (EVFLAG && eflag) {
-          numtyp e = (numtyp)0;
-          energy+=e;
-        }
         if (EVFLAG && vflag) {
           virial[0] += delx*delx*force;
           virial[1] += dely*dely*force;
@@ -198,7 +195,7 @@ __kernel void k_sph_taitwater(const __global numtyp4 *restrict x_,
   } // if ii
   store_answers(f,energy,virial,ii,inum,tid,t_per_atom,offset,eflag,vflag,
                 ans,engv);
-  store_drhoE(drhoEacc,ii,inum,tid,t_per_atom,offset,drhoE);
+  store_drhoE(drhoEacc,ii,inum,tid,t_per_atom,offset,i,drhoE);
 }
 
 __kernel void k_sph_taitwater_fast(const __global numtyp4 *restrict x_,
@@ -210,12 +207,12 @@ __kernel void k_sph_taitwater_fast(const __global numtyp4 *restrict x_,
                                    const __global int * dev_packed,
                                    __global acctyp3 *restrict ans,
                                    __global acctyp *restrict engv,
-                                   __global acctyp2 *restrict drhoE,
+                                   __global acctyp *restrict drhoE,
                                    const int eflag, const int vflag,
                                    const int inum, const int nbor_pitch,
                                    const __global numtyp4 *restrict v_,
                                    const int dimension, const int t_per_atom) {
-  int tid, ii, offset;
+  int tid, ii, offset, i;
   atom_info(t_per_atom,ii,tid,offset);
 
   #ifndef ONETYPE
@@ -245,10 +242,10 @@ __kernel void k_sph_taitwater_fast(const __global numtyp4 *restrict x_,
     for (int i=0; i<6; i++) virial[i]=(acctyp)0;
   }
   acctyp2 drhoEacc;
-  drhoEacc.x = drhoEacc.x = (acctyp)0;
+  drhoEacc.x = drhoEacc.y = (acctyp)0;
 
   if (ii<inum) {
-    int i, numj, nbor, nbor_end;
+    int numj, nbor, nbor_end;
     nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,i,numj,
               n_stride,nbor_end,nbor);
 
@@ -337,8 +334,7 @@ __kernel void k_sph_taitwater_fast(const __global numtyp4 *restrict x_,
         numtyp fvisc = (numtyp)0;
         if (delVdotDelR < (numtyp)0) {
           numtyp mu = h * delVdotDelR / (rsq + (numtyp)0.01 * h * h);
-          fvisc = -coeffx * (soundspeed_itype
-              + soundspeed_jtype) * mu / (rhoi + rhoj);
+          fvisc = -coeffx * (soundspeed_itype + soundspeed_jtype) * mu / (rhoi + rhoj);
         }
 
         // total pair force & thermal energy increment
@@ -349,16 +345,12 @@ __kernel void k_sph_taitwater_fast(const __global numtyp4 *restrict x_,
         f.y+=dely*force;
         f.z+=delz*force;
 
-        // and change in density
+        // and change in density, drho[i]
         drhoEacc.x += mass_jtype * delVdotDelR * wfd;
 
-        // change in thermal energy
+        // change in thermal energy, desph[i]
         drhoEacc.y += deltaE;
 
-        if (EVFLAG && eflag) {
-          numtyp e = (numtyp)0;
-          energy+=e;
-        }
         if (EVFLAG && vflag) {
           virial[0] += delx*delx*force;
           virial[1] += dely*dely*force;
@@ -372,6 +364,6 @@ __kernel void k_sph_taitwater_fast(const __global numtyp4 *restrict x_,
   } // if ii
 
   store_answers(f,energy,virial,ii,inum,tid,t_per_atom,offset,eflag,vflag, ans,engv);
-  store_drhoE(drhoEacc,ii,inum,tid,t_per_atom,offset,drhoE);
+  store_drhoE(drhoEacc,ii,inum,tid,t_per_atom,offset,i,drhoE);
 }
 
diff --git a/src/GPU/pair_sph_heatconduction_gpu.cpp b/src/GPU/pair_sph_heatconduction_gpu.cpp
index a81de53c91..ba68c5fcc5 100644
--- a/src/GPU/pair_sph_heatconduction_gpu.cpp
+++ b/src/GPU/pair_sph_heatconduction_gpu.cpp
@@ -113,7 +113,7 @@ void PairSPHHeatConductionGPU::compute(int eflag, int vflag)
         neighbor->ago, inum, nall, atom->x, atom->type,
         sublo, subhi, atom->tag, atom->nspecial, atom->special, eflag, vflag,
         eflag_atom, vflag_atom, host_start, &ilist, &numneigh,
-        cpu_time, success, atom->v);
+        cpu_time, success, atom->vest);
   } else {
     inum = list->inum;
     ilist = list->ilist;
@@ -122,7 +122,7 @@ void PairSPHHeatConductionGPU::compute(int eflag, int vflag)
     sph_heatconduction_gpu_compute(neighbor->ago, inum, nall, atom->x, atom->type,
                        ilist, numneigh, firstneigh, eflag, vflag,
                        eflag_atom, vflag_atom, host_start, cpu_time, success,
-                       atom->tag, atom->v);
+                       atom->tag, atom->vest);
   }
   if (!success) error->one(FLERR, "Insufficient memory on accelerator");
 
diff --git a/src/GPU/pair_sph_lj_gpu.cpp b/src/GPU/pair_sph_lj_gpu.cpp
index 46d7b38073..d503a26335 100644
--- a/src/GPU/pair_sph_lj_gpu.cpp
+++ b/src/GPU/pair_sph_lj_gpu.cpp
@@ -114,7 +114,7 @@ void PairSPHLJGPU::compute(int eflag, int vflag)
         neighbor->ago, inum, nall, atom->x, atom->type,
         sublo, subhi, atom->tag, atom->nspecial, atom->special, eflag, vflag,
         eflag_atom, vflag_atom, host_start, &ilist, &numneigh,
-        cpu_time, success, atom->v);
+        cpu_time, success, atom->vest);
   } else {
     inum = list->inum;
     ilist = list->ilist;
@@ -123,7 +123,7 @@ void PairSPHLJGPU::compute(int eflag, int vflag)
     sph_lj_gpu_compute(neighbor->ago, inum, nall, atom->x, atom->type,
                        ilist, numneigh, firstneigh, eflag, vflag,
                        eflag_atom, vflag_atom, host_start, cpu_time, success,
-                       atom->tag, atom->v);
+                       atom->tag, atom->vest);
   }
   if (!success) error->one(FLERR, "Insufficient memory on accelerator");
 
@@ -136,21 +136,21 @@ void PairSPHLJGPU::compute(int eflag, int vflag)
   int nlocal = atom->nlocal;
   if (acc_float) {
     auto drhoE_ptr = (float *)drhoE_pinned;
-    int idx = 0;
-    for (int i = 0; i < nlocal; i++) {
-      drho[i] = drhoE_ptr[idx];
-      desph[i] = drhoE_ptr[idx+1];
-      idx += 2;
-    }
+    for (int i = 0; i < nlocal; i++)
+      drho[i] += drhoE_ptr[i];
+
+    drhoE_ptr += nlocal;
+    for (int i = 0; i < nlocal; i++)
+      desph[i] += drhoE_ptr[i];
 
   } else {
     auto drhoE_ptr = (double *)drhoE_pinned;
-    int idx = 0;
-    for (int i = 0; i < nlocal; i++) {
-      drho[i] = drhoE_ptr[idx];
-      desph[i] = drhoE_ptr[idx+1];
-      idx += 2;
-    }
+    for (int i = 0; i < nlocal; i++)
+      drho[i] += drhoE_ptr[i];
+
+    drhoE_ptr += nlocal;
+    for (int i = 0; i < nlocal; i++)
+      desph[i] += drhoE_ptr[i];
   }
 
   if (atom->molecular != Atom::ATOMIC && neighbor->ago == 0)
diff --git a/src/GPU/pair_sph_taitwater_gpu.cpp b/src/GPU/pair_sph_taitwater_gpu.cpp
index 6f2762c144..23252cea8a 100644
--- a/src/GPU/pair_sph_taitwater_gpu.cpp
+++ b/src/GPU/pair_sph_taitwater_gpu.cpp
@@ -18,6 +18,7 @@
 #include "pair_sph_taitwater_gpu.h"
 
 #include "atom.h"
+#include "comm.h"
 #include "domain.h"
 #include "error.h"
 #include "force.h"
@@ -85,6 +86,25 @@ void PairSPHTaitwaterGPU::compute(int eflag, int vflag)
 {
   ev_init(eflag, vflag);
 
+  // check consistency of pair coefficients
+
+  if (first) {
+    for (int i = 1; i <= atom->ntypes; i++) {
+      for (int j = 1; i <= atom->ntypes; i++) {
+        if (cutsq[i][j] > 1.e-32) {
+          if (!setflag[i][i] || !setflag[j][j]) {
+            if (comm->me == 0) {
+              printf(
+                  "SPH particle types %d and %d interact with cutoff=%g, but not all of their single particle properties are set.\n",
+                  i, j, sqrt(cutsq[i][j]));
+            }
+          }
+        }
+      }
+    }
+    first = 0;
+  }
+
   int nall = atom->nlocal + atom->nghost;
   int inum, host_start;
 
@@ -110,7 +130,7 @@ void PairSPHTaitwaterGPU::compute(int eflag, int vflag)
     firstneigh = sph_taitwater_gpu_compute_n(
         neighbor->ago, inum, nall, atom->x, atom->type, sublo, subhi, atom->tag, atom->nspecial,
         atom->special, eflag, vflag, eflag_atom, vflag_atom, host_start, &ilist, &numneigh,
-        cpu_time, success, atom->v);
+        cpu_time, success, atom->vest);
   } else {
     inum = list->inum;
     ilist = list->ilist;
@@ -118,7 +138,7 @@ void PairSPHTaitwaterGPU::compute(int eflag, int vflag)
     firstneigh = list->firstneigh;
     sph_taitwater_gpu_compute(neighbor->ago, inum, nall, atom->x, atom->type, ilist, numneigh, firstneigh,
                        eflag, vflag, eflag_atom, vflag_atom, host_start, cpu_time, success,
-                       atom->tag, atom->v);
+                       atom->tag, atom->vest);
   }
   if (!success) error->one(FLERR, "Insufficient memory on accelerator");
 
@@ -131,21 +151,21 @@ void PairSPHTaitwaterGPU::compute(int eflag, int vflag)
   int nlocal = atom->nlocal;
   if (acc_float) {
     auto drhoE_ptr = (float *)drhoE_pinned;
-    int idx = 0;
-    for (int i = 0; i < nlocal; i++) {
-      drho[i] = drhoE_ptr[idx];
-      desph[i] = drhoE_ptr[idx+1];
-      idx += 2;
-    }
+    for (int i = 0; i < nlocal; i++)
+      drho[i] += drhoE_ptr[i];
+
+    drhoE_ptr += nlocal;
+    for (int i = 0; i < nlocal; i++)
+      desph[i] += drhoE_ptr[i];
 
   } else {
     auto drhoE_ptr = (double *)drhoE_pinned;
-    int idx = 0;
-    for (int i = 0; i < nlocal; i++) {
-      drho[i] = drhoE_ptr[idx];
-      desph[i] = drhoE_ptr[idx+1];
-      idx += 2;
-    }
+    for (int i = 0; i < nlocal; i++)
+      drho[i] += drhoE_ptr[i];
+
+    drhoE_ptr += nlocal;
+    for (int i = 0; i < nlocal; i++)
+      desph[i] += drhoE_ptr[i];
   }
 
   if (atom->molecular != Atom::ATOMIC && neighbor->ago == 0)