git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@4199 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/group.cpp

@@ -1241,7 +1241,7 @@ void Group::angmom(int igroup, double *cm, double *lmom, int iregion)
 }
 
 /* ----------------------------------------------------------------------
-   compute moment of inertia tensor around center-of-mass cm
+   compute moment of inertia tensor around center-of-mass cm of group
    must unwrap atoms to compute itensor correctly
 ------------------------------------------------------------------------- */
 
@@ -1293,6 +1293,60 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3])
   MPI_Allreduce(&ione[0][0],&itensor[0][0],9,MPI_DOUBLE,MPI_SUM,world);
 }
 
+/* ----------------------------------------------------------------------
+   compute moment of inertia tensor around cm of group of atoms in region
+   must unwrap atoms to compute itensor correctly
+------------------------------------------------------------------------- */
+
+void Group::inertia(int igroup, double *cm, double itensor[3][3], int iregion)
+{
+  int i,j;
+
+  int groupbit = bitmask[igroup];
+  Region *region = domain->regions[iregion];
+
+  double **x = atom->x;
+  int *mask = atom->mask;
+  int *type = atom->type;
+  int *image = atom->image;
+  double *mass = atom->mass;
+  double *rmass = atom->rmass;
+  int nlocal = atom->nlocal;
+
+  int xbox,ybox,zbox;
+  double dx,dy,dz,massone;
+  double xprd = domain->xprd;
+  double yprd = domain->yprd;
+  double zprd = domain->zprd;
+  double ione[3][3];
+  for (i = 0; i < 3; i++)
+    for (j = 0; j < 3; j++)
+      ione[i][j] = 0.0;
+
+  for (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];
+      if (rmass) massone = rmass[i];
+      else massone = mass[type[i]];
+      ione[0][0] += massone * (dy*dy + dz*dz);
+      ione[1][1] += massone * (dx*dx + dz*dz);
+      ione[2][2] += massone * (dx*dx + dy*dy);
+      ione[0][1] -= massone * dx*dy;
+      ione[1][2] -= massone * dy*dz;
+      ione[0][2] -= massone * dx*dz;
+    }
+  ione[1][0] = ione[0][1];
+  ione[2][1] = ione[1][2];
+  ione[2][0] = ione[0][2];
+
+  MPI_Allreduce(&ione[0][0],&itensor[0][0],9,MPI_DOUBLE,MPI_SUM,world);
+}
+
 /* ----------------------------------------------------------------------
    compute angular velocity omega from L = Iw, inverting I to solve for w
    really not a group operation, but L and I were computed for a group

src/group.h

@@ -55,6 +55,7 @@ class Group : protected Pointers {
   void angmom(int, double *, double *);    // angular momentum of group
   void angmom(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
 
  private: