diff --git a/src/ASPHERE/fix_nh_asphere.cpp b/src/ASPHERE/fix_nh_asphere.cpp
index 19c18e402c..fb50a7b2ff 100644
--- a/src/ASPHERE/fix_nh_asphere.cpp
+++ b/src/ASPHERE/fix_nh_asphere.cpp
@@ -71,14 +71,14 @@ void FixNHAsphere::richardson(double *q, double *m, double *moments)
   // full update from dq/dt = 1/2 w q
 
   double wq[4];
-  MathExtra::multiply_vec_quat(w,q,wq);
+  MathExtra::vecquat(w,q,wq);
 
   double qfull[4];
   qfull[0] = q[0] + dtq * wq[0];
   qfull[1] = q[1] + dtq * wq[1];
   qfull[2] = q[2] + dtq * wq[2];
   qfull[3] = q[3] + dtq * wq[3];
-  MathExtra::normalize4(qfull);
+  MathExtra::qnormalize(qfull);
 
   // 1st half of update from dq/dt = 1/2 w q
 
@@ -87,13 +87,13 @@ void FixNHAsphere::richardson(double *q, double *m, double *moments)
   qhalf[1] = q[1] + 0.5*dtq * wq[1];
   qhalf[2] = q[2] + 0.5*dtq * wq[2];
   qhalf[3] = q[3] + 0.5*dtq * wq[3];
-  MathExtra::normalize4(qhalf);
+  MathExtra::qnormalize(qhalf);
 
   // re-compute omega at 1/2 step from m at 1/2 step and q at 1/2 step
   // recompute wq
 
   omega_from_mq(qhalf,m,moments,w);
-  MathExtra::multiply_vec_quat(w,qhalf,wq);
+  MathExtra::vecquat(w,qhalf,wq);
 
   // 2nd half of update from dq/dt = 1/2 w q
 
@@ -101,7 +101,7 @@ void FixNHAsphere::richardson(double *q, double *m, double *moments)
   qhalf[1] += 0.5*dtq * wq[1];
   qhalf[2] += 0.5*dtq * wq[2];
   qhalf[3] += 0.5*dtq * wq[3];
-  MathExtra::normalize4(qhalf);
+  MathExtra::qnormalize(qhalf);
 
   // corrected Richardson update
 
@@ -109,7 +109,7 @@ void FixNHAsphere::richardson(double *q, double *m, double *moments)
   q[1] = 2.0*qhalf[1] - qfull[1];
   q[2] = 2.0*qhalf[2] - qfull[2];
   q[3] = 2.0*qhalf[3] - qfull[3];
-  MathExtra::normalize4(q);
+  MathExtra::qnormalize(q);
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/ASPHERE/fix_nve_asphere.cpp b/src/ASPHERE/fix_nve_asphere.cpp
index b444cbf114..9f1769b7f0 100755
--- a/src/ASPHERE/fix_nve_asphere.cpp
+++ b/src/ASPHERE/fix_nve_asphere.cpp
@@ -155,14 +155,14 @@ void FixNVEAsphere::richardson(double *q, double *m, double *moments)
   // full update from dq/dt = 1/2 w q
 
   double wq[4];
-  MathExtra::multiply_vec_quat(w,q,wq);
+  MathExtra::vecquat(w,q,wq);
 
   double qfull[4];
   qfull[0] = q[0] + dtq * wq[0];
   qfull[1] = q[1] + dtq * wq[1];
   qfull[2] = q[2] + dtq * wq[2];
   qfull[3] = q[3] + dtq * wq[3];
-  MathExtra::normalize4(qfull);
+  MathExtra::qnormalize(qfull);
 
   // 1st half of update from dq/dt = 1/2 w q
 
@@ -171,13 +171,13 @@ void FixNVEAsphere::richardson(double *q, double *m, double *moments)
   qhalf[1] = q[1] + 0.5*dtq * wq[1];
   qhalf[2] = q[2] + 0.5*dtq * wq[2];
   qhalf[3] = q[3] + 0.5*dtq * wq[3];
-  MathExtra::normalize4(qhalf);
+  MathExtra::qnormalize(qhalf);
 
   // re-compute omega at 1/2 step from m at 1/2 step and q at 1/2 step
   // recompute wq
 
   omega_from_mq(qhalf,m,moments,w);
-  MathExtra::multiply_vec_quat(w,qhalf,wq);
+  MathExtra::vecquat(w,qhalf,wq);
 
   // 2nd half of update from dq/dt = 1/2 w q
 
@@ -185,7 +185,7 @@ void FixNVEAsphere::richardson(double *q, double *m, double *moments)
   qhalf[1] += 0.5*dtq * wq[1];
   qhalf[2] += 0.5*dtq * wq[2];
   qhalf[3] += 0.5*dtq * wq[3];
-  MathExtra::normalize4(qhalf);
+  MathExtra::qnormalize(qhalf);
 
   // corrected Richardson update
 
@@ -193,7 +193,7 @@ void FixNVEAsphere::richardson(double *q, double *m, double *moments)
   q[1] = 2.0*qhalf[1] - qfull[1];
   q[2] = 2.0*qhalf[2] - qfull[2];
   q[3] = 2.0*qhalf[3] - qfull[3];
-  MathExtra::normalize4(q);
+  MathExtra::qnormalize(q);
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/ASPHERE/pair_gayberne.cpp b/src/ASPHERE/pair_gayberne.cpp
index 3cbcbac1f4..8517042416 100755
--- a/src/ASPHERE/pair_gayberne.cpp
+++ b/src/ASPHERE/pair_gayberne.cpp
@@ -43,8 +43,7 @@ enum{SPHERE_SPHERE,SPHERE_ELLIPSE,ELLIPSE_SPHERE,ELLIPSE_ELLIPSE};
 PairGayBerne::PairGayBerne(LAMMPS *lmp) : Pair(lmp)
 {
   avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
-  if (!avec) 
-    error->all("Pair gayberne requires atom style ellipsoid");
+  if (!avec) error->all("Pair gayberne requires atom style ellipsoid");
 
   single_enable = 0;
 }
diff --git a/src/ASPHERE/pair_resquared.cpp b/src/ASPHERE/pair_resquared.cpp
index 9b779e2e90..65f4b2bf01 100755
--- a/src/ASPHERE/pair_resquared.cpp
+++ b/src/ASPHERE/pair_resquared.cpp
@@ -45,8 +45,7 @@ PairRESquared::PairRESquared(LAMMPS *lmp) : Pair(lmp),
                                             cr60(pow(60.0,1.0/3.0))
 {
   avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
-  if (!avec) 
-    error->all("Pair gayberne requires atom style ellipsoid");
+  if (!avec) error->all("Pair resquared requires atom style ellipsoid");
 
   single_enable = 0;
 
diff --git a/src/SRD/fix_srd.cpp b/src/SRD/fix_srd.cpp
index b589908792..f315bc2ffa 100644
--- a/src/SRD/fix_srd.cpp
+++ b/src/SRD/fix_srd.cpp
@@ -2092,7 +2092,8 @@ void FixSRD::parameterize()
   // big particle must either have radius > 0 or shape > 0 defined
   // apply radfactor at end
 
-  AtomVecEllipsoid::Bonus *ebonus = avec_ellipsoid->bonus;
+  AtomVecEllipsoid::Bonus *ebonus;
+  if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
   double *radius = atom->radius;
   int *ellipsoid = atom->ellipsoid;
   int *mask = atom->mask;
@@ -2406,7 +2407,8 @@ void FixSRD::big_static()
   double rad,arad,brad,crad;
   double *shape;
 
-  AtomVecEllipsoid::Bonus *ebonus = avec_ellipsoid->bonus;
+  AtomVecEllipsoid::Bonus *ebonus;
+  if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
   double *radius = atom->radius;
   int *ellipsoid = atom->ellipsoid;
   int *type = atom->type;
@@ -2449,7 +2451,8 @@ void FixSRD::big_dynamic()
   int i;
   double *shape,*quat;
 
-  AtomVecEllipsoid::Bonus *ebonus = avec_ellipsoid->bonus;
+  AtomVecEllipsoid::Bonus *ebonus;
+  if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
   double **omega = atom->omega;
   double **angmom = atom->angmom;
   double *rmass = atom->rmass;
diff --git a/src/atom_vec_ellipsoid.cpp b/src/atom_vec_ellipsoid.cpp
index 5e00a4f24b..8a3ab90aaf 100755
--- a/src/atom_vec_ellipsoid.cpp
+++ b/src/atom_vec_ellipsoid.cpp
@@ -1198,7 +1198,7 @@ void AtomVecEllipsoid::data_atom_bonus(int m, char **values)
   quat[1] = atof(values[4]);
   quat[2] = atof(values[5]);
   quat[3] = atof(values[6]);
-  MathExtra::normalize4(quat);
+  MathExtra::qnormalize(quat);
 
   // reset ellipsoid mass
   // previously stored density in rmass
diff --git a/src/fix_rigid.cpp b/src/fix_rigid.cpp
index 647162b7eb..989196cea5 100644
--- a/src/fix_rigid.cpp
+++ b/src/fix_rigid.cpp
@@ -34,7 +34,6 @@ using namespace LAMMPS_NS;
 
 #define TOLERANCE 1.0e-6
 #define EPSILON 1.0e-7
-#define MAXJACOBI 50
 
 #define MIN(A,B) ((A) < (B)) ? (A) : (B)
 #define MAX(A,B) ((A) > (B)) ? (A) : (B)
@@ -470,7 +469,8 @@ void FixRigid::init()
 
   extended = dorientflag = qorientflag = 0;
 
-  AtomVecEllipsoid::Bonus *ebonus = avec_ellipsoid->bonus;
+  AtomVecEllipsoid::Bonus *ebonus;
+  if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
   double **mu = atom->mu;
   double *radius = atom->radius;
   double *rmass = atom->rmass;
@@ -589,6 +589,7 @@ void FixRigid::init()
   pre_neighbor();
 
   // compute 6 moments of inertia of each body
+  // stored as 6-vector in Voigt ordering
   // dx,dy,dz = coords relative to center-of-mass
   
   double dx,dy,dz,rad;
@@ -632,8 +633,7 @@ void FixRigid::init()
   // extended particles may contribute extra terms to moments of inertia
 
   if (extended) {
-    double ex[3],ey[3],ez[3],idiag[3];
-    double p[3][3],ptrans[3][3],ispace[3][3],itemp[3][3];
+    double ivec[6];
     double *shape,*quatatom;
 
     for (i = 0; i < nlocal; i++) {
@@ -652,19 +652,13 @@ void FixRigid::init()
       if (eflags[i] & INERTIA_ELLIPSOID) {
 	shape = ebonus[ellipsoid[i]].shape;
 	quatatom = ebonus[ellipsoid[i]].quat;
-	MathExtra::quat_to_mat(quatatom,p);
-	MathExtra::quat_to_mat_trans(quatatom,ptrans);
-	idiag[0] = 0.2*massone * (shape[1]*shape[1]+shape[2]*shape[2]);
-	idiag[1] = 0.2*massone * (shape[0]*shape[0]+shape[2]*shape[2]);
-	idiag[2] = 0.2*massone * (shape[0]*shape[0]+shape[1]*shape[1]);
-	MathExtra::diag_times3(idiag,ptrans,itemp);
-	MathExtra::times3(p,itemp,ispace);
-	sum[ibody][0] += ispace[0][0];
-	sum[ibody][1] += ispace[1][1];
-	sum[ibody][2] += ispace[2][2];
-	sum[ibody][3] += ispace[0][1];
-	sum[ibody][4] += ispace[1][2];
-	sum[ibody][5] += ispace[0][2];
+	MathExtra::inertia_ellipsoid(shape,quatatom,massone,ivec);
+	sum[ibody][0] += ivec[0];
+	sum[ibody][1] += ivec[1];
+	sum[ibody][2] += ivec[2];
+	sum[ibody][3] += ivec[3];
+	sum[ibody][4] += ivec[4];
+	sum[ibody][5] += ivec[5];
       }
     }
   }
@@ -674,12 +668,9 @@ void FixRigid::init()
   // inertia = 3 eigenvalues = principal moments of inertia
   // ex_space,ey_space,ez_space = 3 eigenvectors = principal axes of rigid body
   
-  double **tensor,**evectors;
-  memory->create(tensor,3,3,"fix_rigid:tensor");
-  memory->create(evectors,3,3,"fix_rigid:evectors");
-
-  int ierror;
+  double tensor[3][3],evectors[3][3];
   double ez0,ez1,ez2;
+  int ierror;
 
   for (ibody = 0; ibody < nbody; ibody++) {
     tensor[0][0] = all[ibody][0];
@@ -689,7 +680,7 @@ void FixRigid::init()
     tensor[1][2] = tensor[2][1] = all[ibody][4];
     tensor[0][2] = tensor[2][0] = all[ibody][5];
 
-    ierror = jacobi(tensor,inertia[ibody],evectors);
+    ierror = MathExtra::jacobi(tensor,inertia[ibody],evectors);
     if (ierror) error->all("Insufficient Jacobi rotations for rigid body");
 
     ex_space[ibody][0] = evectors[0][0];
@@ -733,14 +724,10 @@ void FixRigid::init()
 
     // create initial quaternion
   
-    q_from_exyz(ex_space[ibody],ey_space[ibody],ez_space[ibody],quat[ibody]);
+    MathExtra::exyz_to_q(ex_space[ibody],ey_space[ibody],ez_space[ibody],
+			 quat[ibody]);
   }
 
-  // free temporary memory
-  
-  memory->destroy(tensor);
-  memory->destroy(evectors);
-  
   // displace = initial atom coords in basis of principal axes
   // set displace = 0.0 for atoms not in any rigid body
   // for extended particles, set their orientation wrt to rigid body
@@ -795,9 +782,9 @@ void FixRigid::init()
 
       if (eflags[i] & ORIENT_QUAT) {
 	quatatom = ebonus[ellipsoid[i]].quat;
-	qconjugate(quat[ibody],qc);
-	quatquat(qc,quatatom,qorient[i]);
-	qnormalize(qorient[i]);
+	MathExtra::qconjugate(quat[ibody],qc);
+	MathExtra::quatquat(qc,quatatom,qorient[i]);
+	MathExtra::qnormalize(qorient[i]);
       } else if (qorientflag)
 	qorient[i][0] = qorient[i][1] = qorient[i][2] = qorient[i][3] = 0.0;
     }
@@ -805,10 +792,11 @@ void FixRigid::init()
   
   // test for valid principal moments & axes
   // recompute moments of inertia around new axes
+  // stored as 6-vector in Voigt ordering
   // 3 diagonal moments should equal principal moments
   // 3 off-diagonal moments should be 0.0
   // extended particles may contribute extra terms to moments of inertia
-  
+
   for (ibody = 0; ibody < nbody; ibody++)
     for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
   
@@ -830,8 +818,7 @@ void FixRigid::init()
   }
 
   if (extended) {
-    double ex[3],ey[3],ez[3],idiag[3];
-    double p[3][3],ptrans[3][3],ispace[3][3],itemp[3][3];
+    double ivec[6];
     double *shape;
 
     for (i = 0; i < nlocal; i++) {
@@ -849,19 +836,13 @@ void FixRigid::init()
       }
       if (eflags[i] & INERTIA_ELLIPSOID) {
 	shape = ebonus[ellipsoid[i]].shape;
-	MathExtra::quat_to_mat(qorient[i],p);
-	MathExtra::quat_to_mat_trans(qorient[i],ptrans);
-	idiag[0] = 0.2*massone * (shape[1]*shape[1]+shape[2]*shape[2]);
-	idiag[1] = 0.2*massone * (shape[0]*shape[0]+shape[2]*shape[2]);
-	idiag[2] = 0.2*massone * (shape[0]*shape[0]+shape[1]*shape[1]);
-	MathExtra::diag_times3(idiag,ptrans,itemp);
-	MathExtra::times3(p,itemp,ispace);
-	sum[ibody][0] += ispace[0][0];
-	sum[ibody][1] += ispace[1][1];
-	sum[ibody][2] += ispace[2][2];
-	sum[ibody][3] += ispace[0][1];
-	sum[ibody][4] += ispace[1][2];
-	sum[ibody][5] += ispace[0][2];
+	MathExtra::inertia_ellipsoid(shape,qorient[i],massone,ivec);
+	sum[ibody][0] += ivec[0];
+	sum[ibody][1] += ivec[1];
+	sum[ibody][2] += ivec[2];
+	sum[ibody][3] += ivec[3];
+	sum[ibody][4] += ivec[4];
+	sum[ibody][5] += ivec[5];
       }
     }
   }
@@ -1048,8 +1029,8 @@ void FixRigid::setup(int vflag)
   // set velocities from angmom & omega
 
   for (ibody = 0; ibody < nbody; ibody++)
-    omega_from_angmom(angmom[ibody],ex_space[ibody],ey_space[ibody],
-		      ez_space[ibody],inertia[ibody],omega[ibody]);
+    MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody],
+			       ez_space[ibody],inertia[ibody],omega[ibody]);
   set_v();
 
   // guesstimate virial as 2x the set_v contribution
@@ -1209,8 +1190,8 @@ void FixRigid::final_integrate()
     angmom[ibody][1] += dtf * torque[ibody][1] * tflag[ibody][1];
     angmom[ibody][2] += dtf * torque[ibody][2] * tflag[ibody][2];
 
-    omega_from_angmom(angmom[ibody],ex_space[ibody],ey_space[ibody],
-		      ez_space[ibody],inertia[ibody],omega[ibody]);
+    MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody],
+			       ez_space[ibody],inertia[ibody],omega[ibody]);
   }
 
   // set velocity/rotation of atoms in rigid bodies
@@ -1227,20 +1208,19 @@ void FixRigid::richardson(double *q, double *w,
 {
   // compute omega at 1/2 step from m at 1/2 step and q at 0
   
-  omega_from_angmom(m,ex,ey,ez,moments,w);
+  MathExtra::angmom_to_omega(m,ex,ey,ez,moments,w);
 
   // full update from dq/dt = 1/2 w q
 
   double wq[4];
-  vecquat(w,q,wq);
+  MathExtra::vecquat(w,q,wq);
 
   double qfull[4];
   qfull[0] = q[0] + dtq * wq[0];
   qfull[1] = q[1] + dtq * wq[1];
   qfull[2] = q[2] + dtq * wq[2];
   qfull[3] = q[3] + dtq * wq[3];
-
-  qnormalize(qfull);
+  MathExtra::qnormalize(qfull);
 
   // 1st half update from dq/dt = 1/2 w q
 
@@ -1249,16 +1229,15 @@ void FixRigid::richardson(double *q, double *w,
   qhalf[1] = q[1] + 0.5*dtq * wq[1];
   qhalf[2] = q[2] + 0.5*dtq * wq[2];
   qhalf[3] = q[3] + 0.5*dtq * wq[3];
-
-  qnormalize(qhalf);
+  MathExtra::qnormalize(qhalf);
 
   // udpate ex,ey,ez from qhalf
   // re-compute omega at 1/2 step from m at 1/2 step and q at 1/2 step
   // recompute wq
 
-  exyz_from_q(qhalf,ex,ey,ez);
-  omega_from_angmom(m,ex,ey,ez,moments,w);
-  vecquat(w,qhalf,wq);
+  MathExtra::q_to_exyz(qhalf,ex,ey,ez);
+  MathExtra::angmom_to_omega(m,ex,ey,ez,moments,w);
+  MathExtra::vecquat(w,qhalf,wq);
 
   // 2nd half update from dq/dt = 1/2 w q
 
@@ -1266,8 +1245,7 @@ void FixRigid::richardson(double *q, double *w,
   qhalf[1] += 0.5*dtq * wq[1];
   qhalf[2] += 0.5*dtq * wq[2];
   qhalf[3] += 0.5*dtq * wq[3];
-
-  qnormalize(qhalf);
+  MathExtra::qnormalize(qhalf);
 
   // corrected Richardson update
 
@@ -1275,9 +1253,9 @@ void FixRigid::richardson(double *q, double *w,
   q[1] = 2.0*qhalf[1] - qfull[1];
   q[2] = 2.0*qhalf[2] - qfull[2];
   q[3] = 2.0*qhalf[3] - qfull[3];
+  MathExtra::qnormalize(q);
 
-  qnormalize(q);
-  exyz_from_q(q,ex,ey,ez);
+  MathExtra::q_to_exyz(q,ex,ey,ez);
 }
 
 /* ----------------------------------------------------------------------
@@ -1511,304 +1489,6 @@ void FixRigid::deform(int flag)
       domain->lamda2x(xcm[ibody],xcm[ibody]);
 }
 
-/* ----------------------------------------------------------------------
-   compute evalues and evectors of 3x3 real symmetric matrix
-   based on Jacobi rotations
-   adapted from Numerical Recipes jacobi() function
-------------------------------------------------------------------------- */
-
-int FixRigid::jacobi(double **matrix, double *evalues, double **evectors)
-{
-  int i,j,k;
-  double tresh,theta,tau,t,sm,s,h,g,c,b[3],z[3];
-  
-  for (i = 0; i < 3; i++) {
-    for (j = 0; j < 3; j++) evectors[i][j] = 0.0;
-    evectors[i][i] = 1.0;
-  }
-  for (i = 0; i < 3; i++) {
-    b[i] = evalues[i] = matrix[i][i];
-    z[i] = 0.0;
-  }
-  
-  for (int iter = 1; iter <= MAXJACOBI; iter++) {
-    sm = 0.0;
-    for (i = 0; i < 2; i++)
-      for (j = i+1; j < 3; j++)
-	sm += fabs(matrix[i][j]);
-    if (sm == 0.0) return 0;
-    
-    if (iter < 4) tresh = 0.2*sm/(3*3);
-    else tresh = 0.0;
-    
-    for (i = 0; i < 2; i++) {
-      for (j = i+1; j < 3; j++) {
-	g = 100.0*fabs(matrix[i][j]);
-	if (iter > 4 && fabs(evalues[i])+g == fabs(evalues[i])
-	    && fabs(evalues[j])+g == fabs(evalues[j]))
-	  matrix[i][j] = 0.0;
-	else if (fabs(matrix[i][j]) > tresh) {
-	  h = evalues[j]-evalues[i];
-	  if (fabs(h)+g == fabs(h)) t = (matrix[i][j])/h;
-	  else {
-	    theta = 0.5*h/(matrix[i][j]);
-	    t = 1.0/(fabs(theta)+sqrt(1.0+theta*theta));
-	    if (theta < 0.0) t = -t;
-	  }
-	  c = 1.0/sqrt(1.0+t*t);
-	  s = t*c;
-	  tau = s/(1.0+c);
-	  h = t*matrix[i][j];
-	  z[i] -= h;
-	  z[j] += h;
-	  evalues[i] -= h;
-	  evalues[j] += h;
-	  matrix[i][j] = 0.0;
-	  for (k = 0; k < i; k++) rotate(matrix,k,i,k,j,s,tau);
-	  for (k = i+1; k < j; k++) rotate(matrix,i,k,k,j,s,tau);
-	  for (k = j+1; k < 3; k++) rotate(matrix,i,k,j,k,s,tau);
-	  for (k = 0; k < 3; k++) rotate(evectors,k,i,k,j,s,tau);
-	}
-      }
-    }
-    
-    for (i = 0; i < 3; i++) {
-      evalues[i] = b[i] += z[i];
-      z[i] = 0.0;
-    }
-  }
-  return 1;
-}
-
-/* ----------------------------------------------------------------------
-   perform a single Jacobi rotation
-------------------------------------------------------------------------- */
-
-void FixRigid::rotate(double **matrix, int i, int j, int k, int l,
-		      double s, double tau)
-{
-  double g = matrix[i][j];
-  double h = matrix[k][l];
-  matrix[i][j] = g-s*(h+g*tau);
-  matrix[k][l] = h+s*(g-h*tau);
-}
-
-/* ----------------------------------------------------------------------
-   create unit quaternion from space-frame ex,ey,ez
-   ex,ey,ez are columns of a rotation matrix
-------------------------------------------------------------------------- */
-
-void FixRigid::q_from_exyz(double *ex, double *ey, double *ez, double *q)
-{
-  // squares of quaternion components
-  
-  double q0sq = 0.25 * (ex[0] + ey[1] + ez[2] + 1.0);
-  double q1sq = q0sq - 0.5 * (ey[1] + ez[2]);
-  double q2sq = q0sq - 0.5 * (ex[0] + ez[2]);
-  double q3sq = q0sq - 0.5 * (ex[0] + ey[1]);
-  
-  // some component must be greater than 1/4 since they sum to 1
-  // compute other components from it
-  
-  if (q0sq >= 0.25) {
-    q[0] = sqrt(q0sq);
-    q[1] = (ey[2] - ez[1]) / (4.0*q[0]);
-    q[2] = (ez[0] - ex[2]) / (4.0*q[0]);
-    q[3] = (ex[1] - ey[0]) / (4.0*q[0]);
-  } else if (q1sq >= 0.25) {
-    q[1] = sqrt(q1sq);
-    q[0] = (ey[2] - ez[1]) / (4.0*q[1]);
-    q[2] = (ey[0] + ex[1]) / (4.0*q[1]);
-    q[3] = (ex[2] + ez[0]) / (4.0*q[1]);
-  } else if (q2sq >= 0.25) {
-    q[2] = sqrt(q2sq);
-    q[0] = (ez[0] - ex[2]) / (4.0*q[2]);
-    q[1] = (ey[0] + ex[1]) / (4.0*q[2]);
-    q[3] = (ez[1] + ey[2]) / (4.0*q[2]);
-  } else if (q3sq >= 0.25) {
-    q[3] = sqrt(q3sq);
-    q[0] = (ex[1] - ey[0]) / (4.0*q[3]);
-    q[1] = (ez[0] + ex[2]) / (4.0*q[3]);
-    q[2] = (ez[1] + ey[2]) / (4.0*q[3]);
-  } else
-    error->all("Quaternion creation numeric error");
-
-  qnormalize(q);
-}
-
-/* ----------------------------------------------------------------------
-   compute space-frame ex,ey,ez from current quaternion q
-   ex,ey,ez = space-frame coords of 1st,2nd,3rd principal axis
-   operation is ex = q' d q = Q d, where d is (1,0,0) = 1st axis in body frame
-------------------------------------------------------------------------- */
-
-void FixRigid::exyz_from_q(double *q, double *ex, double *ey, double *ez)
-{
-  ex[0] = q[0]*q[0] + q[1]*q[1] - q[2]*q[2] - q[3]*q[3];
-  ex[1] = 2.0 * (q[1]*q[2] + q[0]*q[3]);
-  ex[2] = 2.0 * (q[1]*q[3] - q[0]*q[2]);
-  
-  ey[0] = 2.0 * (q[1]*q[2] - q[0]*q[3]);
-  ey[1] = q[0]*q[0] - q[1]*q[1] + q[2]*q[2] - q[3]*q[3];
-  ey[2] = 2.0 * (q[2]*q[3] + q[0]*q[1]);
-  
-  ez[0] = 2.0 * (q[1]*q[3] + q[0]*q[2]);
-  ez[1] = 2.0 * (q[2]*q[3] - q[0]*q[1]);
-  ez[2] = q[0]*q[0] - q[1]*q[1] - q[2]*q[2] + q[3]*q[3];
-}
-
-/* ----------------------------------------------------------------------
-   vector-quaternion multiply: c = a*b, where a = (0,a)
-------------------------------------------------------------------------- */
-
-void FixRigid::vecquat(double *a, double *b, double *c)
-{
-  c[0] = -a[0]*b[1] - a[1]*b[2] - a[2]*b[3];
-  c[1] = b[0]*a[0] + a[1]*b[3] - a[2]*b[2];
-  c[2] = b[0]*a[1] + a[2]*b[1] - a[0]*b[3];
-  c[3] = b[0]*a[2] + a[0]*b[2] - a[1]*b[1];
-}
-
-/* ----------------------------------------------------------------------
-   quaternion-vector multiply: c = a*b, where b = (0,b)
-------------------------------------------------------------------------- */
-
-void FixRigid::quatvec(double *a, double *b, double *c)
-{
-  c[0] = -a[1]*b[0] - a[2]*b[1] - a[3]*b[2];
-  c[1] = a[0]*b[0] + a[2]*b[2] - a[3]*b[1];
-  c[2] = a[0]*b[1] + a[3]*b[0] - a[1]*b[2];
-  c[3] = a[0]*b[2] + a[1]*b[1] - a[2]*b[0];
-}
-
-/* ----------------------------------------------------------------------
-   quaternion-quaternion multiply: c = a*b
-------------------------------------------------------------------------- */
-
-void FixRigid::quatquat(double *a, double *b, double *c)
-{
-  c[0] = a[0]*b[0] - a[1]*b[1] - a[2]*b[2] - a[3]*b[3];
-  c[1] = a[0]*b[1] + b[0]*a[1] + a[2]*b[3] - a[3]*b[2];
-  c[2] = a[0]*b[2] + b[0]*a[2] + a[3]*b[1] - a[1]*b[3];
-  c[3] = a[0]*b[3] + b[0]*a[3] + a[1]*b[2] - a[2]*b[1];
-}
-
-/* ----------------------------------------------------------------------
-   quaternion multiply: c = inv(a)*b
-   a is a quaternion
-   b is a four component vector
-   c is a three component vector
-------------------------------------------------------------------------- */
-
-void FixRigid::invquatvec(double *a, double *b, double *c)
-{ 
-  c[0] = -a[1]*b[0] + a[0]*b[1] + a[3]*b[2] - a[2]*b[3];
-  c[1] = -a[2]*b[0] - a[3]*b[1] + a[0]*b[2] + a[1]*b[3];
-  c[2] = -a[3]*b[0] + a[2]*b[1] - a[1]*b[2] + a[0]*b[3];
-}
-
-/* ----------------------------------------------------------------------
-   conjugate of a quaternion: qc = conjugate of q
-   assume q is of unit length
-------------------------------------------------------------------------- */
-
-void FixRigid::qconjugate(double *q, double *qc)
-{
-  qc[0] = q[0];
-  qc[1] = -q[1];
-  qc[2] = -q[2];
-  qc[3] = -q[3];
-}
-
-/* ----------------------------------------------------------------------
-   normalize a quaternion
-------------------------------------------------------------------------- */
-
-void FixRigid::qnormalize(double *q)
-{
-  double norm = 1.0 / sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + q[3]*q[3]);
-  q[0] *= norm;
-  q[1] *= norm;
-  q[2] *= norm;
-  q[3] *= norm;
-}
-
-/* ----------------------------------------------------------------------
-  matvec_rows: c = Ab, where rows of A are x, y, z
-------------------------------------------------------------------------- */
-
-void FixRigid::matvec_rows(double *x, double *y, double *z, 
-			   double *b, double *c)
-{
-  c[0] = x[0]*b[0] + x[1]*b[1] + x[2]*b[2];
-  c[1] = y[0]*b[0] + y[1]*b[1] + y[2]*b[2];
-  c[2] = z[0]*b[0] + z[1]*b[1] + z[2]*b[2];
-}
-
-/* ----------------------------------------------------------------------
-  matvec_cols: c = Ab, where columns of A are x, y, z
-------------------------------------------------------------------------- */
-
-void FixRigid::matvec_cols(double *x, double *y, double *z,
-			   double *b, double *c)
-{
-  c[0] = x[0]*b[0] + y[0]*b[1] + z[0]*b[2];
-  c[1] = x[1]*b[0] + y[1]*b[1] + z[1]*b[2];
-  c[2] = x[2]*b[0] + y[2]*b[1] + z[2]*b[2];
-}
-
-/* ----------------------------------------------------------------------
-   compute omega from angular momentum, both in space frame
-   only know Idiag so need to do M = Iw in body frame
-   ex,ey,ez are column vectors of rotation matrix P
-   Mbody = P_transpose Mspace
-   wbody = Mbody / Idiag
-   wspace = P wbody
-   set wbody component to 0.0 if inertia component is 0.0
-     otherwise body can spin easily around that axis
-------------------------------------------------------------------------- */
-
-void FixRigid::omega_from_angmom(double *m, double *ex, double *ey, double *ez,
-				 double *idiag, double *w)
-{
-  double wbody[3];
-
-  if (idiag[0] == 0.0) wbody[0] = 0.0;
-  else wbody[0] = (m[0]*ex[0] + m[1]*ex[1] + m[2]*ex[2]) / idiag[0];
-  if (idiag[1] == 0.0) wbody[1] = 0.0;
-  else wbody[1] = (m[0]*ey[0] + m[1]*ey[1] + m[2]*ey[2]) / idiag[1];
-  if (idiag[2] == 0.0) wbody[2] = 0.0;
-  else wbody[2] = (m[0]*ez[0] + m[1]*ez[1] + m[2]*ez[2]) / idiag[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];
-  w[2] = wbody[0]*ex[2] + wbody[1]*ey[2] + wbody[2]*ez[2];
-}
-
-/* ----------------------------------------------------------------------
-   compute angular momentum from omega, both in space frame
-   only know Idiag so need to do M = Iw in body frame
-   ex,ey,ez are column vectors of rotation matrix P
-   wbody = P_transpose wspace
-   Mbody = Idiag wbody
-   Mspace = P Mbody
-------------------------------------------------------------------------- */
-
-void FixRigid::angmom_from_omega(double *w, 
-				  double *ex, double *ey, double *ez,
-				  double *idiag, double *m)
-{
-  double mbody[3];
-
-  mbody[0] = (w[0]*ex[0] + w[1]*ex[1] + w[2]*ex[2]) * idiag[0];
-  mbody[1] = (w[0]*ey[0] + w[1]*ey[1] + w[2]*ey[2]) * idiag[1];
-  mbody[2] = (w[0]*ez[0] + w[1]*ez[1] + w[2]*ez[2]) * idiag[2];
-
-  m[0] = mbody[0]*ex[0] + mbody[1]*ey[0] + mbody[2]*ez[0];
-  m[1] = mbody[0]*ex[1] + mbody[1]*ey[1] + mbody[2]*ez[1];
-  m[2] = mbody[0]*ex[2] + mbody[1]*ey[2] + mbody[2]*ez[2];
-}
-
 /* ----------------------------------------------------------------------
    set space-frame coords and velocity of each atom in each rigid body
    set orientation and rotation of extended particles
@@ -1933,7 +1613,8 @@ void FixRigid::set_xv()
   if (extended) {
     double *shape,*quatatom;
 
-    AtomVecEllipsoid::Bonus *ebonus = avec_ellipsoid->bonus;
+    AtomVecEllipsoid::Bonus *ebonus;
+    if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
     double **omega_one = atom->omega;
     double **angmom_one = atom->angmom;
     double **mu = atom->mu;
@@ -1950,8 +1631,8 @@ void FixRigid::set_xv()
       }
       if (eflags[i] & ORIENT_QUAT) {
 	quatatom = ebonus[ellipsoid[i]].quat;
-	quatquat(quat[ibody],qorient[i],quatatom);
-	qnormalize(quatatom);
+	MathExtra::quatquat(quat[ibody],qorient[i],quatatom);
+	MathExtra::qnormalize(quatatom);
       }
       if (eflags[i] & OMEGA) {
 	omega_one[i][0] = omega[ibody][0];
@@ -1964,8 +1645,9 @@ void FixRigid::set_xv()
 	ione[0] = 0.2*rmass[i] * (shape[1]*shape[1] + shape[2]*shape[2]);
 	ione[1] = 0.2*rmass[i] * (shape[0]*shape[0] + shape[2]*shape[2]);
 	ione[2] = 0.2*rmass[i] * (shape[0]*shape[0] + shape[1]*shape[1]);
-	exyz_from_q(quatatom,exone,eyone,ezone);
-	angmom_from_omega(omega[ibody],exone,eyone,ezone,ione,angmom_one[i]);
+	MathExtra::q_to_exyz(quatatom,exone,eyone,ezone);
+	MathExtra::omega_to_angmom(omega[ibody],exone,eyone,ezone,ione,
+				   angmom_one[i]);
       }
     }
   }
@@ -2075,7 +1757,8 @@ void FixRigid::set_v()
   if (extended) {
     double *shape,*quatatom;
 
-    AtomVecEllipsoid::Bonus *ebonus = avec_ellipsoid->bonus;
+    AtomVecEllipsoid::Bonus *ebonus;
+    if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
     double **omega_one = atom->omega;
     double **angmom_one = atom->angmom;
     int *ellipsoid = atom->ellipsoid;
@@ -2095,8 +1778,9 @@ void FixRigid::set_v()
 	ione[0] = 0.2*rmass[i] * (shape[1]*shape[1] + shape[2]*shape[2]);
 	ione[1] = 0.2*rmass[i] * (shape[0]*shape[0] + shape[2]*shape[2]);
 	ione[2] = 0.2*rmass[i] * (shape[0]*shape[0] + shape[1]*shape[1]);
-	exyz_from_q(quatatom,exone,eyone,ezone);
-	angmom_from_omega(omega[ibody],exone,eyone,ezone,ione,angmom_one[i]);
+	MathExtra::q_to_exyz(quatatom,exone,eyone,ezone);
+	MathExtra::omega_to_angmom(omega[ibody],exone,eyone,ezone,ione,
+				   angmom_one[i]);
       }
     }
   }
diff --git a/src/fix_rigid.h b/src/fix_rigid.h
index 6428d7c327..0af9175d54 100644
--- a/src/fix_rigid.h
+++ b/src/fix_rigid.h
@@ -102,27 +102,9 @@ class FixRigid : public Fix {
   int ORIENT_DIPOLE,ORIENT_QUAT;
   int OMEGA,ANGMOM,TORQUE;
 
-  int jacobi(double **, double *, double **);
-  void rotate(double **, int, int, int, int, double, double);
-  void q_from_exyz(double *, double *, double *, double *);
-  void exyz_from_q(double *, double *, double *, double *);
-
-  void vecquat(double *, double *, double *);
-  void quatvec(double *, double *, double *);
-  void quatquat(double *, double *, double *);
-  void invquatvec(double *, double *, double *); 
-  void qconjugate(double *, double *);
-  void qnormalize(double *);
-  void matvec_rows(double *, double *, double *, double *, double *);
-  void matvec_cols(double *, double *, double *, double *, double *);
-
   void richardson(double *, double *, double *, double *,
 		  double *, double *, double *);
   void no_squish_rotate(int, double *, double *, double *, double);
-  void omega_from_angmom(double *, double *, double *,
-			 double *, double *, double *);
-  void angmom_from_omega(double *, double *, double *,
-			 double *, double *, double *);
   void set_xv();
   void set_v();
 };
diff --git a/src/fix_rigid_nve.cpp b/src/fix_rigid_nve.cpp
index f99a20ff71..ccd908e8f8 100644
--- a/src/fix_rigid_nve.cpp
+++ b/src/fix_rigid_nve.cpp
@@ -21,6 +21,7 @@
 #include "stdio.h"
 #include "string.h"
 #include "fix_rigid_nve.h"
+#include "math_extra.h"
 #include "atom.h"
 #include "domain.h"
 #include "update.h"
@@ -57,9 +58,9 @@ void FixRigidNVE::setup(int vflag)
   
   double mbody[3];
   for (int ibody = 0; ibody < nbody; ibody++) {
-    matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody],
-		angmom[ibody],mbody);
-    quatvec(quat[ibody],mbody,conjqm[ibody]);
+    MathExtra::matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody],
+			   angmom[ibody],mbody);
+    MathExtra::quatvec(quat[ibody],mbody,conjqm[ibody]);
     conjqm[ibody][0] *= 2.0;
     conjqm[ibody][1] *= 2.0;
     conjqm[ibody][2] *= 2.0;
@@ -98,9 +99,9 @@ void FixRigidNVE::initial_integrate(int vflag)
     torque[ibody][1] *= tflag[ibody][1];
     torque[ibody][2] *= tflag[ibody][2];
     
-    matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody],
-		torque[ibody],tbody);
-    quatvec(quat[ibody],tbody,fquat);
+    MathExtra::matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody],
+			   torque[ibody],tbody);
+    MathExtra::quatvec(quat[ibody],tbody,fquat);
     
     conjqm[ibody][0] += dtf2 * fquat[0];
     conjqm[ibody][1] += dtf2 * fquat[1];
@@ -119,18 +120,18 @@ void FixRigidNVE::initial_integrate(int vflag)
     // transform p back to angmom
     // update angular velocity
     
-    exyz_from_q(quat[ibody],ex_space[ibody],ey_space[ibody],
-      ez_space[ibody]);
-    invquatvec(quat[ibody],conjqm[ibody],mbody);
-    matvec_cols(ex_space[ibody],ey_space[ibody],ez_space[ibody],
-		mbody,angmom[ibody]);
+    MathExtra::q_to_exyz(quat[ibody],ex_space[ibody],ey_space[ibody],
+			 ez_space[ibody]);
+    MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody);
+    MathExtra::matvec_cols(ex_space[ibody],ey_space[ibody],ez_space[ibody],
+			   mbody,angmom[ibody]);
     
     angmom[ibody][0] *= 0.5;
     angmom[ibody][1] *= 0.5;
     angmom[ibody][2] *= 0.5;
     
-    omega_from_angmom(angmom[ibody],ex_space[ibody],ey_space[ibody],
-      ez_space[ibody],inertia[ibody],omega[ibody]);
+    MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody],
+			       ez_space[ibody],inertia[ibody],omega[ibody]);
   }
   
   // virial setup before call to set_xv
@@ -247,25 +248,25 @@ void FixRigidNVE::final_integrate()
     torque[ibody][1] *= tflag[ibody][1];
     torque[ibody][2] *= tflag[ibody][2];
     
-    matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody],
-		torque[ibody],tbody);
-    quatvec(quat[ibody],tbody,fquat);
+    MathExtra::matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody],
+			   torque[ibody],tbody);
+    MathExtra::quatvec(quat[ibody],tbody,fquat);
     
     conjqm[ibody][0] += dtf2 * fquat[0];
     conjqm[ibody][1] += dtf2 * fquat[1];
     conjqm[ibody][2] += dtf2 * fquat[2];
     conjqm[ibody][3] += dtf2 * fquat[3];
     
-    invquatvec(quat[ibody],conjqm[ibody],mbody);
-    matvec_cols(ex_space[ibody],ey_space[ibody],ez_space[ibody],
-		mbody,angmom[ibody]);
-      
+    MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody);
+    MathExtra::matvec_cols(ex_space[ibody],ey_space[ibody],ez_space[ibody],
+			   mbody,angmom[ibody]);
+    
     angmom[ibody][0] *= 0.5;
     angmom[ibody][1] *= 0.5;
     angmom[ibody][2] *= 0.5;  
     
-    omega_from_angmom(angmom[ibody],ex_space[ibody],ey_space[ibody],
-		  ez_space[ibody],inertia[ibody],omega[ibody]);
+    MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody],
+			       ez_space[ibody],inertia[ibody],omega[ibody]);
   }
   
   // set velocity/rotation of atoms in rigid bodies
diff --git a/src/fix_rigid_nvt.cpp b/src/fix_rigid_nvt.cpp
index 34f83de9ee..816ea19a37 100644
--- a/src/fix_rigid_nvt.cpp
+++ b/src/fix_rigid_nvt.cpp
@@ -22,6 +22,7 @@
 #include "stdlib.h"
 #include "string.h"
 #include "fix_rigid_nvt.h"
+#include "math_extra.h"
 #include "atom.h"
 #include "domain.h"
 #include "update.h"
@@ -163,9 +164,9 @@ void FixRigidNVT::setup(int vflag)
   
   double mbody[3];
   for (int ibody = 0; ibody < nbody; ibody++) {
-    matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody],
-		angmom[ibody],mbody);
-    quatvec(quat[ibody],mbody,conjqm[ibody]);
+    MathExtra::matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody],
+			   angmom[ibody],mbody);
+    MathExtra::quatvec(quat[ibody],mbody,conjqm[ibody]);
     conjqm[ibody][0] *= 2.0;
     conjqm[ibody][1] *= 2.0;
     conjqm[ibody][2] *= 2.0;
@@ -224,9 +225,9 @@ void FixRigidNVT::initial_integrate(int vflag)
     
     // step 1.3 - apply torque (body coords) to quaternion momentum
 
-    matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody],
-		torque[ibody],tbody);
-    quatvec(quat[ibody],tbody,fquat);
+    MathExtra::matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody],
+			   torque[ibody],tbody);
+    MathExtra::quatvec(quat[ibody],tbody,fquat);
     
     conjqm[ibody][0] += dtf2 * fquat[0];
     conjqm[ibody][1] += dtf2 * fquat[1];
@@ -249,18 +250,18 @@ void FixRigidNVT::initial_integrate(int vflag)
     // transform p back to angmom
     // update angular velocity
     
-    exyz_from_q(quat[ibody],ex_space[ibody],ey_space[ibody],
-      ez_space[ibody]);
-    invquatvec(quat[ibody],conjqm[ibody],mbody);
-    matvec_cols(ex_space[ibody],ey_space[ibody],ez_space[ibody],
-		mbody,angmom[ibody]);
+    MathExtra::q_to_exyz(quat[ibody],ex_space[ibody],ey_space[ibody],
+			 ez_space[ibody]);
+    MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody);
+    MathExtra::matvec_cols(ex_space[ibody],ey_space[ibody],ez_space[ibody],
+			   mbody,angmom[ibody]);
     
     angmom[ibody][0] *= 0.5;
     angmom[ibody][1] *= 0.5;
     angmom[ibody][2] *= 0.5;
     
-    omega_from_angmom(angmom[ibody],ex_space[ibody],ey_space[ibody],
-      ez_space[ibody],inertia[ibody],omega[ibody]);
+    MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody],
+			       ez_space[ibody],inertia[ibody],omega[ibody]);
     
     akin_r += angmom[ibody][0]*omega[ibody][0] + 
       angmom[ibody][1]*omega[ibody][1] + angmom[ibody][2]*omega[ibody][2];
@@ -397,12 +398,12 @@ void FixRigidNVT::final_integrate()
     
     // convert torque to the body frame 
     
-    matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody],
-		torque[ibody],tbody);
+    MathExtra::matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody],
+			   torque[ibody],tbody);
     
     // compute "force" for quaternion
     
-    quatvec(quat[ibody],tbody,fquat);
+    MathExtra::quatvec(quat[ibody],tbody,fquat);
     
     // update the conjugate quaternion momentum (conjqm)
     
@@ -411,20 +412,21 @@ void FixRigidNVT::final_integrate()
     conjqm[ibody][2] = scale_r * conjqm[ibody][2] + dtf2 * fquat[2];
     conjqm[ibody][3] = scale_r * conjqm[ibody][3] + dtf2 * fquat[3];
     
-    // compute angular momentum in the body frame then convert to the space-fixed frame
+    // compute angular momentum in the body frame
+    // then convert to the space-fixed frame
+    
+    MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody);
+    MathExtra::matvec_cols(ex_space[ibody],ey_space[ibody],ez_space[ibody],
+			   mbody,angmom[ibody]);
     
-    invquatvec(quat[ibody],conjqm[ibody],mbody);
-    matvec_cols(ex_space[ibody],ey_space[ibody],ez_space[ibody],
-		mbody,angmom[ibody]);
-      
     angmom[ibody][0] *= 0.5;
     angmom[ibody][1] *= 0.5;
     angmom[ibody][2] *= 0.5;  
     
     // compute new angular velocity
     
-    omega_from_angmom(angmom[ibody],ex_space[ibody],ey_space[ibody],
-		  ez_space[ibody],inertia[ibody],omega[ibody]);
+    MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody],
+			       ez_space[ibody],inertia[ibody],omega[ibody]);
   }
   
   // set velocity/rotation of atoms in rigid bodies
diff --git a/src/math_extra.h b/src/math_extra.h
index d3e37345c3..26c6c7ed2d 100755
--- a/src/math_extra.h
+++ b/src/math_extra.h
@@ -64,8 +64,8 @@ namespace MathExtra {
 
   void write3(const double mat[3][3]);
   int mldivide3(const double mat[3][3], const double *vec, double *ans);
-  int jacobi(double **matrix, double *evalues, double **evectors);
-  void rotate(double **matrix, int i, int j, int k, int l,
+  int jacobi(double matrix[3][3], double *evalues, double evectors[3][3]);
+  void rotate(double matrix[3][3], int i, int j, int k, int l,
 	      double s, double tau);
   
   // shape matrix operations
@@ -75,13 +75,8 @@ namespace MathExtra {
                                    double *ans);
   // quaternion operations
   
-  inline void normalize4(double *quat);
   inline void axisangle_to_quat(const double *v, const double angle,
                                 double *quat);
-  inline void multiply_quat_quat(const double *one, const double *two,
-                                 double *ans);
-  inline void multiply_vec_quat(const double *one, const double *two,
-                                double *ans);
   inline void vecquat(double *a, double *b, double *c);
   inline void quatvec(double *a, double *b, double *c);
   inline void quatquat(double *a, double *b, double *c);
@@ -93,10 +88,10 @@ namespace MathExtra {
   inline void matvec_cols(double *x, double *y, double *z,
 			  double *b, double *c);
 
-  void omega_from_angmom(double *m, double *ex, double *ey, double *ez,
-			 double *idiag, double *w);
-  void angmom_from_omega(double *w, double *ex, double *ey, double *ez,
-			 double *idiag, double *m);
+  void angmom_to_omega(double *m, double *ex, double *ey, double *ez,
+		       double *idiag, double *w);
+  void omega_to_angmom(double *w, double *ex, double *ey, double *ez,
+		       double *idiag, double *m);
   void exyz_to_q(double *ex, double *ey, double *ez, double *q);
   void q_to_exyz(double *q, double *ex, double *ey, double *ez);
   void quat_to_mat(const double *quat, double mat[3][3]);
@@ -424,20 +419,6 @@ void MathExtra::multiply_shape_shape(const double *one, const double *two,
   ans[5] = one[0]*two[5] + one[5]*two[1];
 }
 
-/* ----------------------------------------------------------------------
-   normalize a quaternion
-------------------------------------------------------------------------- */
-
-void MathExtra::normalize4(double *quat)
-{
-  double scale = 1.0/sqrt(quat[0]*quat[0]+quat[1]*quat[1] +
-			  quat[2]*quat[2]+quat[3]*quat[3]);
-  quat[0] *= scale;
-  quat[1] *= scale;
-  quat[2] *= scale;
-  quat[3] *= scale;
-}
-
 /* ----------------------------------------------------------------------
    compute quaternion from axis-angle rotation
    v MUST be a unit vector
@@ -454,35 +435,6 @@ void MathExtra::axisangle_to_quat(const double *v, const double angle,
   quat[3] = v[2]*sina;
 }
 
-/* ----------------------------------------------------------------------
-   multiply 2 quaternions
-   effectively concatenates rotations
-   NOT a commutative operation
-------------------------------------------------------------------------- */
-
-void MathExtra::multiply_quat_quat(const double *one, const double *two,
-				   double *ans)
-{
-  ans[0] = one[0]*two[0]-one[1]*two[1]-one[2]*two[2]-one[3]*two[3];
-  ans[1] = one[0]*two[1]+one[1]*two[0]+one[2]*two[3]-one[3]*two[2];
-  ans[2] = one[0]*two[2]-one[1]*two[3]+one[2]*two[0]+one[3]*two[1];
-  ans[3] = one[0]*two[3]+one[1]*two[2]-one[2]*two[1]+one[3]*two[0];
-}
-
-/* ----------------------------------------------------------------------
-   multiply 3 vector times quaternion
-   3 vector one is treated as quaternion (0,one)
-------------------------------------------------------------------------- */
-
-void MathExtra::multiply_vec_quat(const double *one, const double *two,
-				  double *ans)
-{
-  ans[0] = -one[0]*two[1]-one[1]*two[2]-one[2]*two[3];
-  ans[1] =  one[0]*two[0]+one[1]*two[3]-one[2]*two[2];
-  ans[2] = -one[0]*two[3]+one[1]*two[0]+one[2]*two[1];
-  ans[3] =  one[0]*two[2]-one[1]*two[1]+one[2]*two[0];
-}
-
 /* ----------------------------------------------------------------------
    vector-quaternion multiply: c = a*b, where a = (0,a)
 ------------------------------------------------------------------------- */
@@ -509,6 +461,7 @@ void MathExtra::quatvec(double *a, double *b, double *c)
 
 /* ----------------------------------------------------------------------
    quaternion-quaternion multiply: c = a*b
+   NOT a commutative operation
 ------------------------------------------------------------------------- */
 
 void MathExtra::quatquat(double *a, double *b, double *c)
diff --git a/src/set.cpp b/src/set.cpp
index 5be7f0beb8..ee7cc89fce 100644
--- a/src/set.cpp
+++ b/src/set.cpp
@@ -431,7 +431,7 @@ void Set::set(int keyword)
       quat[1] = xvalue * sintheta2;
       quat[2] = yvalue * sintheta2;
       quat[3] = zvalue * sintheta2;
-      MathExtra::normalize4(quat);
+      MathExtra::qnormalize(quat);
     }
     count++;
   }