diff --git a/src/fix_nve_sphere.cpp b/src/fix_nve_sphere.cpp
index 9f7b4a9eaa175faf5644a9a8a2b283f7127e336f..42cca31f6b3c9038ff96cc0bdaf3c1891ce83e64 100644
--- a/src/fix_nve_sphere.cpp
+++ b/src/fix_nve_sphere.cpp
@@ -133,12 +133,13 @@ void FixNVESphere::initial_integrate(int vflag)
   
   // update mu for dipoles
   
-
   if (extra == DIPOLE) {
     double **mu = atom->mu;
     if (dlm == NODLM) {
+
       // d_mu/dt = omega cross mu
       // renormalize mu to dipole length
+
       for (int i = 0; i < nlocal; i++)
         if (mask[i] & groupbit)
           if (mu[i][3] > 0.0) {
@@ -152,7 +153,9 @@ void FixNVESphere::initial_integrate(int vflag)
             mu[i][2] = g[2]*scale;
           }
     } else {
-      // Integrate orientation following Dullweber-Leimkuhler-Maclachlan scheme
+
+      // integrate orientation following Dullweber-Leimkuhler-Maclachlan scheme
+
       for (int i = 0; i < nlocal; i++) {
         if (mask[i] & groupbit && mu[i][3] > 0.0) {
           
@@ -160,8 +163,9 @@ void FixNVESphere::initial_integrate(int vflag)
           // Q is the rotation matrix from space frame to body frame
           // i.e. v_b = Q.v_s
           
-          // Define mu to lie along the z axis in the body frame
-          // We take the unit dipole to avoid getting a scaling matrix
+          // define mu to lie along the z axis in the body frame
+          // take the unit dipole to avoid getting a scaling matrix
+
           inv_len_mu = 1.0/mu[i][3];
           a[0] = mu[i][0]*inv_len_mu;
           a[1] = mu[i][1]*inv_len_mu;
@@ -180,9 +184,15 @@ void FixNVESphere::initial_integrate(int vflag)
           if (s2 != 0.0){ // i.e. the vectors are not parallel
             scale = (1.0 - a[2])/s2;
             
-            Q[0][0] = 1.0 - scale*a[0]*a[0]; Q[0][1] = -scale*a[0]*a[1];      Q[0][2] = -a[0];
-            Q[1][0] = -scale*a[0]*a[1];      Q[1][1] = 1.0 - scale*a[1]*a[1]; Q[1][2] = -a[1];
-            Q[2][0] = a[0];                  Q[2][1] = a[1];                  Q[2][2] = 1.0 - scale*(a[0]*a[0] + a[1]*a[1]);
+            Q[0][0] = 1.0 - scale*a[0]*a[0]; 
+            Q[0][1] = -scale*a[0]*a[1];
+            Q[0][2] = -a[0];
+            Q[1][0] = -scale*a[0]*a[1];
+            Q[1][1] = 1.0 - scale*a[1]*a[1];
+            Q[1][2] = -a[1];
+            Q[2][0] = a[0];
+            Q[2][1] = a[1];
+            Q[2][2] = 1.0 - scale*(a[0]*a[0] + a[1]*a[1]);
           } else { // if parallel then we just have I or -I
             Q[0][0] = 1.0/a[2];  Q[0][1] = 0.0;       Q[0][2] = 0.0;
             Q[1][0] = 0.0;       Q[1][1] = 1.0/a[2];  Q[1][2] = 0.0;
@@ -242,7 +252,9 @@ void FixNVESphere::initial_integrate(int vflag)
           
           // Transform w back into space frame w_temp = Q^T.w
           transpose_matvec(Q_temp,w,w_temp);
-          omega[i][0] = w_temp[0]; omega[i][1] = w_temp[1]; omega[i][2] = w_temp[2];
+          omega[i][0] = w_temp[0]; 
+          omega[i][1] = w_temp[1];
+          omega[i][2] = w_temp[2];
           
           // Set dipole according to updated Q: mu = Q^T.[0 0 1] * |mu|
           mu[i][0] = Q_temp[2][0] * mu[i][3];
@@ -289,7 +301,8 @@ void FixNVESphere::final_integrate()
       omega[i][0] += dtirotate * torque[i][0];
       omega[i][1] += dtirotate * torque[i][1];
       omega[i][2] += dtirotate * torque[i][2];
-      rke += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + omega[i][2]*omega[i][2])*radius[i]*radius[i]*rmass[i];
+      rke += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + 
+              omega[i][2]*omega[i][2])*radius[i]*radius[i]*rmass[i];
     }
   
 }
diff --git a/src/group.cpp b/src/group.cpp
index 1322afcebf76c59891cc474405337ec164c37b0e..5ea815174fad476feccd7630691c773fe782d75e 100644
--- a/src/group.cpp
+++ b/src/group.cpp
@@ -1684,6 +1684,7 @@ void Group::omega(double *angmom, double inertia[3][3], double *w)
 
   // non-singular I matrix
   // use L = Iw, inverting I to solve for w
+  // this should give exact zeroing of omega by velocity command
 
   if (determinant > EPSILON) {
 
@@ -1715,7 +1716,7 @@ void Group::omega(double *angmom, double inertia[3][3], double *w)
       inverse[2][2]*angmom[2];
 
   // handle (nearly) singular I matrix
-  // due to 2-atom group or linear molecule
+  // typically due to 2-atom group or linear molecule
   // use jacobi() and angmom_to_omega() to calculate valid omega
 
   } else {