diff --git a/src/ASPHERE/compute_temp_asphere.cpp b/src/ASPHERE/compute_temp_asphere.cpp
index 34b7a4aac3..c58862fc03 100755
--- a/src/ASPHERE/compute_temp_asphere.cpp
+++ b/src/ASPHERE/compute_temp_asphere.cpp
@@ -133,7 +133,6 @@ void ComputeTempAsphere::setup()
 {
   dynamic = 0;
   if (dynamic_user || group->dynamic[igroup]) dynamic = 1;
-  fix_dof = -1;
   dof_compute();
 }
 
@@ -141,7 +140,7 @@ void ComputeTempAsphere::setup()
 
 void ComputeTempAsphere::dof_compute()
 {
-  if (fix_dof) adjust_dof_fix();
+  fix_dof = modify->adjust_dof_fix(igroup);
 
   // 6 dof for 3d, 3 dof for 2d
   // which dof are included also depends on mode
@@ -268,6 +267,8 @@ double ComputeTempAsphere::compute_scalar()
 
   MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
   if (dynamic || tempbias == 2) dof_compute();
+  if (tfactor == 0.0 && scalar != 0.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 7512d48afd..2ac071e8b2 100644
--- a/src/CORESHELL/compute_temp_cs.cpp
+++ b/src/CORESHELL/compute_temp_cs.cpp
@@ -201,7 +201,6 @@ void ComputeTempCS::setup()
 
   // calculate DOF for temperature
 
-  fix_dof = -1;
   dof_compute();
 }
 
@@ -209,7 +208,7 @@ void ComputeTempCS::setup()
 
 void ComputeTempCS::dof_compute()
 {
-  if (fix_dof) adjust_dof_fix();
+  fix_dof = modify->adjust_dof_fix(igroup);
   int nper = domain->dimension;
   double natoms = group->count(igroup);
   dof = nper * natoms;
@@ -257,6 +256,8 @@ double ComputeTempCS::compute_scalar()
 
   MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
   if (dynamic) dof_compute();
+  if (tfactor == 0.0 && scalar != 0.0) 
+    error->all(FLERR,"Temperature compute degrees of freedom < 0");
   scalar *= tfactor;
   return scalar;
 }
diff --git a/src/MISC/fix_deposit.cpp b/src/MISC/fix_deposit.cpp
index 1cf78d2dcc..0b35627ebe 100644
--- a/src/MISC/fix_deposit.cpp
+++ b/src/MISC/fix_deposit.cpp
@@ -249,6 +249,41 @@ void FixDeposit::init()
       error->all(FLERR,"Fix deposit and fix shake not using "
                  "same molecule template ID");
   }
+
+  // for finite size spherical particles:
+  // warn if near < 2 * maxrad of existing and inserted particles
+  //   since may lead to overlaps
+  // if inserted molecule does not define diameters, 
+  //   use AtomVecSphere::create_atom() default radius = 0.5
+
+  if (atom->radius_flag) {
+    double *radius = atom->radius;
+    int nlocal = atom->nlocal;
+
+    double maxrad = 0.0;
+    for (int i = 0; i < nlocal; i++)
+      maxrad = MAX(maxrad,radius[i]);
+
+    double maxradall;
+    MPI_Allreduce(&maxrad,&maxradall,1,MPI_DOUBLE,MPI_MAX,world);
+
+    double maxradinsert = 0.0;
+    if (mode == MOLECULE) {
+      for (int i = 0; i < nmol; i++) {
+        if (onemols[i]->radiusflag)
+          maxradinsert = MAX(maxradinsert,onemols[i]->maxradius);
+        else maxradinsert = MAX(maxradinsert,0.5);
+      }
+    } else maxradinsert = 0.5;
+
+    double separation = MAX(2.0*maxradinsert,maxradall+maxradinsert);
+    if (sqrt(nearsq) < separation && comm->me == 0) {
+      char str[128];
+      sprintf(str,"Fix deposit near setting < possible overlap separation %g",
+              separation);
+      error->warning(FLERR,str);
+    }
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -389,7 +424,8 @@ void FixDeposit::pre_exchange()
       }
     }
 
-    // if distance to any inserted atom is less than near, try again
+    // check distance between any existing atom and any inserted atom
+    // if less than near, try again
     // use minimum_image() to account for PBC
 
     double **x = atom->x;
diff --git a/src/POEMS/fix_poems.cpp b/src/POEMS/fix_poems.cpp
index 7f07accc57..47e96f7192 100644
--- a/src/POEMS/fix_poems.cpp
+++ b/src/POEMS/fix_poems.cpp
@@ -71,6 +71,7 @@ FixPOEMS::FixPOEMS(LAMMPS *lmp, int narg, char **arg) :
   time_integrate = 1;
   rigid_flag = 1;
   virial_flag = 1;
+  dof_flag = 1;
 
   MPI_Comm_rank(world,&me);
 
diff --git a/src/RIGID/fix_rigid.cpp b/src/RIGID/fix_rigid.cpp
index 359e70f851..5cac49be92 100644
--- a/src/RIGID/fix_rigid.cpp
+++ b/src/RIGID/fix_rigid.cpp
@@ -66,6 +66,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
   rigid_flag = 1;
   virial_flag = 1;
   create_attribute = 1;
+  dof_flag = 1;
 
   MPI_Comm_rank(world,&me);
   MPI_Comm_size(world,&nprocs);
diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp
index cced1764f7..42f9089170 100644
--- a/src/RIGID/fix_rigid_small.cpp
+++ b/src/RIGID/fix_rigid_small.cpp
@@ -78,6 +78,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
   rigid_flag = 1;
   virial_flag = 1;
   create_attribute = 1;
+  dof_flag = 1;
 
   MPI_Comm_rank(world,&me);
   MPI_Comm_size(world,&nprocs);
@@ -2191,7 +2192,7 @@ void FixRigidSmall::setup_bodies_static()
 
 /* ----------------------------------------------------------------------
    one-time initialization of dynamic rigid body attributes
-   Vcm and angmom, computed explicitly from constituent particles
+   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
diff --git a/src/RIGID/fix_shake.cpp b/src/RIGID/fix_shake.cpp
index bf41b17a5b..bf0c6d625c 100644
--- a/src/RIGID/fix_shake.cpp
+++ b/src/RIGID/fix_shake.cpp
@@ -55,6 +55,7 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
 
   virial_flag = 1;
   create_attribute = 1;
+  dof_flag = 1;
 
   // error check
 
diff --git a/src/USER-CUDA/fix_shake_cuda.cpp b/src/USER-CUDA/fix_shake_cuda.cpp
index 79e3e82f67..070274cedd 100644
--- a/src/USER-CUDA/fix_shake_cuda.cpp
+++ b/src/USER-CUDA/fix_shake_cuda.cpp
@@ -63,6 +63,7 @@ FixShakeCuda::FixShakeCuda(LAMMPS* lmp, int narg, char** arg) :
 
   virial_flag = 1;
   create_attribute = 1;
+  dof_flag = 1;
 
   // error check
 
diff --git a/src/USER-EFF/compute_temp_deform_eff.cpp b/src/USER-EFF/compute_temp_deform_eff.cpp
index 30729097a4..cb0cf0f7f2 100644
--- a/src/USER-EFF/compute_temp_deform_eff.cpp
+++ b/src/USER-EFF/compute_temp_deform_eff.cpp
@@ -93,7 +93,6 @@ void ComputeTempDeformEff::setup()
 {
   dynamic = 0;
   if (dynamic_user || group->dynamic[igroup]) dynamic = 1;
-  fix_dof = -1;
   dof_compute();
 }
 
@@ -101,6 +100,7 @@ void ComputeTempDeformEff::setup()
 
 void ComputeTempDeformEff::dof_compute()
 {
+  fix_dof = modify->adjust_dof_fix(igroup);
   double natoms = group->count(igroup);
   dof = domain->dimension * natoms;
   dof -= extra_dof + fix_dof;
@@ -120,6 +120,7 @@ void ComputeTempDeformEff::dof_compute()
   MPI_Allreduce(&one,&nelectrons,1,MPI_INT,MPI_SUM,world);
 
   // Assume 3/2 k T per nucleus
+
   dof -= domain->dimension * nelectrons;
 
   if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
@@ -173,6 +174,8 @@ double ComputeTempDeformEff::compute_scalar()
 
   MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
   if (dynamic) dof_compute();
+  if (tfactor == 0.0 && scalar != 0.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 634bd83c06..3dcd0170db 100644
--- a/src/USER-EFF/compute_temp_eff.cpp
+++ b/src/USER-EFF/compute_temp_eff.cpp
@@ -59,7 +59,6 @@ void ComputeTempEff::setup()
 {
   dynamic = 0;
   if (dynamic_user || group->dynamic[igroup]) dynamic = 1;
-  fix_dof = -1;
   dof_compute();
 }
 
@@ -67,7 +66,7 @@ void ComputeTempEff::setup()
 
 void ComputeTempEff::dof_compute()
 {
-  if (fix_dof) adjust_dof_fix();
+  fix_dof = adjust_dof_fix(igroup);
   double natoms = group->count(igroup);
   dof = domain->dimension * natoms;
   dof -= extra_dof + fix_dof;
@@ -121,6 +120,8 @@ double ComputeTempEff::compute_scalar()
 
   MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
   if (dynamic) dof_compute();
+  if (tfactor == 0.0 && scalar != 0.0) 
+    error->all(FLERR,"Temperature compute degrees of freedom < 0");
   scalar *= tfactor;
   return scalar;
 }
diff --git a/src/USER-LB/fix_lb_rigid_pc_sphere.cpp b/src/USER-LB/fix_lb_rigid_pc_sphere.cpp
index c031d1b782..d97ba671b5 100755
--- a/src/USER-LB/fix_lb_rigid_pc_sphere.cpp
+++ b/src/USER-LB/fix_lb_rigid_pc_sphere.cpp
@@ -51,6 +51,7 @@ FixLbRigidPCSphere::FixLbRigidPCSphere(LAMMPS *lmp, int narg, char **arg) :
   rigid_flag = 1;
   create_attribute = 1;
   virial_flag = 1;
+  dof_flag = 1;
 
   // perform initial allocation of atom-based arrays
   // register with Atom class
diff --git a/src/USER-MISC/compute_temp_rotate.cpp b/src/USER-MISC/compute_temp_rotate.cpp
index 7c85fc1d15..17ffc0ce71 100644
--- a/src/USER-MISC/compute_temp_rotate.cpp
+++ b/src/USER-MISC/compute_temp_rotate.cpp
@@ -70,7 +70,6 @@ void ComputeTempRotate::setup()
 {
   dynamic = 0;
   if (dynamic_user || group->dynamic[igroup]) dynamic = 1;
-  fix_dof = -1;
   dof_compute();
 }
 
@@ -78,7 +77,7 @@ void ComputeTempRotate::setup()
 
 void ComputeTempRotate::dof_compute()
 {
-  if (fix_dof) adjust_dof_fix();
+  fix_dof = adjust_dof_fix(igroup);
   double natoms = group->count(igroup);
   int nper = domain->dimension;
   dof = nper * natoms;
@@ -144,6 +143,8 @@ double ComputeTempRotate::compute_scalar()
 
   MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
   if (dynamic) dof_compute();
+  if (tfactor == 0.0 && scalar != 0.0) 
+    error->all(FLERR,"Temperature compute degrees of freedom < 0");
   scalar *= tfactor;
   return scalar;
 }
diff --git a/src/compute.cpp b/src/compute.cpp
index 4beefac23c..026dbe4103 100644
--- a/src/compute.cpp
+++ b/src/compute.cpp
@@ -143,17 +143,6 @@ void Compute::modify_params(int narg, char **arg)
   }
 }
 
-/* ----------------------------------------------------------------------
-   calculate adjustment in DOF due to fixes
-------------------------------------------------------------------------- */
-
-void Compute::adjust_dof_fix()
-{
-  fix_dof = 0;
-  for (int i = 0; i < modify->nfix; i++)
-    fix_dof += modify->fix[i]->dof(igroup);
-}
-
 /* ----------------------------------------------------------------------
    reset extra_dof to its default value
 ------------------------------------------------------------------------- */
diff --git a/src/compute.h b/src/compute.h
index d914652be2..281010c1ee 100644
--- a/src/compute.h
+++ b/src/compute.h
@@ -92,7 +92,6 @@ class Compute : protected Pointers {
   Compute(class LAMMPS *, int, char **);
   virtual ~Compute();
   void modify_params(int, char **);
-  void adjust_dof_fix();
   void reset_extra_dof();
 
   virtual void init() = 0;
diff --git a/src/compute_temp.cpp b/src/compute_temp.cpp
index 6c0fb09c08..f5f28df80f 100644
--- a/src/compute_temp.cpp
+++ b/src/compute_temp.cpp
@@ -18,6 +18,8 @@
 #include "update.h"
 #include "force.h"
 #include "domain.h"
+#include "modify.h"
+#include "comm.h"
 #include "group.h"
 #include "error.h"
 
@@ -52,7 +54,6 @@ void ComputeTemp::setup()
 {
   dynamic = 0;
   if (dynamic_user || group->dynamic[igroup]) dynamic = 1;
-  fix_dof = -1;
   dof_compute();
 }
 
@@ -60,7 +61,7 @@ void ComputeTemp::setup()
 
 void ComputeTemp::dof_compute()
 {
-  if (fix_dof) adjust_dof_fix();
+  fix_dof = modify->adjust_dof_fix(igroup);
   double natoms = group->count(igroup);
   dof = domain->dimension * natoms;
   dof -= extra_dof + fix_dof;
@@ -96,6 +97,8 @@ double ComputeTemp::compute_scalar()
 
   MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
   if (dynamic) dof_compute();
+  if (tfactor == 0.0 && scalar != 0.0) 
+    error->all(FLERR,"Temperature compute degrees of freedom < 0");
   scalar *= tfactor;
   return scalar;
 }
diff --git a/src/compute_temp_com.cpp b/src/compute_temp_com.cpp
index 43591a92b2..8d8b3fe1fe 100644
--- a/src/compute_temp_com.cpp
+++ b/src/compute_temp_com.cpp
@@ -20,6 +20,7 @@
 #include "force.h"
 #include "group.h"
 #include "domain.h"
+#include "modify.h"
 #include "lattice.h"
 #include "error.h"
 
@@ -62,7 +63,6 @@ void ComputeTempCOM::setup()
 {
   dynamic = 0;
   if (dynamic_user || group->dynamic[igroup]) dynamic = 1;
-  fix_dof = -1;
   dof_compute();
 }
 
@@ -70,7 +70,7 @@ void ComputeTempCOM::setup()
 
 void ComputeTempCOM::dof_compute()
 {
-  if (fix_dof) adjust_dof_fix();
+  fix_dof = modify->adjust_dof_fix(igroup);
   double natoms = group->count(igroup);
   int nper = domain->dimension;
   dof = nper * natoms;
@@ -113,6 +113,8 @@ double ComputeTempCOM::compute_scalar()
 
   MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
   if (dynamic) dof_compute();
+  if (tfactor == 0.0 && scalar != 0.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 7acce83342..50b7dd7344 100644
--- a/src/compute_temp_deform.cpp
+++ b/src/compute_temp_deform.cpp
@@ -88,7 +88,6 @@ void ComputeTempDeform::setup()
 {
   dynamic = 0;
   if (dynamic_user || group->dynamic[igroup]) dynamic = 1;
-  fix_dof = -1;
   dof_compute();
 }
 
@@ -96,7 +95,7 @@ void ComputeTempDeform::setup()
 
 void ComputeTempDeform::dof_compute()
 {
-  if (fix_dof) adjust_dof_fix();
+  fix_dof = modify->adjust_dof_fix(igroup);
   double natoms = group->count(igroup);
   dof = domain->dimension * natoms;
   dof -= extra_dof + fix_dof;
@@ -149,6 +148,8 @@ double ComputeTempDeform::compute_scalar()
 
   MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
   if (dynamic) dof_compute();
+  if (tfactor == 0.0 && scalar != 0.0) 
+    error->all(FLERR,"Temperature compute degrees of freedom < 0");
   scalar *= tfactor;
   return scalar;
 }
diff --git a/src/compute_temp_partial.cpp b/src/compute_temp_partial.cpp
index a9468103ce..7fdd0c30d0 100644
--- a/src/compute_temp_partial.cpp
+++ b/src/compute_temp_partial.cpp
@@ -18,6 +18,7 @@
 #include "update.h"
 #include "force.h"
 #include "domain.h"
+#include "modify.h"
 #include "group.h"
 #include "memory.h"
 #include "error.h"
@@ -63,7 +64,6 @@ void ComputeTempPartial::setup()
 {
   dynamic = 0;
   if (dynamic_user || group->dynamic[igroup]) dynamic = 1;
-  fix_dof = -1;
   dof_compute();
 }
 
@@ -74,7 +74,7 @@ void ComputeTempPartial::setup()
 
 void ComputeTempPartial::dof_compute()
 {
-  if (fix_dof) adjust_dof_fix();
+  fix_dof = modify->adjust_dof_fix(igroup);
   double natoms = group->count(igroup);
   int nper = xflag+yflag+zflag;
   dof = nper * natoms;
@@ -120,6 +120,8 @@ double ComputeTempPartial::compute_scalar()
 
   MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
   if (dynamic) dof_compute();
+  if (tfactor == 0.0 && scalar != 0.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 58afbb3a5b..b6905946ab 100644
--- a/src/compute_temp_profile.cpp
+++ b/src/compute_temp_profile.cpp
@@ -161,7 +161,6 @@ ComputeTempProfile::~ComputeTempProfile()
 
 void ComputeTempProfile::init()
 {
-  fix_dof = -1;
   dof_compute();
 
   // ptrs to domain data
@@ -189,7 +188,6 @@ void ComputeTempProfile::setup()
 {
   dynamic = 0;
   if (dynamic_user || group->dynamic[igroup]) dynamic = 1;
-  fix_dof = -1;
   dof_compute();
 }
 
@@ -197,7 +195,7 @@ void ComputeTempProfile::setup()
 
 void ComputeTempProfile::dof_compute()
 {
-  if (fix_dof) adjust_dof_fix();
+  fix_dof = modify->adjust_dof_fix(igroup);
   double natoms = group->count(igroup);
   int nper = domain->dimension;
   dof = nper * natoms;
@@ -245,6 +243,8 @@ double ComputeTempProfile::compute_scalar()
 
   MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
   if (dynamic) dof_compute();
+  if (tfactor == 0.0 && scalar != 0.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 c4b1ee4dc1..04f714ad6f 100644
--- a/src/compute_temp_ramp.cpp
+++ b/src/compute_temp_ramp.cpp
@@ -119,7 +119,6 @@ void ComputeTempRamp::setup()
 {
   dynamic = 0;
   if (dynamic_user || group->dynamic[igroup]) dynamic = 1;
-  fix_dof = -1;
   dof_compute();
 }
 
@@ -127,7 +126,7 @@ void ComputeTempRamp::setup()
 
 void ComputeTempRamp::dof_compute()
 {
-  if (fix_dof) adjust_dof_fix();
+  fix_dof = modify->adjust_dof_fix(igroup);
   double natoms = group->count(igroup);
   int nper = domain->dimension;
   dof = nper * natoms;
@@ -173,6 +172,8 @@ double ComputeTempRamp::compute_scalar()
 
   MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
   if (dynamic) dof_compute();
+  if (tfactor == 0.0 && scalar != 0.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 b80836e27d..b7ca03837a 100644
--- a/src/compute_temp_sphere.cpp
+++ b/src/compute_temp_sphere.cpp
@@ -20,6 +20,7 @@
 #include "force.h"
 #include "domain.h"
 #include "modify.h"
+#include "comm.h"
 #include "group.h"
 #include "error.h"
 
@@ -114,7 +115,6 @@ void ComputeTempSphere::setup()
 {
   dynamic = 0;
   if (dynamic_user || group->dynamic[igroup]) dynamic = 1;
-  fix_dof = -1;
   dof_compute();
 }
 
@@ -124,7 +124,7 @@ void ComputeTempSphere::dof_compute()
 {
   int count,count_all;
 
-  if (fix_dof) adjust_dof_fix();
+  fix_dof = modify->adjust_dof_fix(igroup);
 
   // 6 or 3 dof for extended/point particles for 3d
   // 3 or 2 dof for extended/point particles for 2d
@@ -252,6 +252,8 @@ double ComputeTempSphere::compute_scalar()
 
   MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
   if (dynamic || tempbias == 2) dof_compute();
+  if (tfactor == 0.0 && scalar != 0.0)
+    error->all(FLERR,"Temperature compute degrees of freedom < 0");
   scalar *= tfactor;
   return scalar;
 }
diff --git a/src/fix.cpp b/src/fix.cpp
index e670b380cd..308997d4a0 100644
--- a/src/fix.cpp
+++ b/src/fix.cpp
@@ -65,6 +65,7 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
   restart_pbc = 0;
   wd_header = wd_section = 0;
   dynamic_group_allow = 0;
+  dof_flag = 0;
   cudable_comm = 0;
 
   scalar_flag = vector_flag = array_flag = 0;
diff --git a/src/fix.h b/src/fix.h
index 1b53ec2645..bd19b8a971 100644
--- a/src/fix.h
+++ b/src/fix.h
@@ -49,6 +49,7 @@ class Fix : protected Pointers {
   int wd_header;                 // # of header values fix writes to data file
   int wd_section;                // # of sections fix writes to data file
   int dynamic_group_allow;       // 1 if can be used with dynamic group, else 0
+  int dof_flag;                  // 1 if has dof() method (not min_dof())
   int cudable_comm;              // 1 if fix has CUDA-enabled communication
 
   int scalar_flag;               // 0/1 if compute_scalar() function exists
diff --git a/src/modify.cpp b/src/modify.cpp
index 0745c1c6ca..c3079fddd5 100644
--- a/src/modify.cpp
+++ b/src/modify.cpp
@@ -64,6 +64,7 @@ Modify::Modify(LAMMPS *lmp) : Pointers(lmp)
 
   end_of_step_every = NULL;
 
+  list_dofflag = NULL;
   list_timeflag = NULL;
 
   nfix_restart_global = 0;
@@ -135,6 +136,7 @@ Modify::~Modify()
   delete [] list_min_energy;
 
   delete [] end_of_step_every;
+  delete [] list_dofflag;
   delete [] list_timeflag;
 
   restart_deallocate();
@@ -184,6 +186,8 @@ void Modify::init()
   list_init(MIN_POST_FORCE,n_min_post_force,list_min_post_force);
   list_init(MIN_ENERGY,n_min_energy,list_min_energy);
 
+  list_init_dofflag(n_dofflag,list_dofflag);
+
   // init each fix
   // not sure if now needs to come before compute init
   // used to b/c temperature computes called fix->dof() in their init,
@@ -898,6 +902,19 @@ int Modify::check_package(const char *package_fix_name)
   return 1;
 }
 
+/* ----------------------------------------------------------------------
+   loop over fixes with dof() method
+   accumulate # of DOFs removed by fixes and return it
+   called by temperature computes
+------------------------------------------------------------------------- */
+
+int Modify::adjust_dof_fix(int igroup)
+{
+  int n = 0;
+  for (int ifix = 0; ifix < n_dofflag; ifix++)
+    n += fix[list_dofflag[ifix]]->dof(igroup);
+}
+
 /* ----------------------------------------------------------------------
    add a new compute
 ------------------------------------------------------------------------- */
@@ -1282,6 +1299,27 @@ void Modify::list_init_thermo_energy(int mask, int &n, int *&list)
     if (fmask[i] & mask && fix[i]->thermo_energy) list[n++] = i;
 }
 
+
+/* ----------------------------------------------------------------------
+   create list of fix indices for thermo energy fixes
+   only added to list if fix has THERMO_ENERGY mask
+   and its thermo_energy flag was set via fix_modify
+------------------------------------------------------------------------- */
+
+void Modify::list_init_dofflag(int &n, int *&list)
+{
+  delete [] list;
+
+  n = 0;
+  for (int i = 0; i < nfix; i++)
+    if (fix[i]->dof_flag) n++;
+  list = new int[n];
+
+  n = 0;
+  for (int i = 0; i < nfix; i++)
+    if (fix[i]->dof_flag) list[n++] = i;
+}
+
 /* ----------------------------------------------------------------------
    create list of compute indices for computes which store invocation times
 ------------------------------------------------------------------------- */
diff --git a/src/modify.h b/src/modify.h
index 4ba61d79ea..fb4a8f5057 100644
--- a/src/modify.h
+++ b/src/modify.h
@@ -88,6 +88,7 @@ class Modify : protected Pointers {
   void delete_fix(const char *);
   int find_fix(const char *);
   int check_package(const char *);
+  int adjust_dof_fix(int);
 
   void add_compute(int, char **, int trysuffix=0);
   void modify_compute(int, char **);
@@ -121,6 +122,9 @@ class Modify : protected Pointers {
 
   int *end_of_step_every;
 
+  int n_dofflag;             // list of fixes with dof() method
+  int *list_dofflag;
+
   int n_timeflag;            // list of computes that store time invocation
   int *list_timeflag;
 
@@ -137,6 +141,7 @@ class Modify : protected Pointers {
   void list_init(int, int &, int *&);
   void list_init_end_of_step(int, int &, int *&);
   void list_init_thermo_energy(int, int &, int *&);
+  void list_init_dofflag(int &, int *&);
   void list_init_compute();
 
  private:
diff --git a/src/molecule.cpp b/src/molecule.cpp
index f72bebe802..b7c9984892 100644
--- a/src/molecule.cpp
+++ b/src/molecule.cpp
@@ -536,10 +536,12 @@ void Molecule::charges(char *line)
 void Molecule::diameters(char *line)
 {
   int tmp;
+  maxradius = 0.0;
   for (int i = 0; i < natoms; i++) {
     readline(line);
     sscanf(line,"%d %lg",&tmp,&radius[i]);
     radius[i] *= 0.5;
+    maxradius = MAX(maxradius,radius[i]);
   }
 
   for (int i = 0; i < natoms; i++)
diff --git a/src/molecule.h b/src/molecule.h
index b8325d6d81..8ba3faaaec 100644
--- a/src/molecule.h
+++ b/src/molecule.h
@@ -90,7 +90,8 @@ class Molecule : protected Pointers {
   double ex[3],ey[3],ez[3]; // principal axes of molecule in space coords
   double quat[4];           // quaternion for orientation of molecule
 
-  double molradius;    // radius of molecule from COM,
+  double maxradius;    // max radius of any atom in molecule
+  double molradius;    // radius of molecule from geometric center
                        // including finite-size particle radii 
   int comatom;         // index (1-Natom) of atom closest to COM
   double maxextent;    // furthest any atom in molecule is from comatom