--- /home/sjplimp/oldlammps/src/fix_rigid.cpp	2005-01-18 13:29:05.000000000 -0700
+++ /home/sjplimp/lammps/src/fix_rigid.cpp	2005-03-09 16:22:53.000000000 -0700
@@ -28,7 +28,7 @@
 #include "error.h"
 
 #define TOLERANCE 1.0e-6
-#define EPSILON 1.0e-9
+#define EPSILON 1.0e-7
 #define MAXJACOBI 50
 
 #define MIN(A,B) ((A) < (B)) ? (A) : (B)
@@ -40,7 +40,8 @@
 {
   int i,ibody;
 
-  // can't use with granular since mass arrays are different
+  // can't use with pure granular style since mass arrays are different
+  // hybrid granular style would be OK if fix were on non-granular atoms
 
   if (strcmp(atom->style,"granular") == 0)
     error->all("Cannot use fix rigid with atom style granular");
@@ -382,7 +383,7 @@
   
     ierror = jacobi(tensor,inertia[ibody],evectors);
     if (ierror) error->all("Insufficient Jacobi rotations for rigid body");
-  
+
     ex_space[ibody][0] = evectors[0][0];
     ex_space[ibody][1] = evectors[1][0];
     ex_space[ibody][2] = evectors[2][0];
@@ -395,16 +396,16 @@
     ez_space[ibody][1] = evectors[1][2];
     ez_space[ibody][2] = evectors[2][2];
     
-    // if any principal moment < scaled EPSILON, set to EPSILON
+    // if any principal moment < scaled EPSILON, set to 0.0
   
     double max;
     max = MAX(inertia[ibody][0],inertia[ibody][1]);
     max = MAX(max,inertia[ibody][2]);
   
-    inertia[ibody][0] = MAX(EPSILON*max,inertia[ibody][0]);
-    inertia[ibody][1] = MAX(EPSILON*max,inertia[ibody][1]);
-    inertia[ibody][2] = MAX(EPSILON*max,inertia[ibody][2]);
-  
+    if (inertia[ibody][0] < EPSILON*max) inertia[ibody][0] = 0.0;
+    if (inertia[ibody][1] < EPSILON*max) inertia[ibody][1] = 0.0;
+    if (inertia[ibody][2] < EPSILON*max) inertia[ibody][2] = 0.0;
+
     // enforce 3 evectors as a right-handed coordinate system
     // flip 3rd evector if needed
   
@@ -421,7 +422,7 @@
       ez_space[ibody][1] = -ez_space[ibody][1];
       ez_space[ibody][2] = -ez_space[ibody][2];
     }
-  
+
     // create initial quaternion
   
     qcreate(evectors,quat[ibody]);
@@ -598,12 +599,8 @@
 
 void FixRigid::initial_integrate()
 {
-  double wq[4],ex[4],ey[4],ez[4],qfull[4],qhalf[4];
   double dtfm;
 
-  fcm[0][0] = fcm[0][1] = fcm[0][2] = 0.0;
-  torque[0][0] = torque[0][1] = torque[0][2] = 0.0;
-
   for (int ibody = 0; ibody < nbody; ibody++) {
 
     // update vcm by 1/2 step
@@ -1000,6 +997,8 @@
    compute omega from angular momentum
    w = omega = angular velocity in space frame
    wbody = angular velocity in body frame
+   set wbody component to 0.0 if inertia component is 0.0
+     otherwise body can spin easily around that axis
    project space-frame angular momentum onto body axes
      and divide by principal moments
 ------------------------------------------------------------------------- */
@@ -1009,9 +1008,12 @@
 {
   double wbody[3];
 
-  wbody[0] = (m[0]*ex[0] + m[1]*ex[1] + m[2]*ex[2]) / inertia[0];
-  wbody[1] = (m[0]*ey[0] + m[1]*ey[1] + m[2]*ey[2]) / inertia[1];
-  wbody[2] = (m[0]*ez[0] + m[1]*ez[1] + m[2]*ez[2]) / inertia[2];
+  if (inertia[0] == 0.0) wbody[0] = 0.0;
+  else wbody[0] = (m[0]*ex[0] + m[1]*ex[1] + m[2]*ex[2]) / inertia[0];
+  if (inertia[1] == 0.0) wbody[1] = 0.0;
+  else wbody[1] = (m[0]*ey[0] + m[1]*ey[1] + m[2]*ey[2]) / inertia[1];
+  if (inertia[2] == 0.0) wbody[2] = 0.0;
+  else wbody[2] = (m[0]*ez[0] + m[1]*ez[1] + m[2]*ez[2]) / inertia[2];
 
   w[0] = wbody[0]*ex[0] + wbody[1]*ey[0] + wbody[2]*ez[0];
   w[1] = wbody[0]*ex[1] + wbody[1]*ey[1] + wbody[2]*ez[1];