diff --git a/src/group.cpp b/src/group.cpp
index 885ca76bac..de64c6615b 100644
--- a/src/group.cpp
+++ b/src/group.cpp
@@ -1240,6 +1240,87 @@ void Group::angmom(int igroup, double *cm, double *lmom, int iregion)
   MPI_Allreduce(p,lmom,3,MPI_DOUBLE,MPI_SUM,world);
 }
 
+/* ----------------------------------------------------------------------
+   compute the torque T (tq) on group
+   around center-of-mass cm
+   must unwrap atoms to compute T correctly
+------------------------------------------------------------------------- */
+
+void Group::torque(int igroup, double *cm, double *tq)
+{
+  int groupbit = bitmask[igroup];
+
+  double **x = atom->x;
+  double **f = atom->f;
+  int *mask = atom->mask;
+  int *image = atom->image;
+  int nlocal = atom->nlocal;
+
+  int xbox,ybox,zbox;
+  double dx,dy,dz;
+  double xprd = domain->xprd;
+  double yprd = domain->yprd;
+  double zprd = domain->zprd;
+  double tlocal[3];
+  tlocal[0] = tlocal[1] = tlocal[2] = 0.0;
+
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      xbox = (image[i] & 1023) - 512;
+      ybox = (image[i] >> 10 & 1023) - 512;
+      zbox = (image[i] >> 20) - 512;
+      dx = (x[i][0] + xbox*xprd) - cm[0];
+      dy = (x[i][1] + ybox*yprd) - cm[1];
+      dz = (x[i][2] + zbox*zprd) - cm[2];
+      tlocal[0] += dy*f[i][2] - dz*f[i][1];
+      tlocal[1] += dz*f[i][0] - dx*f[i][2];
+      tlocal[2] += dx*f[i][1] - dy*f[i][0];
+    }
+
+  MPI_Allreduce(tlocal,tq,3,MPI_DOUBLE,MPI_SUM,world);
+}
+
+/* ----------------------------------------------------------------------
+   compute the torque T (tq) on group of atoms in region
+   around center-of-mass cm
+   must unwrap atoms to compute T correctly
+------------------------------------------------------------------------- */
+
+void Group::torque(int igroup, double *cm, double *tq, int iregion)
+{
+  int groupbit = bitmask[igroup];
+  Region *region = domain->regions[iregion];
+
+  double **x = atom->x;
+  double **f = atom->f;
+  int *mask = atom->mask;
+  int *image = atom->image;
+  int nlocal = atom->nlocal;
+
+  int xbox,ybox,zbox;
+  double dx,dy,dz;
+  double xprd = domain->xprd;
+  double yprd = domain->yprd;
+  double zprd = domain->zprd;
+  double tlocal[3];
+  tlocal[0] = tlocal[1] = tlocal[2] = 0.0;
+
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
+      xbox = (image[i] & 1023) - 512;
+      ybox = (image[i] >> 10 & 1023) - 512;
+      zbox = (image[i] >> 20) - 512;
+      dx = (x[i][0] + xbox*xprd) - cm[0];
+      dy = (x[i][1] + ybox*yprd) - cm[1];
+      dz = (x[i][2] + zbox*zprd) - cm[2];
+      tlocal[0] += dy*f[i][2] - dz*f[i][1];
+      tlocal[1] += dz*f[i][0] - dx*f[i][2];
+      tlocal[2] += dx*f[i][1] - dy*f[i][0];
+    }
+
+  MPI_Allreduce(tlocal,tq,3,MPI_DOUBLE,MPI_SUM,world);
+}
+
 /* ----------------------------------------------------------------------
    compute moment of inertia tensor around center-of-mass cm of group
    must unwrap atoms to compute itensor correctly
diff --git a/src/group.h b/src/group.h
index 35e800280e..46ed702cad 100644
--- a/src/group.h
+++ b/src/group.h
@@ -54,6 +54,8 @@ class Group : protected Pointers {
   double gyration(int, double, double *, int);
   void angmom(int, double *, double *);    // angular momentum of group
   void angmom(int, double *, double *, int);
+  void torque(int, double *, double *);    // torque on group
+  void torque(int, double *, double *, int);
   void inertia(int, double *, double [3][3]);     // inertia tensor
   void inertia(int, double *, double [3][3], int);
   void omega(double *, double [3][3], double *);  // angular velocity
diff --git a/src/variable.cpp b/src/variable.cpp
index 126d11d5d8..1136593216 100644
--- a/src/variable.cpp
+++ b/src/variable.cpp
@@ -2423,8 +2423,8 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
    customize by adding a group function with optional region arg:
      count(group),mass(group),charge(group),
      xcm(group,dim),vcm(group,dim),fcm(group,dim),
-     bound(group,xmin),gyration(group),ke(group),angmom(group),
-     inertia(group,dim),omega(group,dim)
+     bound(group,xmin),gyration(group),ke(group),angmom(group,dim),
+     torque(group,dim),inertia(group,dim),omega(group,dim)
 ------------------------------------------------------------------------- */
 
 int Variable::group_function(char *word, char *contents, Tree **tree,
@@ -2438,7 +2438,8 @@ int Variable::group_function(char *word, char *contents, Tree **tree,
       strcmp(word,"vcm") && strcmp(word,"fcm") &&
       strcmp(word,"bound") && strcmp(word,"gyration") &&
       strcmp(word,"ke") && strcmp(word,"angmom") &&
-      strcmp(word,"inertia") && strcmp(word,"omega"))
+      strcmp(word,"torque") && strcmp(word,"inertia") && 
+      strcmp(word,"omega"))
     return 0;
 
   // parse contents for arg1,arg2,arg3 separated by commas
@@ -2588,6 +2589,24 @@ int Variable::group_function(char *word, char *contents, Tree **tree,
     else if (strcmp(arg2,"z") == 0) value = lmom[2];
     else error->all("Invalid group function in variable formula");
 
+  } else if (strcmp(word,"torque") == 0) {
+    atom->check_mass();
+    double xcm[3],tq[3];
+    if (narg == 2) {
+      double masstotal = group->mass(igroup);
+      group->xcm(igroup,masstotal,xcm);
+      group->torque(igroup,xcm,tq);
+    } else if (narg == 3) {
+      int iregion = region_function(arg3);
+      double masstotal = group->mass(igroup,iregion);
+      group->xcm(igroup,masstotal,xcm,iregion);
+      group->torque(igroup,xcm,tq,iregion);
+    } else error->all("Invalid group function in variable formula");
+    if (strcmp(arg2,"x") == 0) value = tq[0];
+    else if (strcmp(arg2,"y") == 0) value = tq[1];
+    else if (strcmp(arg2,"z") == 0) value = tq[2];
+    else error->all("Invalid group function in variable formula");
+
   } else if (strcmp(word,"inertia") == 0) {
     atom->check_mass();
     double xcm[3],inertia[3][3];