diff --git a/src/USER-OMP/dihedral_charmm_omp.cpp b/src/USER-OMP/dihedral_charmm_omp.cpp
index 00bb72c061..242a0a1d86 100644
--- a/src/USER-OMP/dihedral_charmm_omp.cpp
+++ b/src/USER-OMP/dihedral_charmm_omp.cpp
@@ -64,7 +64,7 @@ void DihedralCharmmOMP::compute(int eflag, int vflag)
     loop_setup_thr(ifrom, ito, tid, inum, nthreads);
     ThrData *thr = fix->get_thr(tid);
     thr->timer(Timer::START);
-    ev_setup_thr(eflag, vflag, nall, eatom, vatom, NULL, thr);
+    ev_setup_thr(eflag, vflag, nall, eatom, vatom, cvatom, thr);
 
     if (inum > 0) {
       if (evflag) {
diff --git a/src/USER-OMP/dihedral_class2_omp.cpp b/src/USER-OMP/dihedral_class2_omp.cpp
index 12b9303c2f..a54603cc39 100644
--- a/src/USER-OMP/dihedral_class2_omp.cpp
+++ b/src/USER-OMP/dihedral_class2_omp.cpp
@@ -58,7 +58,7 @@ void DihedralClass2OMP::compute(int eflag, int vflag)
     loop_setup_thr(ifrom, ito, tid, inum, nthreads);
     ThrData *thr = fix->get_thr(tid);
     thr->timer(Timer::START);
-    ev_setup_thr(eflag, vflag, nall, eatom, vatom, NULL, thr);
+    ev_setup_thr(eflag, vflag, nall, eatom, vatom, cvatom, thr);
 
     if (inum > 0) {
       if (evflag) {
diff --git a/src/USER-OMP/dihedral_cosine_shift_exp_omp.cpp b/src/USER-OMP/dihedral_cosine_shift_exp_omp.cpp
index 7161e2bb28..f42121f8a9 100644
--- a/src/USER-OMP/dihedral_cosine_shift_exp_omp.cpp
+++ b/src/USER-OMP/dihedral_cosine_shift_exp_omp.cpp
@@ -58,7 +58,7 @@ void DihedralCosineShiftExpOMP::compute(int eflag, int vflag)
     loop_setup_thr(ifrom, ito, tid, inum, nthreads);
     ThrData *thr = fix->get_thr(tid);
     thr->timer(Timer::START);
-    ev_setup_thr(eflag, vflag, nall, eatom, vatom, NULL, thr);
+    ev_setup_thr(eflag, vflag, nall, eatom, vatom, cvatom, thr);
 
     if (inum > 0) {
       if (evflag) {
diff --git a/src/USER-OMP/dihedral_fourier_omp.cpp b/src/USER-OMP/dihedral_fourier_omp.cpp
index c68c7987d0..cd12b3630e 100644
--- a/src/USER-OMP/dihedral_fourier_omp.cpp
+++ b/src/USER-OMP/dihedral_fourier_omp.cpp
@@ -57,7 +57,7 @@ void DihedralFourierOMP::compute(int eflag, int vflag)
     loop_setup_thr(ifrom, ito, tid, inum, nthreads);
     ThrData *thr = fix->get_thr(tid);
     thr->timer(Timer::START);
-    ev_setup_thr(eflag, vflag, nall, eatom, vatom, NULL, thr);
+    ev_setup_thr(eflag, vflag, nall, eatom, vatom, cvatom, thr);
 
     if (inum > 0) {
       if (evflag) {
diff --git a/src/USER-OMP/dihedral_harmonic_omp.cpp b/src/USER-OMP/dihedral_harmonic_omp.cpp
index 67963ec0ca..c3adb113e2 100644
--- a/src/USER-OMP/dihedral_harmonic_omp.cpp
+++ b/src/USER-OMP/dihedral_harmonic_omp.cpp
@@ -58,7 +58,7 @@ void DihedralHarmonicOMP::compute(int eflag, int vflag)
     loop_setup_thr(ifrom, ito, tid, inum, nthreads);
     ThrData *thr = fix->get_thr(tid);
     thr->timer(Timer::START);
-    ev_setup_thr(eflag, vflag, nall, eatom, vatom, NULL, thr);
+    ev_setup_thr(eflag, vflag, nall, eatom, vatom, cvatom, thr);
 
     if (inum > 0) {
       if (evflag) {
diff --git a/src/USER-OMP/dihedral_helix_omp.cpp b/src/USER-OMP/dihedral_helix_omp.cpp
index 13a112f887..b38ff2739a 100644
--- a/src/USER-OMP/dihedral_helix_omp.cpp
+++ b/src/USER-OMP/dihedral_helix_omp.cpp
@@ -61,7 +61,7 @@ void DihedralHelixOMP::compute(int eflag, int vflag)
     loop_setup_thr(ifrom, ito, tid, inum, nthreads);
     ThrData *thr = fix->get_thr(tid);
     thr->timer(Timer::START);
-    ev_setup_thr(eflag, vflag, nall, eatom, vatom, NULL, thr);
+    ev_setup_thr(eflag, vflag, nall, eatom, vatom, cvatom, thr);
 
     if (inum > 0) {
       if (evflag) {
diff --git a/src/USER-OMP/dihedral_multi_harmonic_omp.cpp b/src/USER-OMP/dihedral_multi_harmonic_omp.cpp
index 453d8b0fdd..7b79a63722 100644
--- a/src/USER-OMP/dihedral_multi_harmonic_omp.cpp
+++ b/src/USER-OMP/dihedral_multi_harmonic_omp.cpp
@@ -58,7 +58,7 @@ void DihedralMultiHarmonicOMP::compute(int eflag, int vflag)
     loop_setup_thr(ifrom, ito, tid, inum, nthreads);
     ThrData *thr = fix->get_thr(tid);
     thr->timer(Timer::START);
-    ev_setup_thr(eflag, vflag, nall, eatom, vatom, NULL, thr);
+    ev_setup_thr(eflag, vflag, nall, eatom, vatom, cvatom, thr);
 
     if (inum > 0) {
       if (evflag) {
diff --git a/src/USER-OMP/dihedral_nharmonic_omp.cpp b/src/USER-OMP/dihedral_nharmonic_omp.cpp
index 0f63b52d51..f3d8471c95 100644
--- a/src/USER-OMP/dihedral_nharmonic_omp.cpp
+++ b/src/USER-OMP/dihedral_nharmonic_omp.cpp
@@ -58,7 +58,7 @@ void DihedralNHarmonicOMP::compute(int eflag, int vflag)
     loop_setup_thr(ifrom, ito, tid, inum, nthreads);
     ThrData *thr = fix->get_thr(tid);
     thr->timer(Timer::START);
-    ev_setup_thr(eflag, vflag, nall, eatom, vatom, NULL, thr);
+    ev_setup_thr(eflag, vflag, nall, eatom, vatom, cvatom, thr);
 
     if (inum > 0) {
       if (evflag) {
diff --git a/src/USER-OMP/dihedral_opls_omp.cpp b/src/USER-OMP/dihedral_opls_omp.cpp
index 3b9d84751b..24cc4cd064 100644
--- a/src/USER-OMP/dihedral_opls_omp.cpp
+++ b/src/USER-OMP/dihedral_opls_omp.cpp
@@ -59,7 +59,7 @@ void DihedralOPLSOMP::compute(int eflag, int vflag)
     loop_setup_thr(ifrom, ito, tid, inum, nthreads);
     ThrData *thr = fix->get_thr(tid);
     thr->timer(Timer::START);
-    ev_setup_thr(eflag, vflag, nall, eatom, vatom, NULL, thr);
+    ev_setup_thr(eflag, vflag, nall, eatom, vatom, cvatom, thr);
 
     if (inum > 0) {
       if (evflag) {
diff --git a/src/USER-OMP/dihedral_quadratic_omp.cpp b/src/USER-OMP/dihedral_quadratic_omp.cpp
index 387b50827d..6f82c1e6b0 100644
--- a/src/USER-OMP/dihedral_quadratic_omp.cpp
+++ b/src/USER-OMP/dihedral_quadratic_omp.cpp
@@ -59,7 +59,7 @@ void DihedralQuadraticOMP::compute(int eflag, int vflag)
     loop_setup_thr(ifrom, ito, tid, inum, nthreads);
     ThrData *thr = fix->get_thr(tid);
     thr->timer(Timer::START);
-    ev_setup_thr(eflag, vflag, nall, eatom, vatom, NULL, thr);
+    ev_setup_thr(eflag, vflag, nall, eatom, vatom, cvatom, thr);
 
     if (inum > 0) {
       if (evflag) {
diff --git a/src/USER-OMP/dihedral_table_omp.cpp b/src/USER-OMP/dihedral_table_omp.cpp
index a20f7ebe13..a760fc6094 100644
--- a/src/USER-OMP/dihedral_table_omp.cpp
+++ b/src/USER-OMP/dihedral_table_omp.cpp
@@ -121,7 +121,7 @@ void DihedralTableOMP::compute(int eflag, int vflag)
     loop_setup_thr(ifrom, ito, tid, inum, nthreads);
     ThrData *thr = fix->get_thr(tid);
     thr->timer(Timer::START);
-    ev_setup_thr(eflag, vflag, nall, eatom, vatom, NULL, thr);
+    ev_setup_thr(eflag, vflag, nall, eatom, vatom, cvatom, thr);
 
     if (inum > 0) {
       if (evflag) {
diff --git a/src/USER-OMP/thr_data.h b/src/USER-OMP/thr_data.h
index 8000eb2563..6d8d09b6c9 100644
--- a/src/USER-OMP/thr_data.h
+++ b/src/USER-OMP/thr_data.h
@@ -98,6 +98,7 @@ class ThrData {
   double **vatom_imprp;
   double **vatom_kspce;
   double **cvatom_angle;
+  double **cvatom_dihed;
 
   // per thread segments of various force or similar arrays
 
diff --git a/src/USER-OMP/thr_omp.cpp b/src/USER-OMP/thr_omp.cpp
index 7d8f522f02..105aa80928 100644
--- a/src/USER-OMP/thr_omp.cpp
+++ b/src/USER-OMP/thr_omp.cpp
@@ -125,6 +125,11 @@ void ThrOMP::ev_setup_thr(int eflag, int vflag, int nall, double *eatom,
       if (nall > 0)
         memset(&(thr->vatom_dihed[0][0]),0,nall*6*sizeof(double));
     }
+    if (vflag & 8) {
+      thr->cvatom_dihed = cvatom + tid*nall;
+      if (nall > 0)
+        memset(&(thr->cvatom_dihed[0][0]),0,nall*9*sizeof(double));
+    }
   }
 
   if (thr_style & THR_IMPROPER) {
@@ -317,6 +322,9 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag,
       if (vflag & 4) {
         data_reduce_thr(&(dihedral->vatom[0][0]), nall, nthreads, 6, tid);
       }
+      if (vflag & 8) {
+        data_reduce_thr(&(dihedral->cvatom[0][0]), nall, nthreads, 9, tid);
+      }
 
     }
     break;
@@ -357,6 +365,9 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag,
         data_reduce_thr(&(dihedral->vatom[0][0]), nall, nthreads, 6, tid);
         data_reduce_thr(&(pair->vatom[0][0]), nall, nthreads, 6, tid);
       }
+      if (vflag & 8) {
+        data_reduce_thr(&(dihedral->cvatom[0][0]), nall, nthreads, 9, tid);
+      }
     }
     break;
 
@@ -1140,6 +1151,96 @@ void ThrOMP::ev_tally_thr(Dihedral * const dihed, const int i1, const int i2,
       }
     }
   }
+
+  // per-atom centroid virial
+  if (dihed->cvflag_atom) {
+    double f2[3], v1[9], v2[9], v3[9], v4[9];
+    double a1[3], a2[3], a3[3], a4[3];
+
+    // r0 = (r1+r2+r3+r4)/4
+    // rij = ri-rj
+    // total virial = r10*f1 + r20*f2 + r30*f3 + r40*f4
+    // vb1: r12
+    // vb2: r32
+    // vb3: r43
+
+    // a1 = r10 = (3*r12 - 2*r32 -   r43)/4
+    a1[0] = 0.25*(3*vb1x - 2*vb2x - vb3x);
+    a1[1] = 0.25*(3*vb1y - 2*vb2y - vb3y);
+    a1[2] = 0.25*(3*vb1z - 2*vb2z - vb3z);
+
+    // a2 = r20 = ( -r12 - 2*r32 -   r43)/4
+    a2[0] = 0.25*(-vb1x - 2*vb2x - vb3x);
+    a2[1] = 0.25*(-vb1y - 2*vb2y - vb3y);
+    a2[2] = 0.25*(-vb1z - 2*vb2z - vb3z);
+
+    // a3 = r30 = ( -r12 + 2*r32 -   r43)/4
+    a3[0] = 0.25*(-vb1x + 2*vb2x - vb3x);
+    a3[1] = 0.25*(-vb1y + 2*vb2y - vb3y);
+    a3[2] = 0.25*(-vb1z + 2*vb2z - vb3z);
+
+    // a4 = r40 = ( -r12 + 2*r32 + 3*r43)/4
+    a4[0] = 0.25*(-vb1x + 2*vb2x + 3*vb3x);
+    a4[1] = 0.25*(-vb1y + 2*vb2y + 3*vb3y);
+    a4[2] = 0.25*(-vb1z + 2*vb2z + 3*vb3z);
+
+    f2[0] = - f1[0] - f3[0] - f4[0];
+    f2[1] = - f1[1] - f3[1] - f4[1];
+    f2[2] = - f1[2] - f3[2] - f4[2];
+
+    v1[0] = a1[0]*f1[0];
+    v1[1] = a1[1]*f1[1];
+    v1[2] = a1[2]*f1[2];
+    v1[3] = a1[0]*f1[1];
+    v1[4] = a1[0]*f1[2];
+    v1[5] = a1[1]*f1[2];
+    v1[6] = a1[1]*f1[0];
+    v1[7] = a1[2]*f1[0];
+    v1[8] = a1[2]*f1[1];
+
+    v2[0] = a2[0]*f2[0];
+    v2[1] = a2[1]*f2[1];
+    v2[2] = a2[2]*f2[2];
+    v2[3] = a2[0]*f2[1];
+    v2[4] = a2[0]*f2[2];
+    v2[5] = a2[1]*f2[2];
+    v2[6] = a2[1]*f2[0];
+    v2[7] = a2[2]*f2[0];
+    v2[8] = a2[2]*f2[1];
+
+    v3[0] = a3[0]*f3[0];
+    v3[1] = a3[1]*f3[1];
+    v3[2] = a3[2]*f3[2];
+    v3[3] = a3[0]*f3[1];
+    v3[4] = a3[0]*f3[2];
+    v3[5] = a3[1]*f3[2];
+    v3[6] = a3[1]*f3[0];
+    v3[7] = a3[2]*f3[0];
+    v3[8] = a3[2]*f3[1];
+
+    v4[0] = a4[0]*f4[0];
+    v4[1] = a4[1]*f4[1];
+    v4[2] = a4[2]*f4[2];
+    v4[3] = a4[0]*f4[1];
+    v4[4] = a4[0]*f4[2];
+    v4[5] = a4[1]*f4[2];
+    v4[6] = a4[1]*f4[0];
+    v4[7] = a4[2]*f4[0];
+    v4[8] = a4[2]*f4[1];
+
+    if (newton_bond) {
+      v_tally9(thr->cvatom_dihed[i1],v1);
+      v_tally9(thr->cvatom_dihed[i2],v2);
+      v_tally9(thr->cvatom_dihed[i3],v3);
+      v_tally9(thr->cvatom_dihed[i4],v4);
+    } else {
+      if (i1 < nlocal) v_tally9(thr->cvatom_dihed[i1],v1);
+      if (i2 < nlocal) v_tally9(thr->cvatom_dihed[i2],v2);
+      if (i3 < nlocal) v_tally9(thr->cvatom_dihed[i3],v3);
+      if (i4 < nlocal) v_tally9(thr->cvatom_dihed[i4],v4);
+    }
+  }
+
 }
 
 /* ----------------------------------------------------------------------