diff --git a/src/RIGID/fix_rigid.cpp b/src/RIGID/fix_rigid.cpp
index c4051d495f..f4c89ea721 100644
--- a/src/RIGID/fix_rigid.cpp
+++ b/src/RIGID/fix_rigid.cpp
@@ -639,10 +639,14 @@ void FixRigid::init()
     step_respa = ((Respa *) update->integrate)->step;
 
   // one-time initialization of rigid body attributes
-  // extended flags, masstotal, COM, inertia tensor
+  // setup_bodies_static = extended flags, masstotal, COM, inertia tensor
+  // setup_bodies_dynamic = vcm and angmom
 
-  if (firstflag) setup_bodies();
-  firstflag = 0;
+  if (firstflag) {
+    firstflag = 0;
+    setup_bodies_static();
+    setup_bodies_dynamic();
+  }
 
   // temperature scale factor
 
@@ -655,20 +659,19 @@ void FixRigid::init()
   else tfactor = 0.0;
 }
 
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+   compute initial fcm and torque on bodies, also initial virial
+   reset all particle velocities to be consistent with vcm and omega
+------------------------------------------------------------------------- */
 
 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
 
-  double **v = atom->v;
   double **f = atom->f;
-  double *rmass = atom->rmass;
-  double *mass = atom->mass;
   int *type = atom->type;
   int nlocal = atom->nlocal;
 
@@ -678,29 +681,19 @@ void FixRigid::setup(int vflag)
   for (i = 0; i < nlocal; i++) {
     if (body[i] < 0) continue;
     ibody = body[i];
-    if (rmass) massone = rmass[i];
-    else massone = mass[type[i]];
-
-    sum[ibody][0] += v[i][0] * massone;
-    sum[ibody][1] += v[i][1] * massone;
-    sum[ibody][2] += v[i][2] * massone;
-    sum[ibody][3] += f[i][0];
-    sum[ibody][4] += f[i][1];
-    sum[ibody][5] += f[i][2];
+    sum[ibody][0] += f[i][0];
+    sum[ibody][1] += f[i][1];
+    sum[ibody][2] += f[i][2];
   }
 
   MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
 
   for (ibody = 0; ibody < nbody; ibody++) {
-    vcm[ibody][0] = all[ibody][0]/masstotal[ibody];
-    vcm[ibody][1] = all[ibody][1]/masstotal[ibody];
-    vcm[ibody][2] = all[ibody][2]/masstotal[ibody];
-    fcm[ibody][0] = all[ibody][3];
-    fcm[ibody][1] = all[ibody][4];
-    fcm[ibody][2] = all[ibody][5];
+    fcm[ibody][0] = all[ibody][0];
+    fcm[ibody][1] = all[ibody][1];
+    fcm[ibody][2] = all[ibody][2];
   }
 
-  // angmom = angular momentum of each rigid body
   // torque = torque on each rigid body
 
   tagint *image = atom->image;
@@ -721,52 +714,23 @@ void FixRigid::setup(int vflag)
     dy = unwrap[1] - xcm[ibody][1];
     dz = unwrap[2] - xcm[ibody][2];
 
-    if (rmass) massone = rmass[i];
-    else massone = mass[type[i]];
-
-    sum[ibody][0] += dy * massone*v[i][2] - dz * massone*v[i][1];
-    sum[ibody][1] += dz * massone*v[i][0] - dx * massone*v[i][2];
-    sum[ibody][2] += dx * massone*v[i][1] - dy * massone*v[i][0];
-    sum[ibody][3] += dy * f[i][2] - dz * f[i][1];
-    sum[ibody][4] += dz * f[i][0] - dx * f[i][2];
-    sum[ibody][5] += dx * f[i][1] - dy * f[i][0];
+    sum[ibody][0] += dy * f[i][2] - dz * f[i][1];
+    sum[ibody][1] += dz * f[i][0] - dx * f[i][2];
+    sum[ibody][2] += dx * f[i][1] - dy * f[i][0];
   }
 
-  // extended particles add their rotation/torque to angmom/torque of body
+  // extended particles add their torque to torque of body
 
   if (extended) {
-    AtomVecLine::Bonus *lbonus;
-    if (avec_line) lbonus = avec_line->bonus;
-    double **omega_one = atom->omega;
-    double **angmom_one = atom->angmom;
     double **torque_one = atom->torque;
-    double *radius = atom->radius;
-    int *line = atom->line;
 
     for (i = 0; i < nlocal; i++) {
       if (body[i] < 0) continue;
       ibody = body[i];
-
-      if (eflags[i] & OMEGA) {
-        if (eflags[i] & SPHERE) {
-          radone = radius[i];
-          sum[ibody][0] += SINERTIA*rmass[i] * radone*radone * omega_one[i][0];
-          sum[ibody][1] += SINERTIA*rmass[i] * radone*radone * omega_one[i][1];
-          sum[ibody][2] += SINERTIA*rmass[i] * radone*radone * omega_one[i][2];
-        } else if (eflags[i] & LINE) {
-          radone = lbonus[line[i]].length;
-          sum[ibody][2] += LINERTIA*rmass[i] * 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];
+        sum[ibody][0] += torque_one[i][0];
+        sum[ibody][1] += torque_one[i][1];
+        sum[ibody][2] += torque_one[i][2];
       }
     }
   }
@@ -774,12 +738,9 @@ void FixRigid::setup(int vflag)
   MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
 
   for (ibody = 0; ibody < nbody; ibody++) {
-    angmom[ibody][0] = all[ibody][0];
-    angmom[ibody][1] = all[ibody][1];
-    angmom[ibody][2] = all[ibody][2];
-    torque[ibody][0] = all[ibody][3];
-    torque[ibody][1] = all[ibody][4];
-    torque[ibody][2] = all[ibody][5];
+    torque[ibody][0] = all[ibody][0];
+    torque[ibody][1] = all[ibody][1];
+    torque[ibody][2] = all[ibody][2];
   }
 
   // zero langextra in case Langevin thermostat not used
@@ -1150,7 +1111,7 @@ void FixRigid::pre_neighbor()
 
 int FixRigid::dof(int tgroup)
 {
-  // cannot count DOF correctly unless setup_bodies() has been called
+  // cannot count DOF correctly unless setup_bodies_static() has been called
 
   if (firstflag) {
     if (comm->me == 0) 
@@ -1561,13 +1522,13 @@ void FixRigid::set_v()
 }
 
 /* ----------------------------------------------------------------------
-   one-time initialization of rigid body attributes
+   one-time initialization of static rigid body attributes
    extended flags, masstotal, center-of-mass
    Cartesian and diagonalized inertia tensor
    read per-body attributes from infile if specified
 ------------------------------------------------------------------------- */
 
-void FixRigid::setup_bodies()
+void FixRigid::setup_bodies_static()
 {
   int i,ibody;
 
@@ -2059,6 +2020,104 @@ void FixRigid::setup_bodies()
   if (infile) memory->destroy(inbody);
 }
 
+/* ----------------------------------------------------------------------
+   one-time initialization of dynamic rigid body attributes
+   Vcm and angmom, computed explicitly from constituent particles
+   even if wrong for overlapping particles, is OK,
+     since is just setting initial time=0 Vcm and angmom of the body
+     which can be estimated value
+------------------------------------------------------------------------- */
+
+void FixRigid::setup_bodies_dynamic()
+{
+  int i,n,ibody;
+  double massone,radone;
+
+  // vcm = velocity of center-of-mass of each rigid body
+  // angmom = angular momentum of each rigid body
+
+  double **x = atom->x;
+  double **v = atom->v;
+  double *rmass = atom->rmass;
+  double *mass = atom->mass;
+  tagint *image = atom->image;
+  int *type = atom->type;
+  int nlocal = atom->nlocal;
+
+  double dx,dy,dz;
+  double unwrap[3];
+
+  for (ibody = 0; ibody < nbody; ibody++)
+    for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
+
+  for (i = 0; i < nlocal; i++) {
+    if (body[i] < 0) continue;
+    ibody = body[i];
+
+    if (rmass) massone = rmass[i];
+    else massone = mass[type[i]];
+
+    sum[ibody][0] += v[i][0] * massone;
+    sum[ibody][1] += v[i][1] * massone;
+    sum[ibody][2] += v[i][2] * massone;
+
+    domain->unmap(x[i],image[i],unwrap);
+    dx = unwrap[0] - xcm[ibody][0];
+    dy = unwrap[1] - xcm[ibody][1];
+    dz = unwrap[2] - xcm[ibody][2];
+
+    sum[ibody][3] += dy * massone*v[i][2] - dz * massone*v[i][1];
+    sum[ibody][4] += dz * massone*v[i][0] - dx * massone*v[i][2];
+    sum[ibody][5] += dx * massone*v[i][1] - dy * massone*v[i][0];
+  }
+
+  // extended particles add their rotation to angmom of body
+
+  if (extended) {
+    AtomVecLine::Bonus *lbonus;
+    if (avec_line) lbonus = avec_line->bonus;
+    double **omega_one = atom->omega;
+    double **angmom_one = atom->angmom;
+    double *radius = atom->radius;
+    int *line = atom->line;
+
+    for (i = 0; i < nlocal; i++) {
+      if (body[i] < 0) continue;
+      ibody = body[i];
+
+      if (eflags[i] & OMEGA) {
+        if (eflags[i] & SPHERE) {
+          radone = radius[i];
+          sum[ibody][3] += SINERTIA*rmass[i] * radone*radone * omega_one[i][0];
+          sum[ibody][4] += SINERTIA*rmass[i] * radone*radone * omega_one[i][1];
+          sum[ibody][5] += SINERTIA*rmass[i] * radone*radone * omega_one[i][2];
+        } else if (eflags[i] & LINE) {
+          radone = lbonus[line[i]].length;
+          sum[ibody][5] += LINERTIA*rmass[i] * radone*radone * omega_one[i][2];
+        }
+      }
+      if (eflags[i] & ANGMOM) {
+        sum[ibody][3] += angmom_one[i][0];
+        sum[ibody][4] += angmom_one[i][1];
+        sum[ibody][5] += angmom_one[i][2];
+      }
+    }
+  }
+
+  MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
+
+  // normalize velocity of COM
+
+  for (ibody = 0; ibody < nbody; ibody++) {
+    vcm[ibody][0] = all[ibody][0]/masstotal[ibody];
+    vcm[ibody][1] = all[ibody][1]/masstotal[ibody];
+    vcm[ibody][2] = all[ibody][2]/masstotal[ibody];
+    angmom[ibody][0] = all[ibody][3];
+    angmom[ibody][1] = all[ibody][4];
+    angmom[ibody][2] = all[ibody][5];
+  }
+}
+
 /* ----------------------------------------------------------------------
    read per rigid body info from user-provided file
    which = 0 to read total mass and center-of-mass, store in vec and array
@@ -2341,200 +2400,32 @@ void FixRigid::reset_dt()
 
 /* ----------------------------------------------------------------------
    zero linear momentum of each rigid body
-   called by velocity zero linear command with its vgroup
-   only atoms in velocity command group contribute to vcm of body,
-     and only their velocities are adjusted, odd if partial bodies specified
+   set Vcm to 0.0, then reset velocities of particles via set_v()
 ------------------------------------------------------------------------- */
 
-void FixRigid::zero_momentum(int vgroup)
+void FixRigid::zero_momentum()
 {
-  int i,ibody;
-  double massone;
+  for (int ibody = 0; ibody < nbody; ibody++)
+    vcm[ibody][0] = vcm[ibody][1] = vcm[ibody][2] = 0.0;
 
-  int vgroupbit = group->bitmask[vgroup];
-
-  // vcm = velocity of center-of-mass of each rigid body
-
-  double **v = atom->v;
-  double *rmass = atom->rmass;
-  double *mass = atom->mass;
-  int *mask = atom->mask;
-  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;
-
-  for (i = 0; i < nlocal; i++) {
-    if (!(mask[i] & vgroupbit)) continue;
-    if (body[i] < 0) continue;
-    ibody = body[i];
-    if (rmass) massone = rmass[i];
-    else massone = mass[type[i]];
-    sum[ibody][0] += v[i][0] * massone;
-    sum[ibody][1] += v[i][1] * massone;
-    sum[ibody][2] += v[i][2] * massone;
-  }
-
-  MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
-
-  for (ibody = 0; ibody < nbody; ibody++) {
-    vcm[ibody][0] = all[ibody][0]/masstotal[ibody];
-    vcm[ibody][1] = all[ibody][1]/masstotal[ibody];
-    vcm[ibody][2] = all[ibody][2]/masstotal[ibody];
-  }
-
-  // adjust velocities by vcm to zero linear momentum
-
-  for (i = 0; i < nlocal; i++) {
-    if (!(mask[i] & vgroupbit)) continue;
-    if (body[i] < 0) continue;
-    ibody = body[i];
-    v[i][0] -= vcm[ibody][0];
-    v[i][1] -= vcm[ibody][1];
-    v[i][2] -= vcm[ibody][2];
-  }
+  evflag = 0;
+  set_v();
 }
 
 /* ----------------------------------------------------------------------
    zero angular momentum of each rigid body
-   called by velocity zero linear command with its vgroup
-   only atoms in velocity command group contribute to angmom/omega of body,
-     and only their velocities are adjusted, odd if partial bodies specified
+   set angmom/omega to 0.0, then reset velocities of particles via set_v()
 ------------------------------------------------------------------------- */
 
-void FixRigid::zero_rotation(int vgroup)
+void FixRigid::zero_rotation()
 {
-  int i,ibody;
-  double massone,radone;
-
-  int vgroupbit = group->bitmask[vgroup];
-
-  // angmom = angular momentum of each rigid body
-
-  double **x = atom->x;
-  double **v = atom->v;
-  double *rmass = atom->rmass;
-  double *mass = atom->mass;
-  tagint *image = atom->image;
-  int *mask = atom->mask;
-  int *type = atom->type;
-  int nlocal = atom->nlocal;
-
-  double dx,dy,dz;
-  double unwrap[3];
-
-  for (ibody = 0; ibody < nbody; ibody++)
-    for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
-
-  for (i = 0; i < nlocal; i++) {
-    if (!(mask[i] & vgroupbit)) continue;
-    if (body[i] < 0) continue;
-    ibody = body[i];
-
-    domain->unmap(x[i],image[i],unwrap);
-    dx = unwrap[0] - xcm[ibody][0];
-    dy = unwrap[1] - xcm[ibody][1];
-    dz = unwrap[2] - xcm[ibody][2];
-
-    if (rmass) massone = rmass[i];
-    else massone = mass[type[i]];
-
-    sum[ibody][0] += dy * massone*v[i][2] - dz * massone*v[i][1];
-    sum[ibody][1] += dz * massone*v[i][0] - dx * massone*v[i][2];
-    sum[ibody][2] += dx * massone*v[i][1] - dy * massone*v[i][0];
+  for (int ibody = 0; ibody < nbody; ibody++) {
+    angmom[ibody][0] = angmom[ibody][1] = angmom[ibody][2] = 0.0;
+    omega[ibody][0] = omega[ibody][1] = omega[ibody][2] = 0.0;
   }
 
-  // extended particles add their rotation to angmom of body
-
-  if (extended) {
-    AtomVecLine::Bonus *lbonus;
-    if (avec_line) lbonus = avec_line->bonus;
-    double **omega_one = atom->omega;
-    double **angmom_one = atom->angmom;
-    double *radius = atom->radius;
-    int *line = atom->line;
-
-    for (i = 0; i < nlocal; i++) {
-      if (body[i] < 0) continue;
-      ibody = body[i];
-
-      if (eflags[i] & OMEGA) {
-        if (eflags[i] & SPHERE) {
-          radone = radius[i];
-          sum[ibody][0] += SINERTIA*rmass[i] * radone*radone * omega_one[i][0];
-          sum[ibody][1] += SINERTIA*rmass[i] * radone*radone * omega_one[i][1];
-          sum[ibody][2] += SINERTIA*rmass[i] * radone*radone * omega_one[i][2];
-        } else if (eflags[i] & LINE) {
-          radone = lbonus[line[i]].length;
-          sum[ibody][2] += LINERTIA*rmass[i] * 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];
-      }
-    }
-  }
-
-  MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
-
-  for (ibody = 0; ibody < nbody; ibody++) {
-    angmom[ibody][0] = all[ibody][0];
-    angmom[ibody][1] = all[ibody][1];
-    angmom[ibody][2] = all[ibody][2];
-  }
-
-  // convert angmom to omega
-
-  for (ibody = 0; ibody < nbody; ibody++)
-    MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody],
-                               ez_space[ibody],inertia[ibody],omega[ibody]);
-
-  // adjust velocities to zero omega
-  // vnew_i = v_i - w x r_i
-  // must use unwrapped coords to compute r_i correctly
-
-  for (i = 0; i < nlocal; i++) {
-    if (!(mask[i] & vgroupbit)) continue;
-    if (body[i] < 0) continue;
-    ibody = body[i];
-
-    domain->unmap(x[i],image[i],unwrap);
-    dx = unwrap[0] - xcm[ibody][0];
-    dy = unwrap[1] - xcm[ibody][1];
-    dz = unwrap[2] - xcm[ibody][2];
-
-    v[i][0] -= omega[ibody][1]*dz - omega[ibody][2]*dy;
-    v[i][1] -= omega[ibody][2]*dx - omega[ibody][0]*dz;
-    v[i][2] -= omega[ibody][0]*dy - omega[ibody][1]*dx;
-  }
-
-  // also explicitly zero omega and angmom of extended particles
-
-  if (extended) {
-    double **omega_one = atom->omega;
-    double **angmom_one = atom->angmom;
-
-    for (i = 0; i < nlocal; i++) {
-      if (!(mask[i] & vgroupbit)) continue;
-      if (body[i] < 0) continue;
-
-      if (eflags[i] & OMEGA) {
-        if (eflags[i] & SPHERE) {
-          omega_one[i][0] = 0.0;
-          omega_one[i][1] = 0.0;
-          omega_one[i][2] = 0.0;
-        } else if (eflags[i] & LINE) omega_one[i][2] = 0.0;
-      }
-      if (eflags[i] & ANGMOM) {
-        angmom_one[i][0] = 0.0;
-        angmom_one[i][1] = 0.0;
-        angmom_one[i][2] = 0.0;
-      }
-    }
-  }
+  evflag = 0;
+  set_v();
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/RIGID/fix_rigid.h b/src/RIGID/fix_rigid.h
index a08d9c87d5..913faf5434 100644
--- a/src/RIGID/fix_rigid.h
+++ b/src/RIGID/fix_rigid.h
@@ -51,8 +51,8 @@ class FixRigid : public Fix {
   int dof(int);
   void deform(int);
   void reset_dt();
-  void zero_momentum(int);
-  void zero_rotation(int);
+  void zero_momentum();
+  void zero_rotation();
   virtual void *extract(const char*, int &);
   double extract_ke();
   double extract_erotational();
@@ -136,7 +136,8 @@ class FixRigid : public Fix {
   void no_squish_rotate(int, double *, double *, double *, double) const;
   void set_xv();
   void set_v();
-  void setup_bodies();
+  void setup_bodies_static();
+  void setup_bodies_dynamic();
   void readfile(int, double *, double **, int *);
 };
 
diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp
index 98a6224d30..e91eccfd78 100644
--- a/src/RIGID/fix_rigid_small.cpp
+++ b/src/RIGID/fix_rigid_small.cpp
@@ -345,11 +345,15 @@ void FixRigidSmall::setup_pre_neighbor()
 {
   if (firstflag) {
     firstflag = 0;
-    setup_bodies();
+    setup_bodies_static();
+    setup_bodies_dynamic();
   } else pre_neighbor();
 }
 
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+   compute initial fcm and torque on bodies, also initial virial
+   reset all particle velocities to be consistent with vcm and omega
+------------------------------------------------------------------------- */
 
 void FixRigidSmall::setup(int vflag)
 {
@@ -358,32 +362,23 @@ void FixRigidSmall::setup(int vflag)
 
   //check(1);
 
-  // sum vcm, fcm, angmom, torque across all rigid bodies
-  // vcm = velocity of COM
+  // sum fcm, torque across all rigid bodies
   // fcm = force on COM
-  // angmom = angular momentum around COM
   // torque = torque around COM
 
   double **x = atom->x;
-  double **v = atom->v;
   double **f = atom->f;
-  double *rmass = atom->rmass;
-  double *mass = atom->mass;
   int *type = atom->type;
   int *image = atom->image;
   int nlocal = atom->nlocal;
 
-  double *xcm,*vcm,*fcm,*acm,*tcm;
+  double *xcm,*fcm,*tcm;
   double dx,dy,dz;
   double unwrap[3];
 
   for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++) {
-    vcm = body[ibody].vcm;
-    vcm[0] = vcm[1] = vcm[2] = 0.0;
     fcm = body[ibody].fcm;
     fcm[0] = fcm[1] = fcm[2] = 0.0;
-    acm = body[ibody].angmom;
-    acm[0] = acm[1] = acm[2] = 0.0;
     tcm = body[ibody].torque;
     tcm[0] = tcm[1] = tcm[2] = 0.0;
   }
@@ -392,13 +387,6 @@ void FixRigidSmall::setup(int vflag)
     if (atom2body[i] < 0) continue;
     Body *b = &body[atom2body[i]];
 
-    if (rmass) massone = rmass[i];
-    else massone = mass[type[i]];
-
-    vcm = b->vcm;
-    vcm[0] += v[i][0] * massone;
-    vcm[1] += v[i][1] * massone;
-    vcm[2] += v[i][2] * massone;
     fcm = b->fcm;
     fcm[0] += f[i][0];
     fcm[1] += f[i][1];
@@ -410,10 +398,6 @@ void FixRigidSmall::setup(int vflag)
     dy = unwrap[1] - xcm[1];
     dz = unwrap[2] - xcm[2];
 
-    acm = b->angmom;
-    acm[0] += dy * massone*v[i][2] - dz * massone*v[i][1];
-    acm[1] += dz * massone*v[i][0] - dx * massone*v[i][2];
-    acm[2] += dx * massone*v[i][1] - dy * massone*v[i][0];
     tcm = b->torque;
     tcm[0] += dy * f[i][2] - dz * f[i][1];
     tcm[1] += dz * f[i][0] - dx * f[i][2];
@@ -423,36 +407,11 @@ void FixRigidSmall::setup(int vflag)
   // extended particles add their rotation/torque to angmom/torque of body
 
   if (extended) {
-    AtomVecLine::Bonus *lbonus;
-    if (avec_line) lbonus = avec_line->bonus;
-    double **omega = atom->omega;
-    double **angmom = atom->angmom;
     double **torque = atom->torque;
-    double *radius = atom->radius;
-    int *line = atom->line;
 
     for (i = 0; i < nlocal; i++) {
       if (atom2body[i] < 0) continue;
       Body *b = &body[atom2body[i]];
-
-      if (eflags[i] & OMEGA) {
-        if (eflags[i] & SPHERE) {
-          radone = radius[i];
-          acm = b->angmom;
-          acm[0] += SINERTIA*rmass[i] * radone*radone * omega[i][0];
-          acm[1] += SINERTIA*rmass[i] * radone*radone * omega[i][1];
-          acm[2] += SINERTIA*rmass[i] * radone*radone * omega[i][2];
-        } else if (eflags[i] & LINE) {
-          radone = lbonus[line[i]].length;
-          b->angmom[2] += LINERTIA*rmass[i] * radone*radone * omega[i][2];
-        }
-      }
-      if (eflags[i] & ANGMOM) {
-        acm = b->angmom;
-        acm[0] += angmom[i][0];
-        acm[1] += angmom[i][1];
-        acm[2] += angmom[i][2];
-      }
       if (eflags[i] & TORQUE) {
         tcm = b->torque;
         tcm[0] += torque[i][0];
@@ -462,28 +421,17 @@ void FixRigidSmall::setup(int vflag)
     }
   }
 
-  // reverse communicate vcm, fcm, angmom, torque of all bodies
+  // reverse communicate fcm, torque of all bodies
 
-  commflag = VCM_ANGMOM;
-  comm->reverse_comm_variable_fix(this);
   commflag = FORCE_TORQUE;
   comm->reverse_comm_variable_fix(this);
 
-  // normalize velocity of COM
-
-  for (ibody = 0; ibody < nlocal_body; ibody++) {
-    vcm = body[ibody].vcm;
-    vcm[0] /= body[ibody].mass;
-    vcm[1] /= body[ibody].mass;
-    vcm[2] /= body[ibody].mass;
-  }
-
   // virial setup before call to set_v
 
   if (vflag) v_setup(vflag);
   else evflag = 0;
 
-  // compute and forward communicate updated info of all bodies
+  // compute and forward communicate vcm and omega of all bodies
 
   for (ibody = 0; ibody < nlocal_body; ibody++) {
     Body *b = &body[ibody];
@@ -855,7 +803,7 @@ int FixRigidSmall::dof(int tgroup)
 {
   int i,j;
 
-  // cannot count DOF correctly unless setup_bodies() has been called
+  // cannot count DOF correctly unless setup_bodies_static() has been called
 
   if (firstflag) {
     if (comm->me == 0) 
@@ -1559,7 +1507,7 @@ void FixRigidSmall::ring_farthest(int n, char *cbuf)
    read per-body attributes from infile if specified
 ------------------------------------------------------------------------- */
 
-void FixRigidSmall::setup_bodies()
+void FixRigidSmall::setup_bodies_static()
 {
   int i,itype,ibody;
 
@@ -2047,6 +1995,116 @@ void FixRigidSmall::setup_bodies()
   if (infile) memory->destroy(inbody);
 }
 
+/* ----------------------------------------------------------------------
+   one-time initialization of dynamic rigid body attributes
+   Vcm and angmom, computed explicitly from constituent particles
+   even if wrong for overlapping particles, is OK,
+     since is just setting initial time=0 Vcm and angmom of the body
+     which can be estimated value
+------------------------------------------------------------------------- */
+
+void FixRigidSmall::setup_bodies_dynamic()
+{
+  int i,n,ibody;
+  double massone,radone;
+
+  // sum vcm, angmom across all rigid bodies
+  // vcm = velocity of COM
+  // angmom = angular momentum around COM
+
+  double **x = atom->x;
+  double **v = atom->v;
+  double *rmass = atom->rmass;
+  double *mass = atom->mass;
+  int *type = atom->type;
+  int *image = atom->image;
+  int nlocal = atom->nlocal;
+
+  double *xcm,*vcm,*acm;
+  double dx,dy,dz;
+  double unwrap[3];
+
+  for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++) {
+    vcm = body[ibody].vcm;
+    vcm[0] = vcm[1] = vcm[2] = 0.0;
+    acm = body[ibody].angmom;
+    acm[0] = acm[1] = acm[2] = 0.0;
+  }
+
+  for (i = 0; i < nlocal; i++) {
+    if (atom2body[i] < 0) continue;
+    Body *b = &body[atom2body[i]];
+
+    if (rmass) massone = rmass[i];
+    else massone = mass[type[i]];
+
+    vcm = b->vcm;
+    vcm[0] += v[i][0] * massone;
+    vcm[1] += v[i][1] * massone;
+    vcm[2] += v[i][2] * massone;
+
+    domain->unmap(x[i],image[i],unwrap);
+    xcm = b->xcm;
+    dx = unwrap[0] - xcm[0];
+    dy = unwrap[1] - xcm[1];
+    dz = unwrap[2] - xcm[2];
+
+    acm = b->angmom;
+    acm[0] += dy * massone*v[i][2] - dz * massone*v[i][1];
+    acm[1] += dz * massone*v[i][0] - dx * massone*v[i][2];
+    acm[2] += dx * massone*v[i][1] - dy * massone*v[i][0];
+  }
+
+  // extended particles add their rotation to angmom of body
+
+  if (extended) {
+    AtomVecLine::Bonus *lbonus;
+    if (avec_line) lbonus = avec_line->bonus;
+    double **omega = atom->omega;
+    double **angmom = atom->angmom;
+    double *radius = atom->radius;
+    int *line = atom->line;
+
+    for (i = 0; i < nlocal; i++) {
+      if (atom2body[i] < 0) continue;
+      Body *b = &body[atom2body[i]];
+
+      if (eflags[i] & OMEGA) {
+        if (eflags[i] & SPHERE) {
+          radone = radius[i];
+          acm = b->angmom;
+          acm[0] += SINERTIA*rmass[i] * radone*radone * omega[i][0];
+          acm[1] += SINERTIA*rmass[i] * radone*radone * omega[i][1];
+          acm[2] += SINERTIA*rmass[i] * radone*radone * omega[i][2];
+        } else if (eflags[i] & LINE) {
+          radone = lbonus[line[i]].length;
+          b->angmom[2] += LINERTIA*rmass[i] * radone*radone * omega[i][2];
+        }
+      }
+      if (eflags[i] & ANGMOM) {
+        acm = b->angmom;
+        acm[0] += angmom[i][0];
+        acm[1] += angmom[i][1];
+        acm[2] += angmom[i][2];
+      }
+    }
+  }
+
+  // reverse communicate vcm, angmom of all bodies
+
+  commflag = VCM_ANGMOM;
+  comm->reverse_comm_variable_fix(this);
+
+  // normalize velocity of COM
+
+  for (ibody = 0; ibody < nlocal_body; ibody++) {
+    vcm = body[ibody].vcm;
+    vcm[0] /= body[ibody].mass;
+    vcm[1] /= body[ibody].mass;
+    vcm[2] /= body[ibody].mass;
+  }
+}
+
 /* ----------------------------------------------------------------------
    read per rigid body info from user-provided file
    which = 0 to read total mass and center-of-mass
@@ -2800,234 +2858,52 @@ void FixRigidSmall::reset_dt()
 
 /* ----------------------------------------------------------------------
    zero linear momentum of each rigid body
-   called by velocity zero linear command with its vgroup
-   only atoms in velocity command group contribute to vcm of body,
-     and only their velocities are adjusted, odd if partial bodies specified
+   set Vcm to 0.0, then reset velocities of particles via set_v()
 ------------------------------------------------------------------------- */
 
-void FixRigidSmall::zero_momentum(int vgroup)
+void FixRigidSmall::zero_momentum()
 {
-  int i,ibody;
-  double massone;
-
-  int vgroupbit = group->bitmask[vgroup];
-
-  // sum vcm across all rigid bodies
-
-  double **v = atom->v;
-  double *rmass = atom->rmass;
-  double *mass = atom->mass;
-  int *mask = atom->mask;
-  int *type = atom->type;
-  int nlocal = atom->nlocal;
-
-  double *vcm,*acm;
-
-  for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++) {
+  double *vcm;
+  for (int ibody = 0; ibody < nlocal_body+nghost_body; ibody++) {
     vcm = body[ibody].vcm;
     vcm[0] = vcm[1] = vcm[2] = 0.0;
-    acm = body[ibody].angmom;
-    acm[0] = acm[1] = acm[2] = 0.0;
   }
 
-  for (i = 0; i < nlocal; i++) {
-    if (!(mask[i] & vgroupbit)) continue;
-    if (atom2body[i] < 0) continue;
-    Body *b = &body[atom2body[i]];
-    if (rmass) massone = rmass[i];
-    else massone = mass[type[i]];
-    vcm = b->vcm;
-    vcm[0] += v[i][0] * massone;
-    vcm[1] += v[i][1] * massone;
-    vcm[2] += v[i][2] * massone;
-  }
-
-  // reverse communicate vcm (and angmom) of all bodies
-
-  commflag = VCM_ANGMOM;
-  comm->reverse_comm_variable_fix(this);
-
-  // normalize velocity of COM
-
-  for (ibody = 0; ibody < nlocal_body; ibody++) {
-    vcm = body[ibody].vcm;
-    vcm[0] /= body[ibody].mass;
-    vcm[1] /= body[ibody].mass;
-    vcm[2] /= body[ibody].mass;
-  }
-
-  // adjust velocities by vcm to zero linear momentum
-
-  for (i = 0; i < nlocal; i++) {
-    if (!(mask[i] & vgroupbit)) continue;
-    if (atom2body[i] < 0) continue;
-    Body *b = &body[atom2body[i]];
-    vcm = b->vcm;
-    v[i][0] -= vcm[0];
-    v[i][1] -= vcm[1];
-    v[i][2] -= vcm[2];
-  }
-}
-
-/* ----------------------------------------------------------------------
-   zero angular momentum of each rigid body
-   called by velocity zero linear command with its vgroup
-   only atoms in velocity command group contribute to angmom/omega of body,
-     and only their velocities are adjusted, odd if partial bodies specified
-------------------------------------------------------------------------- */
-
-void FixRigidSmall::zero_rotation(int vgroup)
-{
-  int i,ibody;
-  double massone,radone;
-
-  int vgroupbit = group->bitmask[vgroup];
-
-  // sum angmom across all rigid bodies
-
-  double **x = atom->x;
-  double **v = atom->v;
-  double *rmass = atom->rmass;
-  double *mass = atom->mass;
-  int *type = atom->type;
-  int *image = atom->image;
-  int *mask = atom->mask;
-  int nlocal = atom->nlocal;
-
-  double *xcm,*vcm,*acm;
-  double dx,dy,dz;
-  double unwrap[3];
-
-  for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++) {
-    vcm = body[ibody].vcm;
-    vcm[0] = vcm[1] = vcm[2] = 0.0;
-    acm = body[ibody].angmom;
-    acm[0] = acm[1] = acm[2] = 0.0;
-  }
-
-  for (i = 0; i < nlocal; i++) {
-    if (!(mask[i] & vgroupbit)) continue;
-    if (atom2body[i] < 0) continue;
-    Body *b = &body[atom2body[i]];
-
-    if (rmass) massone = rmass[i];
-    else massone = mass[type[i]];
-
-    domain->unmap(x[i],image[i],unwrap);
-    xcm = b->xcm;
-    dx = unwrap[0] - xcm[0];
-    dy = unwrap[1] - xcm[1];
-    dz = unwrap[2] - xcm[2];
-
-    acm = b->angmom;
-    acm[0] += dy * massone*v[i][2] - dz * massone*v[i][1];
-    acm[1] += dz * massone*v[i][0] - dx * massone*v[i][2];
-    acm[2] += dx * massone*v[i][1] - dy * massone*v[i][0];
-  }
-
-  // extended particles add their rotation to angmom of body
-
-  if (extended) {
-    AtomVecLine::Bonus *lbonus;
-    if (avec_line) lbonus = avec_line->bonus;
-    double **omega = atom->omega;
-    double **angmom = atom->angmom;
-    double **torque = atom->torque;
-    double *radius = atom->radius;
-    int *line = atom->line;
-
-    for (i = 0; i < nlocal; i++) {
-      if (!(mask[i] & vgroupbit)) continue;
-      if (atom2body[i] < 0) continue;
-      Body *b = &body[atom2body[i]];
-
-      if (eflags[i] & OMEGA) {
-        if (eflags[i] & SPHERE) {
-          radone = radius[i];
-          acm = b->angmom;
-          acm[0] += SINERTIA*rmass[i] * radone*radone * omega[i][0];
-          acm[1] += SINERTIA*rmass[i] * radone*radone * omega[i][1];
-          acm[2] += SINERTIA*rmass[i] * radone*radone * omega[i][2];
-        } else if (eflags[i] & LINE) {
-          radone = lbonus[line[i]].length;
-          b->angmom[2] += LINERTIA*rmass[i] * radone*radone * omega[i][2];
-        }
-      }
-      if (eflags[i] & ANGMOM) {
-        acm = b->angmom;
-        acm[0] += angmom[i][0];
-        acm[1] += angmom[i][1];
-        acm[2] += angmom[i][2];
-      }
-    }
-  }
-
-  // reverse communicate angmom (and vcm) of all bodies
-
-  commflag = VCM_ANGMOM;
-  comm->reverse_comm_variable_fix(this);
-
-  // convert angmom to omega and communicate it
-
-  for (ibody = 0; ibody < nlocal_body; ibody++) {
-    Body *b = &body[ibody];
-    MathExtra::angmom_to_omega(b->angmom,b->ex_space,b->ey_space,
-                               b->ez_space,b->inertia,b->omega);
-  }
+  // forward communicate of vcm to all ghost copies
 
   commflag = FINAL;
   comm->forward_comm_variable_fix(this);
 
-  // adjust velocities to zero omega
-  // vnew_i = v_i - w x r_i
-  // must use unwrapped coords to compute r_i correctly
+  // set velocity of atoms in rigid bodues
 
-  double *omega;
+  evflag = 0;
+  set_v();
+}
 
-  for (i = 0; i < nlocal; i++) {
-    if (!(mask[i] & vgroupbit)) continue;
-    if (atom2body[i] < 0) continue;
-    Body *b = &body[atom2body[i]];
+/* ----------------------------------------------------------------------
+   zero angular momentum of each rigid body
+   set angmom/omega to 0.0, then reset velocities of particles via set_v()
+------------------------------------------------------------------------- */
 
-    if (rmass) massone = rmass[i];
-    else massone = mass[type[i]];
-
-    domain->unmap(x[i],image[i],unwrap);
-    xcm = b->xcm;
-    dx = unwrap[0] - xcm[0];
-    dy = unwrap[1] - xcm[1];
-    dz = unwrap[2] - xcm[2];
-
-    omega = b->omega;
-    v[i][0] -= omega[1]*dz - omega[2]*dy;
-    v[i][1] -= omega[2]*dx - omega[0]*dz;
-    v[i][2] -= omega[0]*dy - omega[1]*dx;
+void FixRigidSmall::zero_rotation()
+{
+  double *angmom,*omega;
+  for (int ibody = 0; ibody < nlocal_body+nghost_body; ibody++) {
+    angmom = body[ibody].angmom;
+    angmom[0] = angmom[1] = angmom[2] = 0.0;
+    omega = body[ibody].omega;
+    omega[0] = omega[1] = omega[2] = 0.0;
   }
 
-  // also explicitly zero omega and angmom of extended particles
+  // forward communicate of omega to all ghost copies
 
-  if (extended) {
-    double **omega_one = atom->omega;
-    double **angmom_one = atom->angmom;
+  commflag = FINAL;
+  comm->forward_comm_variable_fix(this);
 
-    for (i = 0; i < nlocal; i++) {
-      if (!(mask[i] & vgroupbit)) continue;
-      if (atom2body[i] < 0) continue;
+  // set velocity of atoms in rigid bodues
 
-      if (eflags[i] & OMEGA) {
-        if (eflags[i] & SPHERE) {
-          omega_one[i][0] = 0.0;
-          omega_one[i][1] = 0.0;
-          omega_one[i][2] = 0.0;
-        } else if (eflags[i] & LINE) omega_one[i][2] = 0.0;
-      }
-      if (eflags[i] & ANGMOM) {
-        angmom_one[i][0] = 0.0;
-        angmom_one[i][1] = 0.0;
-        angmom_one[i][2] = 0.0;
-      }
-    }
-  }
+  evflag = 0;
+  set_v();
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/RIGID/fix_rigid_small.h b/src/RIGID/fix_rigid_small.h
index 650eb34306..0bdb686e9b 100644
--- a/src/RIGID/fix_rigid_small.h
+++ b/src/RIGID/fix_rigid_small.h
@@ -60,8 +60,8 @@ class FixRigidSmall : public Fix {
   int dof(int);
   void deform(int);
   void reset_dt();
-  void zero_momentum(int);
-  void zero_rotation(int);
+  void zero_momentum();
+  void zero_rotation();
   void *extract(const char*, int &);
   double extract_ke();
   double extract_erotational();
@@ -161,7 +161,8 @@ class FixRigidSmall : public Fix {
   void set_xv();
   void set_v();
   void create_bodies();
-  void setup_bodies();
+  void setup_bodies_static();
+  void setup_bodies_dynamic();
   void readfile(int, double **, int *);
   void grow_body();
   void reset_atom2body();
diff --git a/src/velocity.cpp b/src/velocity.cpp
index 501a128c6b..484809d65d 100644
--- a/src/velocity.cpp
+++ b/src/velocity.cpp
@@ -635,9 +635,9 @@ void Velocity::zero(int narg, char **arg)
       if (strcmp(modify->fix[rfix]->style,"rigid/small") == 0) {
         lmp->init();
         ((FixRigidSmall *) modify->fix[rfix])->setup_pre_neighbor();
-        ((FixRigidSmall *) modify->fix[rfix])->zero_momentum(igroup);
+        ((FixRigidSmall *) modify->fix[rfix])->zero_momentum();
       } else if (strstr(modify->fix[rfix]->style,"rigid")) {
-        ((FixRigid *) modify->fix[rfix])->zero_momentum(igroup);
+        ((FixRigid *) modify->fix[rfix])->zero_momentum();
       } else error->all(FLERR,"Velocity rigid used with non-rigid fix-ID");
     }
 
@@ -647,9 +647,9 @@ void Velocity::zero(int narg, char **arg)
       if (strcmp(modify->fix[rfix]->style,"rigid/small") == 0) {
         lmp->init();
         ((FixRigidSmall *) modify->fix[rfix])->setup_pre_neighbor();
-        ((FixRigidSmall *) modify->fix[rfix])->zero_rotation(igroup);
+        ((FixRigidSmall *) modify->fix[rfix])->zero_rotation();
       } else if (strstr(modify->fix[rfix]->style,"rigid")) {
-        ((FixRigid *) modify->fix[rfix])->zero_rotation(igroup);
+        ((FixRigid *) modify->fix[rfix])->zero_rotation();
       } else error->all(FLERR,"Velocity rigid used with non-rigid fix-ID");
     }