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

src/group.cpp

@@ -1150,7 +1150,8 @@ double Group::gyration(int igroup, double masstotal, double *cm, int iregion)
 }
 
 /* ----------------------------------------------------------------------
-   compute the angular momentum L (lmom) of group around center-of-mass cm
+   compute the angular momentum L (lmom) of group
+   around center-of-mass cm
    must unwrap atoms to compute L correctly
 ------------------------------------------------------------------------- */
 
@@ -1193,6 +1194,52 @@ void Group::angmom(int igroup, double *cm, double *lmom)
   MPI_Allreduce(p,lmom,3,MPI_DOUBLE,MPI_SUM,world);
 }
 
+/* ----------------------------------------------------------------------
+   compute the angular momentum L (lmom) of group of atoms in region
+   around center-of-mass cm
+   must unwrap atoms to compute L correctly
+------------------------------------------------------------------------- */
+
+void Group::angmom(int igroup, double *cm, double *lmom, int iregion)
+{
+  int groupbit = bitmask[igroup];
+  Region *region = domain->regions[iregion];
+
+  double **x = atom->x;
+  double **v = atom->v;
+  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 p[3];
+  p[0] = p[1] = p[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];
+      if (rmass) massone = rmass[i];
+      else massone = mass[type[i]];
+      p[0] += massone * (dy*v[i][2] - dz*v[i][1]);
+      p[1] += massone * (dz*v[i][0] - dx*v[i][2]);
+      p[2] += massone * (dx*v[i][1] - dy*v[i][0]);
+    }
+
+  MPI_Allreduce(p,lmom,3,MPI_DOUBLE,MPI_SUM,world);
+}
+
 /* ----------------------------------------------------------------------
    compute moment of inertia tensor around center-of-mass cm
    must unwrap atoms to compute itensor correctly

src/group.h

@@ -53,6 +53,7 @@ class Group : protected Pointers {
   double gyration(int, double, double *);  // radius-of-gyration of group
   double gyration(int, double, double *, int);
   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 omega(double *, double [3][3], double *);  // angular velocity
 

src/variable.cpp

@@ -1587,10 +1587,10 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
    word = group function
    contents = str bewteen parentheses with one,two,three args
    return 0 if not a match, 1 if successfully processed
-   customize by adding a group function:
+   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)
+     bound(group,xmin),gyration(group),ke(group),angmom(group)
 ------------------------------------------------------------------------- */
 
 int Variable::group_function(char *word, char *contents, Tree **tree,
@@ -1732,6 +1732,24 @@ int Variable::group_function(char *word, char *contents, Tree **tree,
     else if (narg == 2) value = group->ke(igroup,region_function(arg2));
     else error->all("Invalid group function in variable formula");
 
+  } else if (strcmp(word,"angmom") == 0) {
+    atom->check_mass();
+    double xcm[3],lmom[3];
+    if (narg == 2) {
+      double masstotal = group->mass(igroup);
+      group->xcm(igroup,masstotal,xcm);
+      group->angmom(igroup,xcm,lmom);
+    } else if (narg == 3) {
+      int iregion = region_function(arg3);
+      double masstotal = group->mass(igroup,iregion);
+      group->xcm(igroup,masstotal,xcm,iregion);
+      group->angmom(igroup,xcm,lmom,iregion);
+    } else error->all("Invalid group function in variable formula");
+    if (strcmp(arg2,"x") == 0) value = lmom[0];
+    else if (strcmp(arg2,"y") == 0) value = lmom[1];
+    else if (strcmp(arg2,"z") == 0) value = lmom[2];
+    else error->all("Invalid group function in variable formula");
+
   } else return 0;
     
   delete [] arg1;
@@ -1743,17 +1761,7 @@ int Variable::group_function(char *word, char *contents, Tree **tree,
   return 1;
 }
 
-/* ----------------------------------------------------------------------
-   process a group function in formula with optional region arg
-   push result onto tree or arg stack
-   word = group function
-   contents = str bewteen parentheses with one,two,three args
-   return 0 if not a match, 1 if successfully processed
-   customize by adding a group function:
-     count(group),mass(group),charge(group),
-     xcm(group,dim),vcm(group,dim),fcm(group,dim),
-     bound(group,xmin),gyration(group)
-------------------------------------------------------------------------- */
+/* ---------------------------------------------------------------------- */
 
 int Variable::region_function(char *id)
 {