diff --git a/src/angle.cpp b/src/angle.cpp
index 1b9532ea32..ee7e855782 100644
--- a/src/angle.cpp
+++ b/src/angle.cpp
@@ -34,9 +34,10 @@ Angle::Angle(LAMMPS *lmp) : Pointers(lmp)
   allocated = 0;
   suffix_flag = Suffix::NONE;
 
-  maxeatom = maxvatom = 0;
+  maxeatom = maxvatom = maxcvatom = 0;
   eatom = NULL;
   vatom = NULL;
+  cvatom = NULL;
   setflag = NULL;
 
   execution_space = Host;
@@ -54,6 +55,7 @@ Angle::~Angle()
 
   memory->destroy(eatom);
   memory->destroy(vatom);
+  memory->destroy(cvatom);
 }
 
 /* ----------------------------------------------------------------------
@@ -85,9 +87,10 @@ void Angle::ev_setup(int eflag, int vflag, int alloc)
   eflag_global = eflag % 2;
   eflag_atom = eflag / 2;
 
-  vflag_either = vflag;
   vflag_global = vflag % 4;
-  vflag_atom = vflag / 4;
+  vflag_atom = vflag & 4;
+  cvflag_atom = vflag & 8;
+  vflag_either = vflag_global || vflag_atom;
 
   // reallocate per-atom arrays if necessary
 
@@ -105,6 +108,13 @@ void Angle::ev_setup(int eflag, int vflag, int alloc)
       memory->create(vatom,comm->nthreads*maxvatom,6,"angle:vatom");
     }
   }
+  if (cvflag_atom && atom->nmax > maxcvatom) {
+    maxcvatom = atom->nmax;
+    if (alloc) {
+      memory->destroy(cvatom);
+      memory->create(cvatom,comm->nthreads*maxcvatom,9,"angle:cvatom");
+    }
+  }
 
   // zero accumulators
 
@@ -127,6 +137,22 @@ void Angle::ev_setup(int eflag, int vflag, int alloc)
       vatom[i][5] = 0.0;
     }
   }
+  if (cvflag_atom && alloc) {
+    n = atom->nlocal;
+    if (force->newton_bond) n += atom->nghost;
+    for (i = 0; i < n; i++) {
+      cvatom[i][0] = 0.0;
+      cvatom[i][1] = 0.0;
+      cvatom[i][2] = 0.0;
+      cvatom[i][3] = 0.0;
+      cvatom[i][4] = 0.0;
+      cvatom[i][5] = 0.0;
+      cvatom[i][6] = 0.0;
+      cvatom[i][7] = 0.0;
+      cvatom[i][8] = 0.0;
+      cvatom[i][9] = 0.0;
+    }
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -230,6 +256,76 @@ void Angle::ev_tally(int i, int j, int k, int nlocal, int newton_bond,
       }
     }
   }
+
+  // per-atom centroid virial
+  if (cvflag_atom) {
+
+    // r0 = (r1+r2+r3)/3
+    // rij = ri-rj
+    // total virial = r10*f1 + r20*f2 + r30*f3
+    // del1: r12
+    // del2: r32
+
+    if (newton_bond || i < nlocal) {
+      double a1[3];
+
+      // a1 = r10 = (2*r12 -   r32)/3
+      a1[0] = THIRD*(2*delx1-delx2);
+      a1[1] = THIRD*(2*dely1-dely2);
+      a1[2] = THIRD*(2*delz1-delz2);
+
+      cvatom[i][0] += a1[0]*f1[0];
+      cvatom[i][1] += a1[1]*f1[1];
+      cvatom[i][2] += a1[2]*f1[2];
+      cvatom[i][3] += a1[0]*f1[1];
+      cvatom[i][4] += a1[0]*f1[2];
+      cvatom[i][5] += a1[1]*f1[2];
+      cvatom[i][6] += a1[1]*f1[0];
+      cvatom[i][7] += a1[2]*f1[0];
+      cvatom[i][8] += a1[2]*f1[1];
+    }
+    if (newton_bond || j < nlocal) {
+      double a2[3];
+      double f2[3];
+
+      // a2 = r20 = ( -r12 -   r32)/3
+      a2[0] = THIRD*(-delx1-delx2);
+      a2[1] = THIRD*(-dely1-dely2);
+      a2[2] = THIRD*(-delz1-delz2);
+
+      f2[0] = - f1[0] - f3[0];
+      f2[1] = - f1[1] - f3[1];
+      f2[2] = - f1[2] - f3[2];
+
+      cvatom[j][0] += a2[0]*f2[0];
+      cvatom[j][1] += a2[1]*f2[1];
+      cvatom[j][2] += a2[2]*f2[2];
+      cvatom[j][3] += a2[0]*f2[1];
+      cvatom[j][4] += a2[0]*f2[2];
+      cvatom[j][5] += a2[1]*f2[2];
+      cvatom[j][6] += a2[1]*f2[0];
+      cvatom[j][7] += a2[2]*f2[0];
+      cvatom[j][8] += a2[2]*f2[1];
+    }
+    if (newton_bond || k < nlocal) {
+      double a3[3];
+
+      // a3 = r30 = ( -r12 + 2*r32)/3
+      a3[0] = THIRD*(-delx1+2*delx2);
+      a3[1] = THIRD*(-dely1+2*dely2);
+      a3[2] = THIRD*(-delz1+2*delz2);
+
+      cvatom[k][0] += a3[0]*f3[0];
+      cvatom[k][1] += a3[1]*f3[1];
+      cvatom[k][2] += a3[2]*f3[2];
+      cvatom[k][3] += a3[0]*f3[1];
+      cvatom[k][4] += a3[0]*f3[2];
+      cvatom[k][5] += a3[1]*f3[2];
+      cvatom[k][6] += a3[1]*f3[0];
+      cvatom[k][7] += a3[2]*f3[0];
+      cvatom[k][8] += a3[2]*f3[1];
+    }
+  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -238,5 +334,6 @@ double Angle::memory_usage()
 {
   double bytes = comm->nthreads*maxeatom * sizeof(double);
   bytes += comm->nthreads*maxvatom*6 * sizeof(double);
+  bytes += comm->nthreads*maxcvatom*9 * sizeof(double);
   return bytes;
 }
diff --git a/src/angle.h b/src/angle.h
index c0d1199dcd..139380ff39 100644
--- a/src/angle.h
+++ b/src/angle.h
@@ -28,6 +28,7 @@ class Angle : protected Pointers {
   double energy;                  // accumulated energies
   double virial[6];               // accumulated virial
   double *eatom,**vatom;          // accumulated per-atom energy/virial
+  double **cvatom;                // accumulated per-atom centroid virial
 
   // KOKKOS host/device flag and data masks
 
@@ -56,12 +57,12 @@ class Angle : protected Pointers {
 
   int evflag;
   int eflag_either,eflag_global,eflag_atom;
-  int vflag_either,vflag_global,vflag_atom;
-  int maxeatom,maxvatom;
+  int vflag_either,vflag_global,vflag_atom,cvflag_atom;
+  int maxeatom,maxvatom,maxcvatom;
 
   void ev_init(int eflag, int vflag, int alloc = 1) {
     if (eflag||vflag) ev_setup(eflag, vflag, alloc);
-    else evflag = eflag_either = eflag_global = eflag_atom = vflag_either = vflag_global = vflag_atom = 0;
+    else evflag = eflag_either = eflag_global = eflag_atom = vflag_either = vflag_global = vflag_atom = cvflag_atom = 0;
   }
   void ev_setup(int, int, int alloc = 1);
   void ev_tally(int, int, int, int, int, double, double *, double *,