diff --git a/src/fix_rigid.cpp b/src/fix_rigid.cpp
index e4e3921e17..8c0ed5663a 100644
--- a/src/fix_rigid.cpp
+++ b/src/fix_rigid.cpp
@@ -16,7 +16,9 @@
 #include "stdlib.h"
 #include "string.h"
 #include "fix_rigid.h"
+#include "math_extra.h"
 #include "atom.h"
+#include "atom_vec.h"
 #include "domain.h"
 #include "update.h"
 #include "respa.h"
@@ -50,12 +52,16 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
 
   // perform initial allocation of atom-based arrays
   // register with Atom class
-  
+
+  extended = dorientflag = qorientflag = 0;
   body = NULL;
   displace = NULL;
+  eflags = NULL;
+  dorient = NULL;
+  qorient = NULL;
   grow_arrays(atom->nmax);
   atom->add_callback(0);
-  
+
   // parse command-line args
   // set nbody and body[i] for each atom
 
@@ -299,6 +305,17 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
   for (ibody = 0; ibody < nbody; ibody++)
     image[ibody] = (512 << 20) | (512 << 10) | 512;
 
+  // bitmasks for properties of extended particles
+
+  INERTIA_SPHERE_RADIUS = 1;
+  INERTIA_SPHERE_SHAPE = 2;
+  INERTIA_ELLIPSOID = 4;
+  ORIENT_DIPOLE = 8;
+  ORIENT_QUAT = 16;
+  OMEGA = 32;
+  ANGMOM = 64;
+  TORQUE = 128;
+
   // print statistics
 
   int nsum = 0;
@@ -322,7 +339,10 @@ FixRigid::~FixRigid()
   
   memory->sfree(body);
   memory->destroy_2d_double_array(displace);
-
+  memory->sfree(eflags);
+  memory->destroy_2d_double_array(dorient);
+  memory->destroy_2d_double_array(qorient);
+  
   // delete nbody-length arrays
 
   memory->sfree(nrigid);
@@ -364,7 +384,7 @@ int FixRigid::setmask()
 
 void FixRigid::init()
 {
-  int i,ibody;
+  int i,itype,ibody;
 
   triclinic = domain->triclinic;
 
@@ -396,14 +416,80 @@ void FixRigid::init()
   if (strcmp(update->integrate_style,"respa") == 0)
     step_respa = ((Respa *) update->integrate)->step;
 
+  // extended = 1 if any particle in a rigid body is finite size
+
+  double *radius = atom->radius;
+  double *rmass = atom->rmass;
+  double *mass = atom->mass;
+  double **shape = atom->shape;
+  double *dipole = atom->dipole;
+  int *type = atom->type;
+  int nlocal = atom->nlocal;
+
+  if (atom->radius_flag || atom->avec->shape_type) {
+    int flag = 0;
+    for (i = 0; i < nlocal; i++) {
+      if (body[i] < 0) continue;
+      if (radius && radius[i] > 0.0) flag = 1;
+      else if (shape && shape[type[i]][0] > 0.0) flag = 1;
+    }
+
+    MPI_Allreduce(&flag,&extended,1,MPI_INT,MPI_MAX,world);
+  }
+
+  // grow extended arrays and set extended flags for each particle
+  // vorientflag = 1 if any particles store dipole orientation
+  // qorientflag = 1 if any particles store quat orientation
+
+  if (extended) {
+    if (atom->mu_flag) dorientflag = 1;
+    if (atom->quat_flag) qorientflag = 1;
+    grow_arrays(atom->nmax);
+
+    for (i = 0; i < nlocal; i++) {
+      eflags[i] = 0;
+      if (body[i] < 0) continue;
+
+      // set INERTIA if radius or shape > 0.0
+
+      if (radius) {
+	if (radius[i] > 0.0) eflags[i] |= INERTIA_SPHERE_RADIUS;
+      } else if (shape) {
+	if (shape[type[i]][0] > 0.0) {
+	  if (shape[type[i]][0] == shape[type[i]][1] &&
+	      shape[type[i]][0] == shape[type[i]][2])
+	    eflags[i] |= INERTIA_SPHERE_SHAPE;
+	  else eflags[i] |= INERTIA_ELLIPSOID;
+	}
+      }
+
+      // other flags only set if particle is finite size
+      // set DIPOLE if atom->mu and atom->dipole exist and dipole[itype] > 0.0
+      // set QUAT if atom->quat exists (could be ellipsoid or sphere)
+      // set TORQUE if atom->torque exists
+      // set exactly one of OMEGA and ANGMOM so particle contributes once
+      // set OMEGA if either radius or rmass exists
+      // set ANGMOM if shape and mass exist
+      // set OMEGE if atom->angmom doesn't exist
+
+      if (eflags[i] == 0) continue;
+
+      if (atom->mu_flag && dipole && dipole[type[i]] > 0.0)
+	eflags[i] |= ORIENT_DIPOLE;
+      if (atom->quat_flag) eflags[i] |= ORIENT_QUAT;
+      if (atom->torque_flag) eflags[i] |= TORQUE;
+      if ((radius || rmass) && atom->omega_flag) eflags[i] |= OMEGA;
+      else if (shape && mass && atom->angmom_flag) eflags[i] |= ANGMOM;
+      else if (atom->omega_flag) eflags[i] != OMEGA;
+      else error->one("Could not set finite-size particle attribute "
+		      "in fix rigid");
+    }
+  }
+
   // compute masstotal & center-of-mass of each rigid body
   
   double **x = atom->x;
-  double *rmass = atom->rmass;
-  double *mass = atom->mass;
   int *image = atom->image;
-  int *type = atom->type;
-  int nlocal = atom->nlocal;
   
   double xprd = domain->xprd;
   double yprd = domain->yprd;
@@ -426,7 +512,7 @@ void FixRigid::init()
     zbox = (image[i] >> 20) - 512;
     if (rmass) massone = rmass[i];
     else massone = mass[type[i]];
-    
+
     if (triclinic == 0) {
       xunwrap = x[i][0] + xbox*xprd;
       yunwrap = x[i][1] + ybox*yprd;
@@ -460,10 +546,11 @@ void FixRigid::init()
   // compute 6 moments of inertia of each body
   // dx,dy,dz = coords relative to center-of-mass
   
+  double dx,dy,dz,rad;
+
   for (ibody = 0; ibody < nbody; ibody++)
     for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
-  double dx,dy,dz;
-  
+
   for (i = 0; i < nlocal; i++) {
     if (body[i] < 0) continue;
     ibody = body[i];
@@ -485,9 +572,10 @@ void FixRigid::init()
     dx = xunwrap - xcm[ibody][0];
     dy = yunwrap - xcm[ibody][1];
     dz = zunwrap - xcm[ibody][2];
+
     if (rmass) massone = rmass[i];
     else massone = mass[type[i]];
-    
+
     sum[ibody][0] += massone * (dy*dy + dz*dz);
     sum[ibody][1] += massone * (dx*dx + dz*dz);
     sum[ibody][2] += massone * (dx*dx + dy*dy);
@@ -495,6 +583,55 @@ void FixRigid::init()
     sum[ibody][4] -= massone * dy*dz;
     sum[ibody][5] -= massone * dx*dz;
   }
+
+  // extended particles 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 *radius = atom->radius;
+    double **shape = atom->shape;
+
+    for (i = 0; i < nlocal; i++) {
+      if (body[i] < 0) continue;
+      ibody = body[i];
+
+      itype = type[i];
+      if (rmass) massone = rmass[i];
+      else massone = mass[itype];
+
+      if (eflags[i] & INERTIA_SPHERE_RADIUS) {
+	sum[ibody][0] += 0.4 * massone * radius[i]*radius[i];
+	sum[ibody][1] += 0.4 * massone * radius[i]*radius[i];
+	sum[ibody][2] += 0.4 * massone * radius[i]*radius[i];
+      }
+      if (eflags[i] & INERTIA_SPHERE_SHAPE) {
+	rad = shape[type[i]][0];
+	sum[ibody][0] += 0.4 * massone * rad*rad;
+	sum[ibody][1] += 0.4 * massone * rad*rad;
+	sum[ibody][2] += 0.4 * massone * rad*rad;
+      }
+      if (eflags[i] & INERTIA_ELLIPSOID) {
+	MathExtra::quat_to_mat(atom->quat[i],p);
+	MathExtra::quat_to_mat_trans(atom->quat[i],ptrans);
+	idiag[0] = 0.2*massone *
+	  (shape[itype][1]*shape[itype][1]+shape[itype][2]*shape[itype][2]);
+	idiag[1] = 0.2*massone *
+	  (shape[itype][0]*shape[itype][0]+shape[itype][2]*shape[itype][2]);
+	idiag[2] = 0.2*massone *
+	  (shape[itype][0]*shape[itype][0]+shape[itype][1]*shape[itype][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];
+      }
+    }
+  }
   
   MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
   
@@ -569,36 +706,65 @@ void FixRigid::init()
   
   // 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
+
+  double qc[4];
+
   for (i = 0; i < nlocal; i++) {
-    if (body[i] < 0) displace[i][0] = displace[i][1] = displace[i][2] = 0.0;
-    else {
-      ibody = body[i];
+    if (body[i] < 0) {
+      displace[i][0] = displace[i][1] = displace[i][2] = 0.0;
+      continue;
+    }
 
-      xbox = (image[i] & 1023) - 512;
-      ybox = (image[i] >> 10 & 1023) - 512;
-      zbox = (image[i] >> 20) - 512;
+    ibody = body[i];
 
-      if (triclinic == 0) {
-	xunwrap = x[i][0] + xbox*xprd;
-	yunwrap = x[i][1] + ybox*yprd;
-	zunwrap = x[i][2] + zbox*zprd;
-      } else {
-	xunwrap = x[i][0] + xbox*xprd + ybox*xy + zbox*xz;
-	yunwrap = x[i][1] + ybox*yprd + zbox*yz;
-	zunwrap = x[i][2] + zbox*zprd;
-      }
+    xbox = (image[i] & 1023) - 512;
+    ybox = (image[i] >> 10 & 1023) - 512;
+    zbox = (image[i] >> 20) - 512;
+    
+    if (triclinic == 0) {
+      xunwrap = x[i][0] + xbox*xprd;
+      yunwrap = x[i][1] + ybox*yprd;
+      zunwrap = x[i][2] + zbox*zprd;
+    } else {
+      xunwrap = x[i][0] + xbox*xprd + ybox*xy + zbox*xz;
+      yunwrap = x[i][1] + ybox*yprd + zbox*yz;
+      zunwrap = x[i][2] + zbox*zprd;
+    }
+    
+    dx = xunwrap - xcm[ibody][0];
+    dy = yunwrap - xcm[ibody][1];
+    dz = zunwrap - xcm[ibody][2];
+    
+    displace[i][0] = dx*ex_space[ibody][0] + dy*ex_space[ibody][1] +
+      dz*ex_space[ibody][2];
+    displace[i][1] = dx*ey_space[ibody][0] + dy*ey_space[ibody][1] +
+      dz*ey_space[ibody][2];
+    displace[i][2] = dx*ez_space[ibody][0] + dy*ez_space[ibody][1] +
+      dz*ez_space[ibody][2];
+    
+    if (extended) {
+      double **mu = atom->mu;
+      double *dipole = atom->dipole;
+      int *type = atom->type;
 
-      dx = xunwrap - xcm[ibody][0];
-      dy = yunwrap - xcm[ibody][1];
-      dz = zunwrap - xcm[ibody][2];
-      
-      displace[i][0] = dx*ex_space[ibody][0] + dy*ex_space[ibody][1] +
-	dz*ex_space[ibody][2];
-      displace[i][1] = dx*ey_space[ibody][0] + dy*ey_space[ibody][1] +
-	dz*ey_space[ibody][2];
-      displace[i][2] = dx*ez_space[ibody][0] + dy*ez_space[ibody][1] +
-	dz*ez_space[ibody][2];
+      if (eflags[i] & ORIENT_DIPOLE) {
+	dorient[i][0] = mu[i][0]*ex_space[ibody][0] + 
+	  mu[i][1]*ex_space[ibody][1] + mu[i][2]*ex_space[ibody][2];
+	dorient[i][1] = mu[i][0]*ey_space[ibody][0] + 
+	  mu[i][1]*ey_space[ibody][1] + mu[i][2]*ey_space[ibody][2];
+	dorient[i][2] = mu[i][0]*ez_space[ibody][0] + 
+	  mu[i][1]*ez_space[ibody][1] + mu[i][2]*ez_space[ibody][2];
+	MathExtra::snormalize3(dipole[type[i]],dorient[i],dorient[i]);
+      } else if (dorientflag)
+	dorient[i][0] = dorient[i][1] = dorient[i][2] = 0.0;
+
+      if (eflags[i] & ORIENT_QUAT) {
+	qconjugate(quat[ibody],qc);
+	quatquat(qc,atom->quat[i],qorient[i]);
+	qnormalize(qorient[i]);
+      } else if (qorientflag)
+	qorient[i][0] = qorient[i][1] = qorient[i][2] = qorient[i][3] = 0.0;
     }
   }
   
@@ -606,6 +772,7 @@ void FixRigid::init()
   // recompute moments of inertia around new axes
   // 3 diagonal moments should equal principal moments
   // 3 off-diagonal moments should be 0.0
+  // extended particles contribute extra terms to moments of inertia
   
   for (ibody = 0; ibody < nbody; ibody++)
     for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
@@ -626,7 +793,54 @@ void FixRigid::init()
     sum[ibody][4] -= massone * displace[i][1]*displace[i][2];
     sum[ibody][5] -= massone * displace[i][0]*displace[i][2];
   }
-  
+
+  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 *radius = atom->radius;
+    double **shape = atom->shape;
+
+    for (i = 0; i < nlocal; i++) {
+      if (body[i] < 0) continue;
+      ibody = body[i];
+
+      itype = type[i];
+      if (rmass) massone = rmass[i];
+      else massone = mass[itype];
+
+      if (eflags[i] & INERTIA_SPHERE_RADIUS) {
+	sum[ibody][0] += 0.4 * massone * radius[i]*radius[i];
+	sum[ibody][1] += 0.4 * massone * radius[i]*radius[i];
+	sum[ibody][2] += 0.4 * massone * radius[i]*radius[i];
+      }
+      if (eflags[i] & INERTIA_SPHERE_SHAPE) {
+	rad = shape[type[i]][0];
+	sum[ibody][0] += 0.4 * massone * rad*rad;
+	sum[ibody][1] += 0.4 * massone * rad*rad;
+	sum[ibody][2] += 0.4 * massone * rad*rad;
+      }
+      if (eflags[i] & INERTIA_ELLIPSOID) {
+	MathExtra::quat_to_mat(qorient[i],p);
+	MathExtra::quat_to_mat_trans(qorient[i],ptrans);
+	idiag[0] = 0.2*massone *
+	  (shape[itype][1]*shape[itype][1]+shape[itype][2]*shape[itype][2]);
+	idiag[1] = 0.2*massone *
+	  (shape[itype][0]*shape[itype][0]+shape[itype][2]*shape[itype][2]);
+	idiag[2] = 0.2*massone *
+	  (shape[itype][0]*shape[itype][0]+shape[itype][1]*shape[itype][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];
+      }
+    }
+  }
+
   MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
 
   for (ibody = 0; ibody < nbody; ibody++) {
@@ -663,20 +877,20 @@ void FixRigid::init()
 void FixRigid::setup(int vflag)
 {
   int i,n,ibody;
+  double massone,radone;
   
   // vcm = velocity of center-of-mass of each rigid body
   // fcm = force on center-of-mass of each rigid body
 
-  int *type = atom->type;
   double **v = atom->v;
   double **f = atom->f;
   double *rmass = atom->rmass;
   double *mass = atom->mass;
+  int *type = atom->type;
   int nlocal = atom->nlocal;
 
   for (ibody = 0; ibody < nbody; ibody++)
     for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
-  double massone;
   
   for (i = 0; i < nlocal; i++) {
     if (body[i] < 0) continue;
@@ -753,6 +967,42 @@ void FixRigid::setup(int vflag)
     sum[ibody][4] += dz * f[i][0] - dx * f[i][2];
     sum[ibody][5] += dx * f[i][1] - dy * f[i][0];
   }
+
+  // extended particles add their rotation/torque to angmom/torque of body
+
+  if (extended) {
+    double **omega_one = atom->omega;
+    double **angmom_one = atom->angmom;
+    double **torque_one = atom->torque;
+    double *radius = atom->radius;
+    double **shape = atom->shape;
+
+    for (i = 0; i < nlocal; i++) {
+      if (body[i] < 0) continue;
+      ibody = body[i];
+
+      if (eflags[i] & OMEGA) {
+	if (rmass) massone = rmass[i];
+	else massone = mass[type[i]];
+	if (radius) radone = radius[i];
+	else radone = shape[type[i]][0];
+	sum[ibody][0] += 0.4 * massone * radone*radone * omega_one[i][0];
+	sum[ibody][1] += 0.4 * massone * radone*radone * omega_one[i][1];
+	sum[ibody][2] += 0.4 * massone * radone*radone * omega_one[i][2];
+      }
+
+      if (eflags[i] & ANGMOM) {
+	sum[ibody][0] += angmom_one[i][0];
+	sum[ibody][1] += angmom_one[i][1];
+	sum[ibody][2] += angmom_one[i][2];
+      }
+      if (eflags[i] & TORQUE) {
+	sum[ibody][3] += torque_one[i][0];
+	sum[ibody][4] += torque_one[i][1];
+	sum[ibody][5] += torque_one[i][2];
+      }
+    }
+  }
   
   MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
 
@@ -777,7 +1027,7 @@ void FixRigid::setup(int vflag)
 		      ez_space[ibody],inertia[ibody],omega[ibody]);
   set_v();
 
-  // guestimate virial as 2x the set_v contribution
+  // guesstimate virial as 2x the set_v contribution
 
   if (vflag_global)
     for (n = 0; n < 6; n++) virial[n] *= 2.0;
@@ -829,7 +1079,7 @@ void FixRigid::initial_integrate(int vflag)
   if (vflag) v_setup(vflag);
   else evflag = 0;
 
-  // set coords and velocities of atoms in rigid bodies
+  // set coords/orient and velocity/rotation of atoms in rigid bodies
   // from quarternion and omega
   
   set_xv();
@@ -856,7 +1106,7 @@ void FixRigid::richardson(double *q, double *w,
   qfull[2] = q[2] + dtq * wq[2];
   qfull[3] = q[3] + dtq * wq[3];
 
-  normalize(qfull);
+  qnormalize(qfull);
 
   // 1st half update from dq/dt = 1/2 w q
 
@@ -866,7 +1116,7 @@ void FixRigid::richardson(double *q, double *w,
   qhalf[2] = q[2] + 0.5*dtq * wq[2];
   qhalf[3] = q[3] + 0.5*dtq * wq[3];
 
-  normalize(qhalf);
+  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
@@ -883,7 +1133,7 @@ void FixRigid::richardson(double *q, double *w,
   qhalf[2] += 0.5*dtq * wq[2];
   qhalf[3] += 0.5*dtq * wq[3];
 
-  normalize(qhalf);
+  qnormalize(qhalf);
 
   // corrected Richardson update
 
@@ -892,7 +1142,7 @@ void FixRigid::richardson(double *q, double *w,
   q[2] = 2.0*qhalf[2] - qfull[2];
   q[3] = 2.0*qhalf[3] - qfull[3];
 
-  normalize(q);
+  qnormalize(q);
   exyz_from_q(q,ex,ey,ez);
 }
 
@@ -904,7 +1154,7 @@ void FixRigid::final_integrate()
   double dtfm;
   double xy,xz,yz;
 
-  // sum forces and torques on atoms in rigid body
+  // sum over atoms to get force and torque on rigid body
   
   int *image = atom->image;
   double **x = atom->x;
@@ -955,7 +1205,24 @@ void FixRigid::final_integrate()
     sum[ibody][4] += dz*f[i][0] - dx*f[i][2];
     sum[ibody][5] += dx*f[i][1] - dy*f[i][0];
   }
-  
+
+  // extended particles add their torque to torque of body
+
+  if (extended) {
+    double **torque_one = atom->torque;
+
+    for (i = 0; i < nlocal; i++) {
+      if (body[i] < 0) continue;
+      ibody = body[i];
+
+      if (eflags[i] & TORQUE) {
+	sum[ibody][3] += torque_one[i][0];
+	sum[ibody][4] += torque_one[i][1];
+	sum[ibody][5] += torque_one[i][2];
+      }
+    }
+  }
+
   MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
   
   for (ibody = 0; ibody < nbody; ibody++) {
@@ -980,7 +1247,7 @@ void FixRigid::final_integrate()
     angmom[ibody][2] += dtf * torque[ibody][2] * tflag[ibody][2];
   }
 
-  // set velocities from angmom & omega
+  // set velocity/rotation of atoms in rigid bodies
   // virial is already setup from initial_integrate
 
   for (ibody = 0; ibody < nbody; ibody++) 
@@ -1238,7 +1505,7 @@ void FixRigid::rotate(double **matrix, int i, int j, int k, int l,
 }
 
 /* ----------------------------------------------------------------------
-   create quaternion from space-frame ex,ey,ez
+   create unit quaternion from space-frame ex,ey,ez
    ex,ey,ez are columns of a rotation matrix
 ------------------------------------------------------------------------- */
 
@@ -1277,7 +1544,7 @@ void FixRigid::q_from_exyz(double *ex, double *ey, double *ez, double *q)
   } else
     error->all("Quaternion creation numeric error");
 
-  normalize(q);
+  qnormalize(q);
 }
 
 /* ----------------------------------------------------------------------
@@ -1331,17 +1598,30 @@ void FixRigid::quatvec(double *a, double *b, double *c)
 
 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[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];
 }
 
+/* ----------------------------------------------------------------------
+   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::normalize(double *q)
+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;
@@ -1355,23 +1635,23 @@ void FixRigid::normalize(double *q)
    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 / Ibody
+   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 *inertia, double *w)
+				 double *idiag, double *w)
 {
   double wbody[3];
 
-  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];
+  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];
@@ -1383,18 +1663,19 @@ void FixRigid::omega_from_angmom(double *m, double *ex, double *ey, double *ez,
    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 = Ibody wbody
+   Mbody = Idiag wbody
    Mspace = P Mbody
 ------------------------------------------------------------------------- */
 
-void FixRigid::angmom_from_omega(double *w, double *ex, double *ey, double *ez,
-				 double *inertia, double *m)
+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]) * inertia[0];
-  mbody[1] = (w[0]*ey[0] + w[1]*ey[1] + w[2]*ey[2]) * inertia[1];
-  mbody[2] = (w[0]*ez[0] + w[1]*ez[1] + w[2]*ez[2]) * inertia[2];
+  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];
@@ -1403,17 +1684,18 @@ void FixRigid::angmom_from_omega(double *w, double *ex, double *ey, double *ez,
 
 /* ----------------------------------------------------------------------
    set space-frame coords and velocity of each atom in each rigid body
+   set orientation and rotation of extended particles
    x = Q displace + Xcm, mapped back to periodic box
    v = Vcm + (W cross (x - Xcm))
 ------------------------------------------------------------------------- */
 
 void FixRigid::set_xv()
 {
-  int ibody;
+  int ibody,itype;
   int xbox,ybox,zbox;
   double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone;
   double xy,xz,yz;
-  double vr[6];
+  double ione[3],exone[3],eyone[3],ezone[3],vr[6],p[3][3];
 
   int *image = atom->image;
   double **x = atom->x;
@@ -1504,7 +1786,6 @@ void FixRigid::set_xv()
     if (evflag) {
       if (rmass) massone = rmass[i];
       else massone = mass[type[i]];
-
       fc0 = massone*(v[i][0] - v0)/dtf - f[i][0];
       fc1 = massone*(v[i][1] - v1)/dtf - f[i][1];
       fc2 = massone*(v[i][2] - v2)/dtf - f[i][2]; 
@@ -1519,25 +1800,66 @@ void FixRigid::set_xv()
       v_tally(1,&i,1.0,vr);
     }
   }
+
+  // set orientation, omega, angmom of each extended particle
+  
+  if (extended) {
+    double **omega_one = atom->omega;
+    double **angmom_one = atom->angmom;
+    double *dipole = atom->dipole;
+    double **shape = atom->shape;
+
+    for (int i = 0; i < nlocal; i++) {
+      if (body[i] < 0) continue;
+      ibody = body[i];
+      
+      if (eflags[i] & ORIENT_DIPOLE) {
+	MathExtra::quat_to_mat(quat[ibody],p);
+	MathExtra::times_column3(p,dorient[i],atom->mu[i]);
+	MathExtra::snormalize3(dipole[type[i]],atom->mu[i],atom->mu[i]);
+      }
+      if (eflags[i] & ORIENT_QUAT) {
+	quatquat(quat[ibody],qorient[i],atom->quat[i]);
+	qnormalize(atom->quat[i]);
+      }
+      if (eflags[i] & OMEGA) {
+	omega_one[i][0] = omega[ibody][0];
+	omega_one[i][1] = omega[ibody][1];
+	omega_one[i][2] = omega[ibody][2];
+      }
+      if (eflags[i] & ANGMOM) {
+	itype = type[i];
+	ione[0] = 0.2*mass[itype] *
+	  (shape[itype][1]*shape[itype][1] + shape[itype][2]*shape[itype][2]);
+	ione[1] = 0.2*mass[itype] *
+	  (shape[itype][0]*shape[itype][0] + shape[itype][2]*shape[itype][2]);
+	ione[2] = 0.2*mass[itype] * 
+	  (shape[itype][0]*shape[itype][0] + shape[itype][1]*shape[itype][1]);
+	exyz_from_q(atom->quat[i],exone,eyone,ezone);
+	angmom_from_omega(omega[ibody],exone,eyone,ezone,ione,angmom_one[i]);
+      }
+    }
+  }
 }
 
 /* ----------------------------------------------------------------------
    set space-frame velocity of each atom in a rigid body
+   set omega and angmom of extended particles
    v = Vcm + (W cross (x - Xcm))
 ------------------------------------------------------------------------- */
 
 void FixRigid::set_v()
 {
-  int ibody;
+  int ibody,itype;
   int xbox,ybox,zbox;
   double dx,dy,dz;
   double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone;
   double xy,xz,yz;
-  double vr[6];
+  double ione[3],exone[3],eyone[3],ezone[3],vr[6];
 
-  double **f = atom->f;
-  double **v = atom->v;
   double **x = atom->x;
+  double **v = atom->v;
+  double **f = atom->f;
   double *rmass = atom->rmass; 
   double *mass = atom->mass; 
   int *type = atom->type;
@@ -1590,7 +1912,6 @@ void FixRigid::set_v()
     if (evflag) {
       if (rmass) massone = rmass[i];
       else massone = mass[type[i]];
-
       fc0 = massone*(v[i][0] - v0)/dtf - f[i][0];
       fc1 = massone*(v[i][1] - v1)/dtf - f[i][1];
       fc2 = massone*(v[i][2] - v2)/dtf - f[i][2]; 
@@ -1619,6 +1940,36 @@ void FixRigid::set_v()
       v_tally(1,&i,1.0,vr);
     }
   }
+
+  // set omega, angmom of each extended particle
+  
+  if (extended) {
+    double **omega_one = atom->omega;
+    double **angmom_one = atom->angmom;
+    double **shape = atom->shape;
+
+    for (int i = 0; i < nlocal; i++) {
+      if (body[i] < 0) continue;
+      ibody = body[i];
+
+      if (eflags[i] & OMEGA) {
+	omega_one[i][0] = omega[ibody][0];
+	omega_one[i][1] = omega[ibody][1];
+	omega_one[i][2] = omega[ibody][2];
+      }
+      if (eflags[i] & ANGMOM) {
+	itype = type[i];
+	ione[0] = 0.2*mass[itype] *
+	  (shape[itype][1]*shape[itype][1] + shape[itype][2]*shape[itype][2]);
+	ione[1] = 0.2*mass[itype] *
+	  (shape[itype][0]*shape[itype][0] + shape[itype][2]*shape[itype][2]);
+	ione[2] = 0.2*mass[itype] * 
+	  (shape[itype][0]*shape[itype][0] + shape[itype][1]*shape[itype][1]);
+	exyz_from_q(atom->quat[i],exone,eyone,ezone);
+	angmom_from_omega(omega[ibody],exone,eyone,ezone,ione,angmom_one[i]);
+      }
+    }
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -1631,6 +1982,11 @@ double FixRigid::memory_usage()
   double bytes = nmax * sizeof(int);
   bytes += nmax*3 * sizeof(double);
   bytes += maxvatom*6 * sizeof(double);
+  if (extended) {
+    bytes += nmax * sizeof(int);
+    if (dorientflag) bytes = nmax*3 * sizeof(double);
+    if (qorientflag) bytes = nmax*4 * sizeof(double);
+  }
   return bytes;
 } 
 
@@ -1642,6 +1998,14 @@ void FixRigid::grow_arrays(int nmax)
 {
   body = (int *) memory->srealloc(body,nmax*sizeof(int),"rigid:body");
   displace = memory->grow_2d_double_array(displace,nmax,3,"rigid:displace");
+  if (extended) {
+    eflags = (int *)
+      memory->srealloc(eflags,nmax*sizeof(int),"rigid:eflags");
+    if (dorientflag)
+      dorient = memory->grow_2d_double_array(dorient,nmax,3,"rigid:dorient");
+    if (qorientflag)
+      qorient = memory->grow_2d_double_array(qorient,nmax,4,"rigid:qorient");
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -1654,6 +2018,20 @@ void FixRigid::copy_arrays(int i, int j)
   displace[j][0] = displace[i][0];
   displace[j][1] = displace[i][1];
   displace[j][2] = displace[i][2];
+  if (extended) {
+    eflags[j] = eflags[i];
+    if (dorientflag) {
+      dorient[j][0] = dorient[i][0];
+      dorient[j][1] = dorient[i][1];
+      dorient[j][2] = dorient[i][2];
+    }
+    if (qorientflag) {
+      qorient[j][0] = qorient[i][0];
+      qorient[j][1] = qorient[i][1];
+      qorient[j][2] = qorient[i][2];
+      qorient[j][3] = qorient[i][3];
+    }
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -1666,7 +2044,22 @@ int FixRigid::pack_exchange(int i, double *buf)
   buf[1] = displace[i][0];
   buf[2] = displace[i][1];
   buf[3] = displace[i][2];
-  return 4;
+  if (!extended) return 4;
+
+  int m = 4;
+  buf[m++] = eflags[i];
+  if (dorientflag) {
+    buf[m++] = dorient[i][0];
+    buf[m++] = dorient[i][1];
+    buf[m++] = dorient[i][2];
+  }
+  if (qorientflag) {
+    buf[m++] = qorient[i][0];
+    buf[m++] = qorient[i][1];
+    buf[m++] = qorient[i][2];
+    buf[m++] = qorient[i][3];
+  }
+  return m;
 }
 
 /* ----------------------------------------------------------------------
@@ -1679,7 +2072,22 @@ int FixRigid::unpack_exchange(int nlocal, double *buf)
   displace[nlocal][0] = buf[1];
   displace[nlocal][1] = buf[2];
   displace[nlocal][2] = buf[3];
-  return 4;
+  if (!extended) return 4;
+
+  int m = 4;
+  eflags[nlocal] = static_cast<int> (buf[m++]);
+  if (dorientflag) {
+    dorient[nlocal][0] = buf[m++];
+    dorient[nlocal][0] = buf[m++];
+    dorient[nlocal][0] = buf[m++];
+  }
+  if (qorientflag) {
+    qorient[nlocal][0] = buf[m++];
+    qorient[nlocal][0] = buf[m++];
+    qorient[nlocal][0] = buf[m++];
+    qorient[nlocal][0] = buf[m++];
+  }
+  return m;
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/fix_rigid.h b/src/fix_rigid.h
index e7cc5ff86f..5bd3ad783f 100644
--- a/src/fix_rigid.h
+++ b/src/fix_rigid.h
@@ -61,14 +61,28 @@ class FixRigid : public Fix {
   double **torque;          // torque on each rigid body in space coords
   double **quat;            // quaternion of each rigid body
   int *image;               // image flags of xcm of each rigid body
-  int *body;                // which body each atom is part of (-1 if none)
-  double **displace;        // displacement of each atom in body coords
   double **fflag;           // flag for on/off of center-of-mass force
   double **tflag;           // flag for on/off of center-of-mass torque
 
+  int *body;                // which body each atom is part of (-1 if none)
+  double **displace;        // displacement of each atom in body coords
+
   double **sum,**all;       // work vectors for each rigid body
   int **remapflag;          // PBC remap flags for each rigid body
 
+  int extended;             // 1 if any particles have extended attributes
+  int dorientflag;          // 1 if particles store dipole orientation
+  int qorientflag;          // 1 if particles store quat orientation
+
+  int *eflags;              // flags for extended particles
+  double **qorient;         // rotation state of ext particle wrt rigid body
+  double **dorient;         // orientation of dipole mu wrt rigid body
+
+                            // bitmasks for eflags
+  int INERTIA_SPHERE_RADIUS,INERTIA_SPHERE_SHAPE,INERTIA_ELLIPSOID;
+  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 *);
@@ -76,7 +90,8 @@ class FixRigid : public Fix {
   void vecquat(double *, double *, double *);
   void quatvec(double *, double *, double *);
   void quatquat(double *, double *, double *);
-  void normalize(double *);
+  void qconjugate(double *, double *);
+  void qnormalize(double *);
   void richardson(double *, double *, double *, double *,
 		  double *, double *, double *);
   void omega_from_angmom(double *, double *, double *,
diff --git a/src/math_extra.h b/src/math_extra.h
index 68469c0904..649d4910a7 100755
--- a/src/math_extra.h
+++ b/src/math_extra.h
@@ -28,6 +28,7 @@ namespace MathExtra {
   // 3 vector operations
 
   inline void normalize3(const double *v, double *ans);
+  inline void snormalize3(const double, const double *v, double *ans);
   inline double dot3(const double *v1, const double *v2);
   inline void cross3(const double *v1, const double *v2, double *ans);
 
@@ -86,10 +87,22 @@ namespace MathExtra {
 
 void MathExtra::normalize3(const double *v, double *ans)
 {
-  double den = sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
-  ans[0] = v[0]/den;
-  ans[1] = v[1]/den;
-  ans[2] = v[2]/den;
+  double scale = 1.0/sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
+  ans[0] = v[0]*scale;
+  ans[1] = v[1]*scale;
+  ans[2] = v[2]*scale;
+}
+
+/* ----------------------------------------------------------------------
+   scale a vector to length
+------------------------------------------------------------------------- */
+
+void MathExtra::snormalize3(const double length, const double *v, double *ans)
+{
+  double scale = length/sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
+  ans[0] = v[0]*scale;
+  ans[1] = v[1]*scale;
+  ans[2] = v[2]*scale;
 }
 
 /* ----------------------------------------------------------------------
@@ -347,12 +360,12 @@ void MathExtra::write3(const double mat[3][3])
 
 void MathExtra::normalize4(double *quat)
 {
-  double den = sqrt(quat[0]*quat[0]+quat[1]*quat[1] +
-		    quat[2]*quat[2]+quat[3]*quat[3]);
-  quat[0] /= den;
-  quat[1] /= den;
-  quat[2] /= den;
-  quat[3] /= den;
+  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;
 }
 
 /* ----------------------------------------------------------------------