diff --git a/src/ASPHERE/compute_temp_asphere.cpp b/src/ASPHERE/compute_temp_asphere.cpp
index 3dd2bbb281..48b443f049 100755
--- a/src/ASPHERE/compute_temp_asphere.cpp
+++ b/src/ASPHERE/compute_temp_asphere.cpp
@@ -179,6 +179,8 @@ void ComputeTempAsphere::dof_compute()
   }
 
   dof -= extra_dof + fix_dof;
+  if (dof < 0.0 && natoms > 0.0) 
+    error->all(FLERR,"Temperature compute degrees of freedom < 0");
   if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
   else tfactor = 0.0;
 }
@@ -267,8 +269,6 @@ double ComputeTempAsphere::compute_scalar()
 
   MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
   if (dynamic || tempbias == 2) dof_compute();
-  if (tfactor == 0.0 && atom->natoms != 0) 
-    error->all(FLERR,"Temperature compute degrees of freedom < 0");
   scalar *= tfactor;
   return scalar;
 }
diff --git a/src/CORESHELL/compute_temp_cs.cpp b/src/CORESHELL/compute_temp_cs.cpp
index 73edbd1313..807a132e94 100644
--- a/src/CORESHELL/compute_temp_cs.cpp
+++ b/src/CORESHELL/compute_temp_cs.cpp
@@ -214,6 +214,8 @@ void ComputeTempCS::dof_compute()
   dof = nper * natoms;
   dof -= nper * nshells;
   dof -= extra_dof + fix_dof;
+  if (dof < 0.0 && natoms > 0.0) 
+    error->all(FLERR,"Temperature compute degrees of freedom < 0");
   if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
   else tfactor = 0.0;
 }
@@ -256,8 +258,6 @@ double ComputeTempCS::compute_scalar()
 
   MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
   if (dynamic) dof_compute();
-  if (tfactor == 0.0 && atom->natoms != 0) 
-    error->all(FLERR,"Temperature compute degrees of freedom < 0");
   scalar *= tfactor;
   return scalar;
 }
diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp
index 8b5697a810..e4cdaaca1c 100644
--- a/src/RIGID/fix_rigid_small.cpp
+++ b/src/RIGID/fix_rigid_small.cpp
@@ -508,7 +508,10 @@ void FixRigidSmall::init()
 /* ----------------------------------------------------------------------
    setup static/dynamic properties of rigid bodies, using current atom info
    allows resetting of atom properties like mass between runs
-   only do static initialization once if using an infile
+   only do static initialization must be done once
+     do not do again if using an infile or mol keyword
+     this is b/c the infile and Molecule class define the COM, etc
+       so should not reset those values
    cannot do this until now, b/c requires comm->setup() to have setup stencil
    invoke pre_neighbor() to insure body xcmimage flags are reset
      needed if Verlet::setup::pbc() has remapped/migrated atoms for 2nd run
@@ -517,9 +520,14 @@ void FixRigidSmall::init()
 
 void FixRigidSmall::setup_pre_neighbor()
 {
-  if (!staticflag || !infile) setup_bodies_static();
+  if (!staticflag || (!infile && !onemols)) setup_bodies_static();
   else pre_neighbor();
-  setup_bodies_dynamic();
+
+  //if (!staticflag || (!infile)) setup_bodies_static();
+  //else pre_neighbor();
+
+  //setup_bodies_dynamic();
+  if (!onemols) setup_bodies_dynamic();
 }
 
 /* ----------------------------------------------------------------------
@@ -880,24 +888,27 @@ void FixRigidSmall::final_integrate_respa(int ilevel, int iloop)
      due to first-time definition of rigid body in setup_bodies_static()
      or due to box flip
    also adjust imagebody = rigid body image flags, due to xcm remap
-   also adjust rgdimage flags of all atoms in bodies for two effects
+   then communicate bodies so other procs will know of changes to body xcm
+   then adjust xcmimage flags of all atoms in bodies via image_shift()
+     for two effects
      (1) change in true image flags due to pbc() call during exchange
      (2) change in imagebody due to xcm remap
-   rgdimage flags are always -1,0,-1 so that body can be unwrapped
+   xcmimage flags are always -1,0,-1 so that body can be unwrapped
      around in-box xcm and stay close to simulation box
-   if don't do this, then a body could end up very far away
-     when unwrapped by true image flags, and set_xv() will
-     compute huge displacements every step to reset coords of 
+   if just inferred unwrapped from atom image flags, 
+     then a body could end up very far away
+     when unwrapped by true image flags
+   then set_xv() will compute huge displacements every step to reset coords of 
      all the body atoms to be back inside the box, ditto for triclinic box flip
      note: so just want to avoid that numeric probem?
 ------------------------------------------------------------------------- */
 
 void FixRigidSmall::pre_neighbor()
 {
-  // acquire ghost bodies via forward comm
-  // also gets new remapflags needed for adjusting atom image flags
-  // reset atom2body for owned atoms
-  // forward comm sets atom2body for ghost atoms
+  for (int ibody = 0; ibody < nlocal_body; ibody++) {
+    Body *b = &body[ibody];
+    domain->remap(b->xcm,b->image);
+  }  
 
   nghost_body = 0;
   commflag = FULL_BODY;
@@ -905,10 +916,6 @@ void FixRigidSmall::pre_neighbor()
   reset_atom2body();
   //check(4);
 
-  for (int ibody = 0; ibody < nlocal_body; ibody++) {
-    Body *b = &body[ibody];
-    domain->remap(b->xcm,b->image);
-  }  
   image_shift();
 }
 
@@ -1654,7 +1661,8 @@ void FixRigidSmall::ring_farthest(int n, char *cbuf)
 /* ----------------------------------------------------------------------
    initialization of rigid body attributes
    called at setup, so body/atom properties can be changed between runs
-   unless reading from infile, in which case only called once
+   unless reading from infile or defined by Molecule class,
+     in which case only called once
    sets extended flags, masstotal, center-of-mass
    sets Cartesian and diagonalized inertia tensor
    sets body image flags, but only on first call
@@ -2180,7 +2188,7 @@ void FixRigidSmall::setup_bodies_static()
    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
+     since is just setting initial time=0 vcm and angmom of the body
      which can be estimated value
 ------------------------------------------------------------------------- */
 
diff --git a/src/RIGID/fix_rigid_small.h b/src/RIGID/fix_rigid_small.h
index 06449da275..5dacdbf9ac 100644
--- a/src/RIGID/fix_rigid_small.h
+++ b/src/RIGID/fix_rigid_small.h
@@ -255,7 +255,7 @@ NPT/NPH fix must be defined in input script after all rigid fixes,
 else the rigid fix contribution to the pressure virial is
 incorrect.
 
-W: Cannot count rigid body degrees-of-freedom before bodies are fully initialized
+W: Cannot count rigid body degrees-of-freedom before bodies are fully initialized h
 
 This means the temperature associated with the rigid bodies may be
 incorrect on this timestep.
diff --git a/src/USER-EFF/compute_temp_deform_eff.cpp b/src/USER-EFF/compute_temp_deform_eff.cpp
index 7472b4d230..92ebdb298a 100644
--- a/src/USER-EFF/compute_temp_deform_eff.cpp
+++ b/src/USER-EFF/compute_temp_deform_eff.cpp
@@ -123,6 +123,8 @@ void ComputeTempDeformEff::dof_compute()
 
   dof -= domain->dimension * nelectrons;
 
+  if (dof < 0.0 && natoms > 0.0) 
+    error->all(FLERR,"Temperature compute degrees of freedom < 0");
   if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
   else tfactor = 0.0;
 }
@@ -174,8 +176,6 @@ double ComputeTempDeformEff::compute_scalar()
 
   MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
   if (dynamic) dof_compute();
-  if (tfactor == 0.0 && atom->natoms != 0) 
-    error->all(FLERR,"Temperature compute degrees of freedom < 0");
   scalar *= tfactor;
   return scalar;
 }
diff --git a/src/USER-EFF/compute_temp_eff.cpp b/src/USER-EFF/compute_temp_eff.cpp
index 9f1bafdd49..9b24e0a707 100644
--- a/src/USER-EFF/compute_temp_eff.cpp
+++ b/src/USER-EFF/compute_temp_eff.cpp
@@ -87,6 +87,8 @@ void ComputeTempEff::dof_compute()
 
   dof -= domain->dimension * nelectrons;
 
+  if (dof < 0.0 && natoms > 0.0) 
+    error->all(FLERR,"Temperature compute degrees of freedom < 0");
   if (dof > 0.0) tfactor = force->mvv2e / (dof * force->boltz);
   else tfactor = 0.0;
 }
@@ -120,8 +122,6 @@ double ComputeTempEff::compute_scalar()
 
   MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
   if (dynamic) dof_compute();
-  if (tfactor == 0.0 && atom->natoms != 0) 
-    error->all(FLERR,"Temperature compute degrees of freedom < 0");
   scalar *= tfactor;
   return scalar;
 }
diff --git a/src/compute_temp.cpp b/src/compute_temp.cpp
index 4ad653b2f9..5a472549ab 100644
--- a/src/compute_temp.cpp
+++ b/src/compute_temp.cpp
@@ -64,6 +64,8 @@ void ComputeTemp::dof_compute()
   double natoms = group->count(igroup);
   dof = domain->dimension * natoms;
   dof -= extra_dof + fix_dof;
+  if (dof < 0.0 && natoms > 0.0) 
+    error->all(FLERR,"Temperature compute degrees of freedom < 0");
   if (dof > 0.0) tfactor = force->mvv2e / (dof * force->boltz);
   else tfactor = 0.0;
 }
@@ -96,8 +98,6 @@ double ComputeTemp::compute_scalar()
 
   MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
   if (dynamic) dof_compute();
-  if (tfactor == 0.0 && atom->natoms != 0) 
-    error->all(FLERR,"Temperature compute degrees of freedom < 0");
   scalar *= tfactor;
   return scalar;
 }
diff --git a/src/compute_temp_chunk.cpp b/src/compute_temp_chunk.cpp
index 2d36bdf0b7..afbe02129a 100644
--- a/src/compute_temp_chunk.cpp
+++ b/src/compute_temp_chunk.cpp
@@ -288,6 +288,8 @@ double ComputeTempChunk::compute_scalar()
   MPI_Allreduce(&rcount,&allcount,1,MPI_DOUBLE,MPI_SUM,world);
 
   double dof = nchunk*cdof + adof*allcount;
+  if (dof < 0.0 && allcount > 0.0)
+    error->all(FLERR,"Temperature compute degrees of freedom < 0");
   double tfactor = 0.0;
   if (dof > 0.0) tfactor = force->mvv2e / (dof * force->boltz);
   scalar *= tfactor;
diff --git a/src/compute_temp_com.cpp b/src/compute_temp_com.cpp
index 8e0db3595d..e9f31dbfe1 100644
--- a/src/compute_temp_com.cpp
+++ b/src/compute_temp_com.cpp
@@ -71,9 +71,10 @@ void ComputeTempCOM::dof_compute()
 {
   adjust_dof_fix();
   double natoms = group->count(igroup);
-  int nper = domain->dimension;
-  dof = nper * natoms;
+  dof = domain->dimension * natoms;
   dof -= extra_dof + fix_dof;
+  if (dof < 0.0 && natoms > 0.0) 
+    error->all(FLERR,"Temperature compute degrees of freedom < 0");
   if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
   else tfactor = 0.0;
 }
@@ -112,8 +113,6 @@ double ComputeTempCOM::compute_scalar()
 
   MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
   if (dynamic) dof_compute();
-  if (tfactor == 0.0 && atom->natoms != 0) 
-    error->all(FLERR,"Temperature compute degrees of freedom < 0");
   scalar *= tfactor;
   return scalar;
 }
diff --git a/src/compute_temp_deform.cpp b/src/compute_temp_deform.cpp
index e608ba3cea..7513f6fb1c 100644
--- a/src/compute_temp_deform.cpp
+++ b/src/compute_temp_deform.cpp
@@ -99,6 +99,8 @@ void ComputeTempDeform::dof_compute()
   double natoms = group->count(igroup);
   dof = domain->dimension * natoms;
   dof -= extra_dof + fix_dof;
+  if (dof < 0.0 && natoms > 0.0) 
+    error->all(FLERR,"Temperature compute degrees of freedom < 0");
   if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
   else tfactor = 0.0;
 }
@@ -148,8 +150,6 @@ double ComputeTempDeform::compute_scalar()
 
   MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
   if (dynamic) dof_compute();
-  if (tfactor == 0.0 && atom->natoms != 0) 
-    error->all(FLERR,"Temperature compute degrees of freedom < 0");
   scalar *= tfactor;
   return scalar;
 }
diff --git a/src/compute_temp_profile.cpp b/src/compute_temp_profile.cpp
index 2b13dd0927..9f0ff3f31e 100644
--- a/src/compute_temp_profile.cpp
+++ b/src/compute_temp_profile.cpp
@@ -196,9 +196,10 @@ void ComputeTempProfile::dof_compute()
 {
   adjust_dof_fix();
   double natoms = group->count(igroup);
-  int nper = domain->dimension;
-  dof = nper * natoms;
+  dof = domain->dimension * natoms;
   dof -= extra_dof + fix_dof;
+  if (dof < 0.0 && natoms > 0.0) 
+    error->all(FLERR,"Temperature compute degrees of freedom < 0");
   if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
   else tfactor = 0.0;
 }
@@ -242,8 +243,6 @@ double ComputeTempProfile::compute_scalar()
 
   MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
   if (dynamic) dof_compute();
-  if (tfactor == 0.0 && atom->natoms != 0) 
-    error->all(FLERR,"Temperature compute degrees of freedom < 0");
   scalar *= tfactor;
   return scalar;
 }
diff --git a/src/compute_temp_ramp.cpp b/src/compute_temp_ramp.cpp
index 76b47a39a2..0c38ed811a 100644
--- a/src/compute_temp_ramp.cpp
+++ b/src/compute_temp_ramp.cpp
@@ -127,9 +127,10 @@ void ComputeTempRamp::dof_compute()
 {
   adjust_dof_fix();
   double natoms = group->count(igroup);
-  int nper = domain->dimension;
-  dof = nper * natoms;
+  dof = domain->dimension * natoms;
   dof -= extra_dof + fix_dof;
+  if (dof < 0.0 && natoms > 0.0) 
+    error->all(FLERR,"Temperature compute degrees of freedom < 0");
   if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
   else tfactor = 0.0;
 }
@@ -171,8 +172,6 @@ double ComputeTempRamp::compute_scalar()
 
   MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
   if (dynamic) dof_compute();
-  if (tfactor == 0.0 && atom->natoms != 0) 
-    error->all(FLERR,"Temperature compute degrees of freedom < 0");
   scalar *= tfactor;
   return scalar;
 }
diff --git a/src/compute_temp_sphere.cpp b/src/compute_temp_sphere.cpp
index 70f13d40b0..050a987933 100644
--- a/src/compute_temp_sphere.cpp
+++ b/src/compute_temp_sphere.cpp
@@ -125,6 +125,7 @@ void ComputeTempSphere::dof_compute()
   int count,count_all;
 
   adjust_dof_fix();
+  double natoms = group->count(igroup);
 
   // 6 or 3 dof for extended/point particles for 3d
   // 3 or 2 dof for extended/point particles for 2d
@@ -165,10 +166,7 @@ void ComputeTempSphere::dof_compute()
   // additional adjustments to dof
 
   if (tempbias == 1) {
-    if (mode == ALL) {
-      double natoms = group->count(igroup);
-      dof -= tbias->dof_remove(-1) * natoms;
-    }
+    if (mode == ALL) dof -= tbias->dof_remove(-1) * natoms;
 
   } else if (tempbias == 2) {
     int *mask = atom->mask;
@@ -208,6 +206,8 @@ void ComputeTempSphere::dof_compute()
   }
 
   dof -= extra_dof + fix_dof;
+  if (dof < 0.0 && natoms > 0.0) 
+    error->all(FLERR,"Temperature compute degrees of freedom < 0");
   if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
   else tfactor = 0.0;
 }
@@ -252,8 +252,6 @@ double ComputeTempSphere::compute_scalar()
 
   MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
   if (dynamic || tempbias == 2) dof_compute();
-  if (tfactor == 0.0 && atom->natoms != 0)
-    error->all(FLERR,"Temperature compute degrees of freedom < 0");
   scalar *= tfactor;
   return scalar;
 }