diff --git a/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp b/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp
index 3a3c2b11e1..76ae162f65 100644
--- a/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp
+++ b/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp
@@ -13,6 +13,7 @@
 
 /* ----------------------------------------------------------------------
    Contributing authors: Amalie Frischknecht and Ahmed Ismail (SNL)
+   simpler force assignment added by Rolf Isele-Holder (Aachen University)
 ------------------------------------------------------------------------- */
 
 #include "math.h"
@@ -75,7 +76,7 @@ PairLJCutCoulLongTIP4P::~PairLJCutCoulLongTIP4P()
 
 void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag)
 {
-  int i,j,ii,jj,inum,jnum,itype,jtype,itable;
+  int i,j,ii,jj,inum,jnum,itype,jtype,itable,key;
   int n,vlist[6];
   int iH1,iH2,jH1,jH2;
   double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul;
@@ -277,6 +278,7 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag)
           // vlist stores 2,4,6 atoms whose forces contribute to virial
 
           n = 0;
+          key = 0;
 
           if (itype != typeO) {
             f[i][0] += delx * cforce;
@@ -290,33 +292,23 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag)
               v[3] = x[i][0] * dely * cforce;
               v[4] = x[i][0] * delz * cforce;
               v[5] = x[i][1] * delz * cforce;
-              vlist[n++] = i;
             }
+            vlist[n++] = i;
 
           } else {
+            key++;
 
             fd[0] = delx*cforce;
             fd[1] = dely*cforce;
             fd[2] = delz*cforce;
 
-            delxOM = x[i][0] - x1[0];
-            delyOM = x[i][1] - x1[1];
-            delzOM = x[i][2] - x1[2];
+            fO[0] = fd[0]*(1 - alpha);
+            fO[1] = fd[1]*(1 - alpha);
+            fO[2] = fd[2]*(1 - alpha);
 
-            ddotf = (delxOM * fd[0] + delyOM * fd[1] + delzOM * fd[2]) /
-              (qdist*qdist);
-
-            f1[0] = ddotf * delxOM;
-            f1[1] = ddotf * delyOM;
-            f1[2] = ddotf * delzOM;
-
-            fO[0] = fd[0] - alpha * (fd[0] - f1[0]);
-            fO[1] = fd[1] - alpha * (fd[1] - f1[1]);
-            fO[2] = fd[2] - alpha * (fd[2] - f1[2]);
-
-            fH[0] = 0.5 * alpha * (fd[0] - f1[0]);
-            fH[1] = 0.5 * alpha * (fd[1] - f1[1]);
-            fH[2] = 0.5 * alpha * (fd[2] - f1[2]);
+            fH[0] = 0.5 * alpha * fd[0];
+            fH[1] = 0.5 * alpha * fd[1];
+            fH[2] = 0.5 * alpha * fd[2];
 
             f[i][0] += fO[0];
             f[i][1] += fO[1];
@@ -340,11 +332,10 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag)
               v[3] = x[i][0]*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1];
               v[4] = x[i][0]*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2];
               v[5] = x[i][1]*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2];
-
-              vlist[n++] = i;
-              vlist[n++] = iH1;
-              vlist[n++] = iH2;
             }
+            vlist[n++] = i;
+            vlist[n++] = iH1;
+            vlist[n++] = iH2;
           }
 
           if (jtype != typeO) {
@@ -359,33 +350,23 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag)
               v[3] -= x[j][0] * dely * cforce;
               v[4] -= x[j][0] * delz * cforce;
               v[5] -= x[j][1] * delz * cforce;
-              vlist[n++] = j;
             }
+            vlist[n++] = j;
 
           } else {
+            key += 2;
 
             fd[0] = -delx*cforce;
             fd[1] = -dely*cforce;
             fd[2] = -delz*cforce;
 
-            delxOM = x[j][0] - x2[0];
-            delyOM = x[j][1] - x2[1];
-            delzOM = x[j][2] - x2[2];
+            fO[0] = fd[0]*(1 - alpha);
+            fO[1] = fd[1]*(1 - alpha);
+            fO[2] = fd[2]*(1 - alpha);
 
-            ddotf = (delxOM * fd[0] + delyOM * fd[1] + delzOM * fd[2]) /
-              (qdist*qdist);
-
-            f1[0] = ddotf * delxOM;
-            f1[1] = ddotf * delyOM;
-            f1[2] = ddotf * delzOM;
-
-            fO[0] = fd[0] - alpha * (fd[0] - f1[0]);
-            fO[1] = fd[1] - alpha * (fd[1] - f1[1]);
-            fO[2] = fd[2] - alpha * (fd[2] - f1[2]);
-
-            fH[0] = 0.5 * alpha * (fd[0] - f1[0]);
-            fH[1] = 0.5 * alpha * (fd[1] - f1[1]);
-            fH[2] = 0.5 * alpha * (fd[2] - f1[2]);
+            fH[0] = 0.5 * alpha * fd[0];
+            fH[1] = 0.5 * alpha * fd[1];
+            fH[2] = 0.5 * alpha * fd[2];
 
             f[j][0] += fO[0];
             f[j][1] += fO[1];
@@ -409,11 +390,10 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag)
               v[3] += x[j][0]*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1];
               v[4] += x[j][0]*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2];
               v[5] += x[j][1]*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2];
-
-              vlist[n++] = j;
-              vlist[n++] = jH1;
-              vlist[n++] = jH2;
             }
+            vlist[n++] = j;
+            vlist[n++] = jH1;
+            vlist[n++] = jH2;
           }
 
           if (eflag) {
@@ -426,7 +406,7 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag)
             if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
           } else ecoul = 0.0;
 
-          if (evflag) ev_tally_list(n,vlist,ecoul,v);
+          if (evflag) ev_tally_list(ecoul,vlist,v,alpha,key);
         }
       }
     }
@@ -614,3 +594,4 @@ double PairLJCutCoulLongTIP4P::memory_usage()
   bytes += 2 * nmax * sizeof(double);
   return bytes;
 }
+
diff --git a/src/KSPACE/pppm_tip4p.cpp b/src/KSPACE/pppm_tip4p.cpp
index d7a8450db6..1146a64d27 100644
--- a/src/KSPACE/pppm_tip4p.cpp
+++ b/src/KSPACE/pppm_tip4p.cpp
@@ -22,8 +22,10 @@
 #include "force.h"
 #include "memory.h"
 #include "error.h"
+#include "math_const.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define OFFSET 16384
 
@@ -38,7 +40,7 @@ using namespace LAMMPS_NS;
 /* ---------------------------------------------------------------------- */
 
 PPPMTIP4P::PPPMTIP4P(LAMMPS *lmp, int narg, char **arg) :
-  PPPMOld(lmp, narg, arg) {}
+  PPPM(lmp, narg, arg) {}
 
 /* ---------------------------------------------------------------------- */
 
@@ -49,7 +51,7 @@ void PPPMTIP4P::init()
   if (force->newton == 0)
     error->all(FLERR,"Kspace style pppm/tip4p requires newton on");
 
-  PPPMOld::init();
+  PPPM::init();
 }
 
 /* ----------------------------------------------------------------------
@@ -162,6 +164,16 @@ void PPPMTIP4P::make_rho()
 ------------------------------------------------------------------------- */
 
 void PPPMTIP4P::fieldforce()
+{
+  if (differentiation_flag == 1) fieldforce_ad();
+  else fieldforce_ik();
+}
+
+/* ----------------------------------------------------------------------
+   interpolate from grid to get electric field & force on my particles for ik
+------------------------------------------------------------------------- */
+
+void PPPMTIP4P::fieldforce_ik()
 {
   int i,l,m,n,nx,ny,nz,mx,my,mz;
   FFT_SCALAR dx,dy,dz,x0,y0,z0;
@@ -170,14 +182,12 @@ void PPPMTIP4P::fieldforce()
   int iH1,iH2;
   double xM[3];
   double fx,fy,fz;
-  double ddotf, rOMx, rOMy, rOMz, f1x, f1y, f1z;
 
   // loop over my charges, interpolate electric field from nearby grid points
   // (nx,ny,nz) = global coords of grid pt to "lower left" of charge
   // (dx,dy,dz) = distance to "lower left" grid pt
   // (mx,my,mz) = global coords of moving stencil pt
   // ek = 3 components of E-field on particle
-
   double *q = atom->q;
   double **x = atom->x;
   double **f = atom->f;
@@ -231,27 +241,245 @@ void PPPMTIP4P::fieldforce()
       fz = qfactor * ekz;
       find_M(i,iH1,iH2,xM);
 
-      rOMx = xM[0] - x[i][0];
-      rOMy = xM[1] - x[i][1];
-      rOMz = xM[2] - x[i][2];
+      f[i][0] += fx*(1 - alpha);
+      f[i][1] += fy*(1 - alpha);
+      f[i][2] += fz*(1 - alpha);
 
-      ddotf = (rOMx * fx + rOMy * fy + rOMz * fz) / (qdist * qdist);
+      f[iH1][0] += 0.5*alpha*fx;
+      f[iH1][1] += 0.5*alpha*fy;
+      f[iH1][2] += 0.5*alpha*fz;
 
-      f1x = ddotf * rOMx;
-      f1y = ddotf * rOMy;
-      f1z = ddotf * rOMz;
+      f[iH2][0] += 0.5*alpha*fx;
+      f[iH2][1] += 0.5*alpha*fy;
+      f[iH2][2] += 0.5*alpha*fz;
+    }
+  }
+}
 
-      f[i][0] += fx - alpha * (fx - f1x);
-      f[i][1] += fy - alpha * (fy - f1y);
-      f[i][2] += fz - alpha * (fz - f1z);
+/* ----------------------------------------------------------------------
+   interpolate from grid to get electric field & force on my particles for ad
+------------------------------------------------------------------------- */
 
-      f[iH1][0] += 0.5*alpha*(fx - f1x);
-      f[iH1][1] += 0.5*alpha*(fy - f1y);
-      f[iH1][2] += 0.5*alpha*(fz - f1z);
+void PPPMTIP4P::fieldforce_ad()
+{
+  int i,l,m,n,nx,ny,nz,mx,my,mz;
+  FFT_SCALAR dx,dy,dz,x0,y0,z0,dx0,dy0,dz0;
+  FFT_SCALAR ekx,eky,ekz;
+  double *xi;
+  int iH1,iH2;
+  double xM[3];
+  double s1,s2,s3;
+  double sf;
+  double *prd;
+  double fx,fy,fz;
 
-      f[iH2][0] += 0.5*alpha*(fx - f1x);
-      f[iH2][1] += 0.5*alpha*(fy - f1y);
-      f[iH2][2] += 0.5*alpha*(fz - f1z);
+  if (triclinic == 0) prd = domain->prd;
+  else prd = domain->prd_lamda;
+
+  double xprd = prd[0];
+  double yprd = prd[1];
+  double zprd = prd[2];
+  double zprd_slab = zprd*slab_volfactor;
+
+  double hx_inv = nx_pppm/xprd;
+  double hy_inv = ny_pppm/yprd;
+  double hz_inv = nz_pppm/zprd;
+  
+
+
+  // loop over my charges, interpolate electric field from nearby grid points
+  // (nx,ny,nz) = global coords of grid pt to "lower left" of charge
+  // (dx,dy,dz) = distance to "lower left" grid pt
+  // (mx,my,mz) = global coords of moving stencil pt
+  // ek = 3 components of E-field on particle
+
+  double *q = atom->q;
+  double **x = atom->x;
+  double **f = atom->f;
+
+  int *type = atom->type;
+  int nlocal = atom->nlocal;
+
+  for (i = 0; i < nlocal; i++) {
+    if (type[i] == typeO) {
+      find_M(i,iH1,iH2,xM);      
+      xi = xM;
+    } else xi = x[i];
+
+    nx = part2grid[i][0];
+    ny = part2grid[i][1];
+    nz = part2grid[i][2];
+    dx = nx+shiftone - (x[i][0]-boxlo[0])*delxinv;
+    dy = ny+shiftone - (x[i][1]-boxlo[1])*delyinv;
+    dz = nz+shiftone - (x[i][2]-boxlo[2])*delzinv;
+    
+    compute_rho1d(dx,dy,dz);
+    compute_drho1d(dx,dy,dz);
+
+    ekx = eky = ekz = ZEROF;
+    for (n = nlower; n <= nupper; n++) {
+      mz = n+nz;
+      for (m = nlower; m <= nupper; m++) {
+        my = m+ny;
+        for (l = nlower; l <= nupper; l++) {
+          mx = l+nx;
+          ekx += drho1d[0][l]*rho1d[1][m]*rho1d[2][n]*u_brick[mz][my][mx];
+          eky += rho1d[0][l]*drho1d[1][m]*rho1d[2][n]*u_brick[mz][my][mx];
+          ekz += rho1d[0][l]*rho1d[1][m]*drho1d[2][n]*u_brick[mz][my][mx];
+        }
+      }
+    }
+
+    ekx *= hx_inv;
+    eky *= hy_inv;
+    ekz *= hz_inv;
+    
+    // convert E-field to force and substract self forces
+    const double qfactor = force->qqrd2e * scale;
+
+    s1 = x[i][0]*hx_inv;
+    s2 = x[i][1]*hy_inv;
+    s3 = x[i][2]*hz_inv;
+    sf = sf_coeff[0]*sin(2*MY_PI*s1);
+    sf += sf_coeff[1]*sin(4*MY_PI*s1);
+    sf *= 2*q[i]*q[i];
+    fx = qfactor*(ekx*q[i] - sf);
+
+    sf = sf_coeff[2]*sin(2*MY_PI*s2);
+    sf += sf_coeff[3]*sin(4*MY_PI*s2);
+    sf *= 2*q[i]*q[i];
+    fy = qfactor*(eky*q[i] - sf);
+
+    sf = sf_coeff[4]*sin(2*MY_PI*s3);
+    sf += sf_coeff[5]*sin(4*MY_PI*s3);
+    sf *= 2*q[i]*q[i];
+    fz = qfactor*(ekz*q[i] - sf);
+
+    if (type[i] != typeO) {
+      f[i][0] += fx;
+      f[i][1] += fy;
+      f[i][2] += fz;
+
+    } else {
+      find_M(i,iH1,iH2,xM);
+
+      f[i][0] += fx*(1 - alpha);
+      f[i][1] += fy*(1 - alpha);
+      f[i][2] += fz*(1 - alpha);
+
+      f[iH1][0] += 0.5*alpha*fx;
+      f[iH1][1] += 0.5*alpha*fy;
+      f[iH1][2] += 0.5*alpha*fz;
+
+      f[iH2][0] += 0.5*alpha*fx;
+      f[iH2][1] += 0.5*alpha*fy;
+      f[iH2][2] += 0.5*alpha*fz;
+    }
+  }
+}
+
+
+/* ----------------------------------------------------------------------
+   interpolate from grid to get electric field & force on my particles 
+------------------------------------------------------------------------- */
+
+void PPPMTIP4P::fieldforce_peratom()
+{
+  int i,l,m,n,nx,ny,nz,mx,my,mz;
+  FFT_SCALAR dx,dy,dz,x0,y0,z0;
+  double *xi;
+  int iH1,iH2;
+  double xM[3];
+  FFT_SCALAR u_pa,v0,v1,v2,v3,v4,v5;
+
+  // loop over my charges, interpolate electric field from nearby grid points
+  // (nx,ny,nz) = global coords of grid pt to "lower left" of charge
+  // (dx,dy,dz) = distance to "lower left" grid pt
+  // (mx,my,mz) = global coords of moving stencil pt
+  // ek = 3 components of E-field on particle
+  double *q = atom->q;
+  double **x = atom->x;
+  double **f = atom->f;
+
+  int *type = atom->type;
+  int nlocal = atom->nlocal;
+
+  for (i = 0; i < nlocal; i++) {
+    if (type[i] == typeO) {
+      find_M(i,iH1,iH2,xM);      
+      xi = xM;
+    } else xi = x[i];
+
+    nx = part2grid[i][0];
+    ny = part2grid[i][1];
+    nz = part2grid[i][2];
+    dx = nx+shiftone - (xi[0]-boxlo[0])*delxinv;
+    dy = ny+shiftone - (xi[1]-boxlo[1])*delyinv;
+    dz = nz+shiftone - (xi[2]-boxlo[2])*delzinv;
+
+    compute_rho1d(dx,dy,dz);
+
+    u_pa = v0 = v1 = v2 = v3 = v4 = v5 = ZEROF;
+    for (n = nlower; n <= nupper; n++) {
+      mz = n+nz;
+      z0 = rho1d[2][n];
+      for (m = nlower; m <= nupper; m++) {
+        my = m+ny;
+        y0 = z0*rho1d[1][m];
+        for (l = nlower; l <= nupper; l++) {
+          mx = l+nx;
+          x0 = y0*rho1d[0][l];
+          if (eflag_atom) u_pa += x0*u_brick[mz][my][mx];	
+          if (vflag_atom) {
+            v0 += x0*v0_brick[mz][my][mx];
+            v1 += x0*v1_brick[mz][my][mx];
+            v2 += x0*v2_brick[mz][my][mx];
+            v3 += x0*v3_brick[mz][my][mx];
+            v4 += x0*v4_brick[mz][my][mx];
+            v5 += x0*v5_brick[mz][my][mx];
+          }
+        }
+      }
+    }
+
+
+    if (eflag_atom) {
+      if (type[i] != typeO) {
+        eatom[i] += q[i]*u_pa;
+      } else {
+        eatom[i] += q[i]*u_pa*(1-alpha);
+        eatom[iH1] += q[i]*u_pa*alpha*0.5;
+        eatom[iH2] += q[i]*u_pa*alpha*0.5;
+      }
+    }
+    if (vflag_atom) {
+      if (type[i] != typeO) {
+        vatom[i][0] += v0*q[i];
+        vatom[i][1] += v1*q[i];
+        vatom[i][2] += v2*q[i];
+        vatom[i][3] += v3*q[i];
+        vatom[i][4] += v4*q[i];
+        vatom[i][5] += v5*q[i];
+      } else {
+        vatom[i][0] += v0*(1-alpha)*q[i];
+        vatom[i][1] += v1*(1-alpha)*q[i];
+        vatom[i][2] += v2*(1-alpha)*q[i];
+        vatom[i][3] += v3*(1-alpha)*q[i];
+        vatom[i][4] += v4*(1-alpha)*q[i];
+        vatom[i][5] += v5*(1-alpha)*q[i];
+        vatom[iH1][0] += v0*alpha*0.5*q[i];
+        vatom[iH1][1] += v1*alpha*0.5*q[i];
+        vatom[iH1][2] += v2*alpha*0.5*q[i];
+        vatom[iH1][3] += v3*alpha*0.5*q[i];
+        vatom[iH1][4] += v4*alpha*0.5*q[i];
+        vatom[iH1][5] += v5*alpha*0.5*q[i];
+        vatom[iH2][0] += v0*alpha*0.5*q[i];
+        vatom[iH2][1] += v1*alpha*0.5*q[i];
+        vatom[iH2][2] += v2*alpha*0.5*q[i];
+        vatom[iH2][3] += v3*alpha*0.5*q[i];
+        vatom[iH2][4] += v4*alpha*0.5*q[i];
+        vatom[iH2][5] += v5*alpha*0.5*q[i];
+      }
     }
   }
 }
diff --git a/src/KSPACE/pppm_tip4p.h b/src/KSPACE/pppm_tip4p.h
index df98002d6c..6726af2ad6 100644
--- a/src/KSPACE/pppm_tip4p.h
+++ b/src/KSPACE/pppm_tip4p.h
@@ -20,11 +20,11 @@ KSpaceStyle(pppm/tip4p,PPPMTIP4P)
 #ifndef LMP_PPPM_TIP4P_H
 #define LMP_PPPM_TIP4P_H
 
-#include "pppm_old.h"
+#include "pppm.h"
 
 namespace LAMMPS_NS {
 
-class PPPMTIP4P : public PPPMOld {
+class PPPMTIP4P : public PPPM {
  public:
   PPPMTIP4P(class LAMMPS *, int, char **);
   virtual ~PPPMTIP4P () {};
@@ -34,6 +34,9 @@ class PPPMTIP4P : public PPPMOld {
   virtual void particle_map();
   virtual void make_rho();
   virtual void fieldforce();
+  void fieldforce_ik();
+  void fieldforce_ad();
+  virtual void fieldforce_peratom();
 
  private:
   void find_M(int, int &, int &, double *);
diff --git a/src/pair.cpp b/src/pair.cpp
index 33a8fcafc6..90732b010a 100644
--- a/src/pair.cpp
+++ b/src/pair.cpp
@@ -755,6 +755,8 @@ void Pair::ev_tally4(int i, int j, int k, int m, double evdwl,
 }
 
 /* ----------------------------------------------------------------------
+   NOTE: this is older version of ev tally for TIP4P
+   NOTE: delete this when all pair TIP4P derivatives have been updated
    tally ecoul and virial into each of n atoms in list
    called by TIP4P potential, newton_pair is always on
    changes v values by dividing by n
@@ -802,6 +804,87 @@ void Pair::ev_tally_list(int n, int *list, double ecoul, double *v)
   }
 }
 
+/* ----------------------------------------------------------------------
+   tally ecoul and virial into each of atoms in list, list length set by key
+   called by TIP4P potential, newton_pair is always on
+ ------------------------------------------------------------------------- */
+
+void Pair::ev_tally_list(double ecoul, int *list, double *v,
+                         double alpha, int key)
+{
+  int i,j;
+
+  if (eflag_either) {
+    if (eflag_global) eng_coul += ecoul;
+    if (eflag_atom) {
+      if (key == 0) {
+        eatom[list[0]] += 0.5*ecoul;
+        eatom[list[1]] += 0.5*ecoul;
+      } else if (key == 1) {
+        eatom[list[0]] += 0.5*ecoul*(1-alpha);
+        eatom[list[1]] += 0.25*ecoul*alpha;
+        eatom[list[2]] += 0.25*ecoul*alpha;
+        eatom[list[3]] += 0.5*ecoul;
+      } else if (key == 2) {
+        eatom[list[0]] += 0.5*ecoul;
+        eatom[list[1]] += 0.5*ecoul*(1-alpha);
+        eatom[list[2]] += 0.25*ecoul*alpha;
+        eatom[list[3]] += 0.25*ecoul*alpha;
+      } else {
+        eatom[list[0]] += 0.5*ecoul*(1-alpha);
+        eatom[list[1]] += 0.25*ecoul*alpha;
+        eatom[list[2]] += 0.25*ecoul*alpha;
+        eatom[list[3]] += 0.5*ecoul*(1-alpha);
+        eatom[list[4]] += 0.25*ecoul*alpha;
+        eatom[list[5]] += 0.25*ecoul*alpha;
+      }
+    }
+  }
+
+  if (vflag_either) {
+    if (vflag_global) {
+      virial[0] += v[0];
+      virial[1] += v[1];
+      virial[2] += v[2];
+      virial[3] += v[3];
+      virial[4] += v[4];
+      virial[5] += v[5];
+    }
+
+    if (vflag_atom) {
+      if (key == 0) {
+        for (i = 0; i <= 5; i++) {
+          vatom[list[0]][i] += 0.5*v[i];
+          vatom[list[1]][i] += 0.5*v[i];
+        }
+      } else if (key == 1) {
+        for (i = 0; i <= 5; i++) {
+          vatom[list[0]][i] += 0.5*v[i]*(1-alpha);
+          vatom[list[1]][i] += 0.25*v[i]*alpha;
+          vatom[list[2]][i] += 0.25*v[i]*alpha;
+          vatom[list[3]][i] += 0.5*v[i];
+        }
+      } else if (key == 2) {
+        for (i = 0; i <= 5; i++) {
+          vatom[list[0]][i] += 0.5*v[i];
+          vatom[list[1]][i] += 0.5*v[i]*(1-alpha);
+          vatom[list[2]][i] += 0.25*v[i]*alpha;
+          vatom[list[3]][i] += 0.25*v[i]*alpha;
+        }
+      } else {
+        for (i = 0; i <= 5; i++) {
+          vatom[list[0]][i] += 0.5*v[i]*(1-alpha);
+          vatom[list[1]][i] += 0.25*v[i]*alpha;
+          vatom[list[2]][i] += 0.25*v[i]*alpha;
+          vatom[list[3]][i] += 0.5*v[i]*(1-alpha);
+          vatom[list[4]][i] += 0.25*v[i]*alpha;
+          vatom[list[5]][i] += 0.25*v[i]*alpha;
+        }
+      }
+    }
+  }
+}
+
 /* ----------------------------------------------------------------------
    tally virial into per-atom accumulators
    called by REAX/C potential, newton_pair is always on
diff --git a/src/pair.h b/src/pair.h
index 011478547e..ac8db364d5 100644
--- a/src/pair.h
+++ b/src/pair.h
@@ -173,6 +173,7 @@ class Pair : protected Pointers {
   void ev_tally4(int, int, int, int, double,
                  double *, double *, double *, double *, double *, double *);
   void ev_tally_list(int, int *, double, double *);
+  void ev_tally_list(double, int *, double *, double, int);
   void v_tally2(int, int, double, double *);
   void v_tally_tensor(int, int, int, int,
                       double, double, double, double, double, double);