diff --git a/src/KSPACE/pair_lj_cut_tip4p_long.cpp b/src/KSPACE/pair_lj_cut_tip4p_long.cpp
index 7c4a239833..7d6d4a8fa2 100644
--- a/src/KSPACE/pair_lj_cut_tip4p_long.cpp
+++ b/src/KSPACE/pair_lj_cut_tip4p_long.cpp
@@ -33,6 +33,7 @@
 
 #include <cmath>
 #include <cstring>
+#include <limits>
 
 using namespace LAMMPS_NS;
 using namespace EwaldConst;
@@ -56,6 +57,9 @@ PairLJCutTIP4PLong::PairLJCutTIP4PLong(LAMMPS *lmp) :
   nmax = 0;
   hneigh = nullptr;
   newsite = nullptr;
+  neigh_stamp = 1;
+  site_stamp = 1;
+  tip4p_cache = nullptr;
 
   // TIP4P cannot compute virial as F dot r
   // due to finding bonded H atoms which are not near O atom
@@ -69,6 +73,7 @@ PairLJCutTIP4PLong::~PairLJCutTIP4PLong()
 {
   memory->destroy(hneigh);
   memory->destroy(newsite);
+  memory->destroy(tip4p_cache);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -77,23 +82,15 @@ void PairLJCutTIP4PLong::compute(int eflag, int vflag)
 {
   ev_init(eflag,vflag);
 
-  // reallocate hneigh & newsite if necessary
-  // initialize hneigh[0] to -1 on steps when reneighboring occurred
-  // initialize hneigh[2] to 0 every step
+  grow_tip4p_cache();
 
-  const int nlocal = atom->nlocal;
-  const int nall = nlocal + atom->nghost;
-
-  if (atom->nmax > nmax) {
-    nmax = atom->nmax;
-    memory->destroy(hneigh);
-    memory->create(hneigh,nmax,3,"pair:hneigh");
-    memory->destroy(newsite);
-    memory->create(newsite,nmax,3,"pair:newsite");
+  if (site_stamp == std::numeric_limits<int>::max()) reset_tip4p_cache_stamps();
+  ++site_stamp;
+
+  if (neighbor->ago == 0) {
+    if (neigh_stamp == std::numeric_limits<int>::max()) reset_tip4p_cache_stamps();
+    ++neigh_stamp;
   }
-  if (neighbor->ago == 0)
-    for (int i = 0; i < nall; i++) hneigh[i][0] = -1;
-  for (int i = 0; i < nall; i++) hneigh[i][2] = 0;
 
   if (!ncoultablebits) {
     if (evflag) {
@@ -134,6 +131,7 @@ void PairLJCutTIP4PLong::eval()
   const int inum = list->inum;
   const int nlocal = atom->nlocal;
   const double cut_coulsqplus = (cut_coul+2.0*qdist) * (cut_coul+2.0*qdist);
+  const bool pure_water = (atom->ntypes == 2);
 
   for (int ii = 0; ii < inum; ii++) {
     const int i = ilist[ii];
@@ -154,21 +152,32 @@ void PairLJCutTIP4PLong::eval()
       const double cut_coul_plus_i_site =
         cut_coul + sqrt(site_delx*site_delx + site_dely*site_dely + site_delz*site_delz);
       const double cut_coulsq_i_h = cut_coul_plus_i_site * cut_coul_plus_i_site;
-      eval_i<CTABLE,EVFLAG,EFLAG,VFLAG,1,0>(i,itype,qtmp,xtmp,ytmp,ztmp,iH1,iH2,x1,jlist,jnum,
-                                            nlocal,cut_coulsq_i_h,cut_coulsqplus);
+      if (pure_water)
+        eval_i<CTABLE,EVFLAG,EFLAG,VFLAG,1,0,1>(i,itype,qtmp,xtmp,ytmp,ztmp,iH1,iH2,x1,jlist,
+                                                jnum,nlocal,cut_coulsq_i_h,
+                                                cut_coulsqplus);
+      else
+        eval_i<CTABLE,EVFLAG,EFLAG,VFLAG,1,0,0>(i,itype,qtmp,xtmp,ytmp,ztmp,iH1,iH2,x1,jlist,
+                                                jnum,nlocal,cut_coulsq_i_h,
+                                                cut_coulsqplus);
     } else if (itype == typeH) {
-      eval_i<CTABLE,EVFLAG,EFLAG,VFLAG,0,1>(i,itype,qtmp,xtmp,ytmp,ztmp,-1,-1,x_d[i],jlist,jnum,
-                                            nlocal,cut_coulsq,cut_coulsqplus);
+      if (pure_water)
+        eval_i<CTABLE,EVFLAG,EFLAG,VFLAG,0,1,1>(i,itype,qtmp,xtmp,ytmp,ztmp,-1,-1,x_d[i],jlist,
+                                                jnum,nlocal,cut_coulsq,cut_coulsqplus);
+      else
+        eval_i<CTABLE,EVFLAG,EFLAG,VFLAG,0,1,0>(i,itype,qtmp,xtmp,ytmp,ztmp,-1,-1,x_d[i],jlist,
+                                                jnum,nlocal,cut_coulsq,cut_coulsqplus);
     } else {
-      eval_i<CTABLE,EVFLAG,EFLAG,VFLAG,0,0>(i,itype,qtmp,xtmp,ytmp,ztmp,-1,-1,x_d[i],jlist,jnum,
-                                            nlocal,cut_coulsq,cut_coulsqplus);
+      eval_i<CTABLE,EVFLAG,EFLAG,VFLAG,0,0,0>(i,itype,qtmp,xtmp,ytmp,ztmp,-1,-1,x_d[i],jlist,
+                                              jnum,nlocal,cut_coulsq,cut_coulsqplus);
     }
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
-template <int CTABLE, int EVFLAG, int EFLAG, int VFLAG, int I_WATER, int I_HYDROGEN>
+template <int CTABLE, int EVFLAG, int EFLAG, int VFLAG, int I_WATER, int I_HYDROGEN,
+          int PURE_WATER>
 void PairLJCutTIP4PLong::eval_i(int i, int itype, double qtmp, double xtmp, double ytmp,
                                 double ztmp, int iH1, int iH2, double *x1, int *jlist,
                                 int jnum, int nlocal, double cut_coulsq_i_h,
@@ -201,6 +210,12 @@ void PairLJCutTIP4PLong::eval_i(int i, int itype, double qtmp, double xtmp, doub
   double *lj3_i = nullptr;
   double *lj4_i = nullptr;
   double *offset_i = nullptr;
+  double cut_ljsq_oo = 0.0;
+  double lj1_oo = 0.0;
+  double lj2_oo = 0.0;
+  double lj3_oo = 0.0;
+  double lj4_oo = 0.0;
+  double offset_oo = 0.0;
 
   if constexpr (!I_HYDROGEN) {
     cut_ljsq_i = cut_ljsq[itype];
@@ -209,12 +224,22 @@ void PairLJCutTIP4PLong::eval_i(int i, int itype, double qtmp, double xtmp, doub
     lj3_i = lj3[itype];
     lj4_i = lj4[itype];
     offset_i = offset[itype];
+
+    if constexpr (PURE_WATER) {
+      cut_ljsq_oo = cut_ljsq_i[typeO];
+      lj1_oo = lj1_i[typeO];
+      lj2_oo = lj2_i[typeO];
+      lj3_oo = lj3_i[typeO];
+      lj4_oo = lj4_i[typeO];
+      offset_oo = offset_i[typeO];
+    }
   }
 
   for (jj = 0; jj < jnum; jj++) {
     j = jlist[jj];
     const int sb = sbmask(j);
-    factor_coul = special_coul[sb];
+    const bool special = (sb != 0);
+    factor_coul = special ? special_coul[sb] : 1.0;
     j &= NEIGHMASK;
 
     delx = xtmp - x[j].x;
@@ -223,11 +248,33 @@ void PairLJCutTIP4PLong::eval_i(int i, int itype, double qtmp, double xtmp, doub
     rsq = delx*delx + dely*dely + delz*delz;
     jtype = type[j];
     const bool j_is_oxygen = (jtype == typeO);
-    const bool j_is_hydrogen = (jtype == typeH);
 
     if constexpr (!I_HYDROGEN) {
-      if (!j_is_hydrogen && rsq < cut_ljsq_i[jtype]) {
-        const double factor_lj = special_lj[sb];
+      if constexpr (PURE_WATER) {
+        if (j_is_oxygen && rsq < cut_ljsq_oo) {
+          const double factor_lj = special ? special_lj[sb] : 1.0;
+          r2inv = 1.0/rsq;
+          r6inv = r2inv*r2inv*r2inv;
+          forcelj = r6inv * (lj1_oo*r6inv - lj2_oo);
+          forcelj *= factor_lj * r2inv;
+
+          f[i].x += delx*forcelj;
+          f[i].y += dely*forcelj;
+          f[i].z += delz*forcelj;
+          f[j].x -= delx*forcelj;
+          f[j].y -= dely*forcelj;
+          f[j].z -= delz*forcelj;
+
+          if constexpr (EFLAG) {
+            evdwl = r6inv*(lj3_oo*r6inv-lj4_oo) - offset_oo;
+            evdwl *= factor_lj;
+          } else evdwl = 0.0;
+
+          if constexpr (EVFLAG)
+            ev_tally(i,j,nlocal,1,evdwl,0.0,forcelj,delx,dely,delz);
+        }
+      } else if (jtype != typeH && rsq < cut_ljsq_i[jtype]) {
+        const double factor_lj = special ? special_lj[sb] : 1.0;
         r2inv = 1.0/rsq;
         r6inv = r2inv*r2inv*r2inv;
         forcelj = r6inv * (lj1_i[jtype]*r6inv - lj2_i[jtype]);
@@ -250,8 +297,7 @@ void PairLJCutTIP4PLong::eval_i(int i, int itype, double qtmp, double xtmp, doub
       }
     }
 
-    double cut_coulsq_precheck = cut_coulsqplus;
-    if (!j_is_oxygen) cut_coulsq_precheck = cut_coulsq_i_h;
+    const double cut_coulsq_precheck = j_is_oxygen ? cut_coulsqplus : cut_coulsq_i_h;
 
     if (rsq < cut_coulsq_precheck) {
       if constexpr (I_WATER) {
@@ -439,10 +485,41 @@ void PairLJCutTIP4PLong::eval_i(int i, int itype, double qtmp, double xtmp, doub
 
 /* ---------------------------------------------------------------------- */
 
+void PairLJCutTIP4PLong::grow_tip4p_cache()
+{
+  if (atom->nmax <= nmax) return;
+
+  nmax = atom->nmax;
+  memory->destroy(tip4p_cache);
+  memory->create(tip4p_cache,nmax,"pair:tip4p_cache");
+
+  for (int i = 0; i < nmax; i++) {
+    tip4p_cache[i].neigh_stamp = 0;
+    tip4p_cache[i].site_stamp = 0;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairLJCutTIP4PLong::reset_tip4p_cache_stamps()
+{
+  for (int i = 0; i < nmax; i++) {
+    tip4p_cache[i].neigh_stamp = 0;
+    tip4p_cache[i].site_stamp = 0;
+  }
+
+  neigh_stamp = 1;
+  site_stamp = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
 double *PairLJCutTIP4PLong::tip4p_site(int i, int &iH1, int &iH2, double **x, tagint *tag,
                                        int *type)
 {
-  if (hneigh[i][0] < 0) {
+  Tip4pSiteCache &cache = tip4p_cache[i];
+
+  if (cache.neigh_stamp != neigh_stamp) {
     iH1 = atom->map(tag[i] + 1);
     iH2 = atom->map(tag[i] + 2);
     if (iH1 == -1 || iH2 == -1)
@@ -452,21 +529,22 @@ double *PairLJCutTIP4PLong::tip4p_site(int i, int &iH1, int &iH2, double **x, ta
 
     iH1 = domain->closest_image(i,iH1);
     iH2 = domain->closest_image(i,iH2);
-    compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);
-    hneigh[i][0] = iH1;
-    hneigh[i][1] = iH2;
-    hneigh[i][2] = 1;
+    cache.iH1 = iH1;
+    cache.iH2 = iH2;
+    cache.neigh_stamp = neigh_stamp;
+    cache.site_stamp = 0;
 
   } else {
-    iH1 = hneigh[i][0];
-    iH2 = hneigh[i][1];
-    if (hneigh[i][2] == 0) {
-      hneigh[i][2] = 1;
-      compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);
-    }
+    iH1 = cache.iH1;
+    iH2 = cache.iH2;
+  }
+
+  if (cache.site_stamp != site_stamp) {
+    compute_newsite(x[i],x[iH1],x[iH2],cache.xM);
+    cache.site_stamp = site_stamp;
   }
 
-  return newsite[i];
+  return cache.xM;
 }
 
 /* ----------------------------------------------------------------------
@@ -679,6 +757,6 @@ double PairLJCutTIP4PLong::memory_usage()
 {
   double bytes = (double)maxeatom * sizeof(double);
   bytes += (double)maxvatom*6 * sizeof(double);
-  bytes += (double)2 * nmax * sizeof(double);
+  bytes += (double)nmax * sizeof(Tip4pSiteCache);
   return bytes;
 }
diff --git a/src/KSPACE/pair_lj_cut_tip4p_long.h b/src/KSPACE/pair_lj_cut_tip4p_long.h
index b97b0319cb..233fbb0dcc 100644
--- a/src/KSPACE/pair_lj_cut_tip4p_long.h
+++ b/src/KSPACE/pair_lj_cut_tip4p_long.h
@@ -39,20 +39,32 @@ class PairLJCutTIP4PLong : public PairLJCutCoulLong {
   double memory_usage() override;
 
  protected:
+  struct Tip4pSiteCache {
+    int iH1;
+    int iH2;
+    int neigh_stamp;
+    int site_stamp;
+    double xM[3];
+  };
+
   std::string typeH_str, typeO_str, typeA_str, typeB_str;
   int typeH, typeO;    // atom types of TIP4P water H and O atoms
   int typeA, typeB;    // angle and bond types of TIP4P water
   double alpha;        // geometric constraint parameter for TIP4P
 
   int nmax;            // info on off-oxygen charge sites
-  int **hneigh;        // 0,1 = indices of 2 H associated with O
-                       // 2 = 0 if site loc not yet computed, 1 if yes
-  double **newsite;    // locations of charge sites
+  int **hneigh;        // legacy OPT cache storage
+  double **newsite;    // legacy OPT cache storage
+  int neigh_stamp;
+  int site_stamp;
+  Tip4pSiteCache *tip4p_cache;
 
   template <int, int, int, int> void eval();
-  template <int, int, int, int, int, int>
+  template <int, int, int, int, int, int, int>
   void eval_i(int, int, double, double, double, double, int, int, double *, int *, int, int,
               double, double);
+  void grow_tip4p_cache();
+  void reset_tip4p_cache_stamps();
   double *tip4p_site(int, int &, int &, double **, tagint *, int *);
   void compute_newsite(double *, double *, double *, double *);
 };
diff --git a/src/KSPACE/pppm_tip4p.cpp b/src/KSPACE/pppm_tip4p.cpp
index ad04d3a89c..4ddee644a8 100644
--- a/src/KSPACE/pppm_tip4p.cpp
+++ b/src/KSPACE/pppm_tip4p.cpp
@@ -27,6 +27,7 @@
 
 #include <cmath>
 #include <cstring>
+#include <limits>
 
 using namespace LAMMPS_NS;
 using namespace MathConst;
@@ -50,20 +51,21 @@ PPPMTIP4P::PPPMTIP4P(LAMMPS *lmp) : PPPM(lmp)
   tip4pflag = 1;
 
   tip4p_cache_nmax = 0;
-  tip4p_h1 = nullptr;
-  tip4p_h2 = nullptr;
-  tip4p_cache_valid = nullptr;
-  tip4p_xm = nullptr;
+  tip4p_cache_stamp = 1;
+  tip4p_oxygen_count = 0;
+  tip4p_nonoxygen_count = 0;
+  tip4p_oxygen = nullptr;
+  tip4p_nonoxygen = nullptr;
+  tip4p_cache = nullptr;
 }
 
 /* ---------------------------------------------------------------------- */
 
 PPPMTIP4P::~PPPMTIP4P()
 {
-  memory->destroy(tip4p_h1);
-  memory->destroy(tip4p_h2);
-  memory->destroy(tip4p_cache_valid);
-  memory->destroy(tip4p_xm);
+  memory->destroy(tip4p_oxygen);
+  memory->destroy(tip4p_nonoxygen);
+  memory->destroy(tip4p_cache);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -86,14 +88,22 @@ void PPPMTIP4P::grow_tip4p_cache()
 
   tip4p_cache_nmax = atom->nmax;
 
-  memory->destroy(tip4p_h1);
-  memory->create(tip4p_h1,tip4p_cache_nmax,"pppm/tip4p:h1");
-  memory->destroy(tip4p_h2);
-  memory->create(tip4p_h2,tip4p_cache_nmax,"pppm/tip4p:h2");
-  memory->destroy(tip4p_cache_valid);
-  memory->create(tip4p_cache_valid,tip4p_cache_nmax,"pppm/tip4p:valid");
-  memory->destroy(tip4p_xm);
-  memory->create(tip4p_xm,tip4p_cache_nmax,3,"pppm/tip4p:xm");
+  memory->destroy(tip4p_oxygen);
+  memory->create(tip4p_oxygen,tip4p_cache_nmax,"pppm/tip4p:oxygen");
+  memory->destroy(tip4p_nonoxygen);
+  memory->create(tip4p_nonoxygen,tip4p_cache_nmax,"pppm/tip4p:nonoxygen");
+  memory->destroy(tip4p_cache);
+  memory->create(tip4p_cache,tip4p_cache_nmax,"pppm/tip4p:cache");
+
+  for (int i = 0; i < tip4p_cache_nmax; i++) tip4p_cache[i].stamp = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PPPMTIP4P::reset_tip4p_cache_stamps()
+{
+  for (int i = 0; i < tip4p_cache_nmax; i++) tip4p_cache[i].stamp = 0;
+  tip4p_cache_stamp = 1;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -102,14 +112,15 @@ void PPPMTIP4P::find_M_cached(int i, int &iH1, int &iH2, double *&xM)
 {
   if (atom->nmax > tip4p_cache_nmax) grow_tip4p_cache();
 
-  if (!tip4p_cache_valid[i]) {
-    find_M(i,tip4p_h1[i],tip4p_h2[i],tip4p_xm[i]);
-    tip4p_cache_valid[i] = 1;
+  Tip4pCache &cache = tip4p_cache[i];
+  if (cache.stamp != tip4p_cache_stamp) {
+    find_M(i,cache.iH1,cache.iH2,cache.xM);
+    cache.stamp = tip4p_cache_stamp;
   }
 
-  iH1 = tip4p_h1[i];
-  iH2 = tip4p_h2[i];
-  xM = tip4p_xm[i];
+  iH1 = cache.iH1;
+  iH2 = cache.iH2;
+  xM = cache.xM;
 }
 
 /* ----------------------------------------------------------------------
@@ -131,10 +142,9 @@ void PPPMTIP4P::particle_map()
   const double boxloy = boxlo[1];
   const double boxloz = boxlo[2];
 
-  if (nlocal > 0) {
-    grow_tip4p_cache();
-    std::memset(tip4p_cache_valid,0,nlocal*sizeof(int));
-  }
+  if (nlocal > 0) grow_tip4p_cache();
+  if (tip4p_cache_stamp == std::numeric_limits<int>::max()) reset_tip4p_cache_stamps();
+  ++tip4p_cache_stamp;
 
   if (!std::isfinite(boxlo[0]) || !std::isfinite(boxlo[1]) || !std::isfinite(boxlo[2]))
     error->one(FLERR,"Non-numeric box dimensions - simulation unstable" + utils::errorurl(6));
@@ -186,9 +196,8 @@ void PPPMTIP4P::particle_map()
 
 void PPPMTIP4P::make_rho()
 {
-  int l,m,n,nx,ny,nz,iH1,iH2;
+  int i,l,m,n,nx,ny,nz;
   FFT_SCALAR dx,dy,dz,x0,y0,z0;
-  double *xM;
 
   // clear 3d density array
 
@@ -200,11 +209,9 @@ void PPPMTIP4P::make_rho()
   // (dx,dy,dz) = distance to "lower left" grid pt
   // (mx,my,mz) = global coords of moving stencil pt
 
-  int *type = atom->type;
   double *q = atom->q;
   const auto * _noalias const x = (dbl3_t *) atom->x[0];
   const auto * _noalias const p2g = (int3_t *) part2grid[0];
-  int nlocal = atom->nlocal;
   const int nx_out = nxhi_out - nxlo_out + 1;
   const int ny_out = nyhi_out - nylo_out + 1;
   const int nyz_out = nx_out * ny_out;
@@ -212,10 +219,13 @@ void PPPMTIP4P::make_rho()
   const double boxloy = boxlo[1];
   const double boxloz = boxlo[2];
 
-  for (int i = 0; i < nlocal; i++) {
+  int *type = atom->type;
+  const int nlocal = atom->nlocal;
+
+  for (i = 0; i < nlocal; i++) {
     double xix, xiy, xiz;
     if (type[i] == typeO) {
-      find_M_cached(i,iH1,iH2,xM);
+      const double *xM = tip4p_cache[i].xM;
       xix = xM[0];
       xiy = xM[1];
       xiz = xM[2];
@@ -256,10 +266,9 @@ void PPPMTIP4P::make_rho()
 
 void PPPMTIP4P::fieldforce_ik()
 {
-  int i,l,m,n,nx,ny,nz,iH1,iH2;
+  int i,l,m,n,nx,ny,nz;
   FFT_SCALAR dx,dy,dz,x0,y0,z0;
   FFT_SCALAR ekx,eky,ekz;
-  double *xM;
   double fx,fy,fz;
 
   // loop over my charges, interpolate electric field from nearby grid points
@@ -278,8 +287,6 @@ void PPPMTIP4P::fieldforce_ik()
   FFT_SCALAR *vdy = &vdy_brick[nzlo_out][nylo_out][nxlo_out];
   FFT_SCALAR *vdz = &vdz_brick[nzlo_out][nylo_out][nxlo_out];
 
-  int *type = atom->type;
-  int nlocal = atom->nlocal;
   const int nx_out = nxhi_out - nxlo_out + 1;
   const int ny_out = nyhi_out - nylo_out + 1;
   const int nyz_out = nx_out * ny_out;
@@ -287,13 +294,17 @@ void PPPMTIP4P::fieldforce_ik()
   const double boxloy = boxlo[1];
   const double boxloz = boxlo[2];
 
+  int *type = atom->type;
+  const int nlocal = atom->nlocal;
+
   for (i = 0; i < nlocal; i++) {
+    const Tip4pCache *cache = nullptr;
     double xix, xiy, xiz;
     if (type[i] == typeO) {
-      find_M_cached(i,iH1,iH2,xM);
-      xix = xM[0];
-      xiy = xM[1];
-      xiz = xM[2];
+      cache = &tip4p_cache[i];
+      xix = cache->xM[0];
+      xiy = cache->xM[1];
+      xiz = cache->xM[2];
     } else {
       xix = x[i].x;
       xiy = x[i].y;
@@ -326,10 +337,8 @@ void PPPMTIP4P::fieldforce_ik()
       }
     }
 
-    // convert E-field to force
-
     const double qfactor = qqrd2e * scale * q[i];
-    if (type[i] != typeO) {
+    if (cache == nullptr) {
       f[i].x += qfactor*ekx;
       f[i].y += qfactor*eky;
       if (zforce_active) f[i].z += qfactor*ekz;
@@ -343,13 +352,13 @@ void PPPMTIP4P::fieldforce_ik()
       f[i].y += fy * one_minus_alpha;
       if (zforce_active) f[i].z += fz * one_minus_alpha;
 
-      f[iH1].x += half_alpha * fx;
-      f[iH1].y += half_alpha * fy;
-      if (zforce_active) f[iH1].z += half_alpha * fz;
+      f[cache->iH1].x += half_alpha * fx;
+      f[cache->iH1].y += half_alpha * fy;
+      if (zforce_active) f[cache->iH1].z += half_alpha * fz;
 
-      f[iH2].x += half_alpha * fx;
-      f[iH2].y += half_alpha * fy;
-      if (zforce_active) f[iH2].z += half_alpha * fz;
+      f[cache->iH2].x += half_alpha * fx;
+      f[cache->iH2].y += half_alpha * fy;
+      if (zforce_active) f[cache->iH2].z += half_alpha * fz;
     }
   }
 }
diff --git a/src/KSPACE/pppm_tip4p.h b/src/KSPACE/pppm_tip4p.h
index abe3418959..d907f96dd3 100644
--- a/src/KSPACE/pppm_tip4p.h
+++ b/src/KSPACE/pppm_tip4p.h
@@ -39,13 +39,23 @@ class PPPMTIP4P : public PPPM {
   void slabcorr() override;
 
  private:
+  struct Tip4pCache {
+    int iH1;
+    int iH2;
+    int stamp;
+    double xM[3];
+  };
+
   int tip4p_cache_nmax;
-  int *tip4p_h1;
-  int *tip4p_h2;
-  int *tip4p_cache_valid;
-  double **tip4p_xm;
+  int tip4p_cache_stamp;
+  int tip4p_oxygen_count;
+  int tip4p_nonoxygen_count;
+  int *tip4p_oxygen;
+  int *tip4p_nonoxygen;
+  Tip4pCache *tip4p_cache;
 
   void grow_tip4p_cache();
+  void reset_tip4p_cache_stamps();
   void find_M_cached(int, int &, int &, double *&);
   void find_M(int, int &, int &, double *);
 };
