diff --git a/doc/src/Section_howto.txt b/doc/src/Section_howto.txt
index aedc2fc64f..e83f5399c0 100644
--- a/doc/src/Section_howto.txt
+++ b/doc/src/Section_howto.txt
@@ -1974,9 +1974,12 @@ The extract functions return a pointer to various global or per-atom
 quantities stored in LAMMPS or to values calculated by a compute, fix,
 or variable.  The pointer returned by the extract_global() function
 can be used as a permanent reference to a value which may change.  For
-the other extract functions, the underlying storage may be reallocated
-as LAMMPS runs, so you need to re-call the function to assure a
-current pointer or returned value(s).
+the extract_atom() method, see the extract() method in the
+src/atom.cpp file for a list of valid per-atom properties.  New names
+could easily be added if the property you want is not listed.  For the
+other extract functions, the underlying storage may be reallocated as
+LAMMPS runs, so you need to re-call the function to assure a current
+pointer or returned value(s).
 
 The lammps_reset_box() function resets the size and shape of the
 simulation box, e.g. as part of restoring a previously extracted and
@@ -1992,11 +1995,15 @@ keyword as a double precision value.
 The lammps_get_natoms() function returns the total number of atoms in
 the system and can be used by the caller to allocate space for the
 lammps_gather_atoms() and lammps_scatter_atoms() functions.  The
-gather function collects atom info of the requested type (atom coords,
-types, forces, etc) from all processors, orders them by atom ID, and
-returns a full list to each calling processor.  The scatter function
-does the inverse.  It distributes the same kinds of values, 
+gather function collects peratom info of the requested type (atom
+coords, types, forces, etc) from all processors, orders them by atom
+ID, and returns a full list to each calling processor.  The scatter
+function does the inverse.  It distributes the same peratom values,
 passed by the caller, to each atom owned by individual processors.
+Both methods are thus a means to extract or assign (overwrite) any
+peratom quantities within LAMMPS.  See the extract() method in the
+src/atom.cpp file for a list of valid per-atom properties.  New names
+could easily be added if the property you want is not listed.
 
 The lammps_create_atoms() function takes a list of N atoms as input
 with atom types and coords (required), an optionally atom IDs and
diff --git a/doc/src/Section_python.txt b/doc/src/Section_python.txt
index 50807e2d95..b994a56402 100644
--- a/doc/src/Section_python.txt
+++ b/doc/src/Section_python.txt
@@ -594,10 +594,10 @@ flag = lmp.set_variable(name,value)       # set existing named string-style vari
 value = lmp.get_thermo(name)              # return current value of a thermo keyword
 
 natoms = lmp.get_natoms()                 # total # of atoms as int
-data = lmp.gather_atoms(name,type,count)  # return atom attribute of all atoms gathered into data, ordered by atom ID
+data = lmp.gather_atoms(name,type,count)  # return per-atom property of all atoms gathered into data, ordered by atom ID
                                           # name = "x", "charge", "type", etc
                                           # count = # of per-atom values, 1 or 3, etc
-lmp.scatter_atoms(name,type,count,data)   # scatter atom attribute of all atoms from data, ordered by atom ID
+lmp.scatter_atoms(name,type,count,data)   # scatter per-atom property to all atoms from data, ordered by atom ID
                                           # name = "x", "charge", "type", etc
                                           # count = # of per-atom values, 1 or 3, etc :pre
 
@@ -656,10 +656,10 @@ argument.
 For extract_atom(), a pointer to internal LAMMPS atom-based data is
 returned, which you can use via normal Python subscripting.  See the
 extract() method in the src/atom.cpp file for a list of valid names.
-Again, new names could easily be added.  A pointer to a vector of
-doubles or integers, or a pointer to an array of doubles (double **)
-or integers (int **) is returned.  You need to specify the appropriate
-data type via the type argument.
+Again, new names could easily be added if the property you want is not
+listed.  A pointer to a vector of doubles or integers, or a pointer to
+an array of doubles (double **) or integers (int **) is returned.  You
+need to specify the appropriate data type via the type argument.
 
 For extract_compute() and extract_fix(), the global, per-atom, or
 local data calculated by the compute or fix can be accessed.  What is
@@ -689,12 +689,16 @@ specified group.
 The get_natoms() method returns the total number of atoms in the
 simulation, as an int.
 
-The gather_atoms() method returns a ctypes vector of ints or doubles
-as specified by type, of length count*natoms, for the property of all
-the atoms in the simulation specified by name, ordered by count and
-then by atom ID.  The vector can be used via normal Python
-subscripting.  If atom IDs are not consecutively ordered within
-LAMMPS, a None is returned as indication of an error.
+The gather_atoms() method allows any per-atom property (coordinates,
+velocities, etc) to be extracted from LAMMPS.  It returns a ctypes
+vector of ints or doubles as specified by type, of length
+count*natoms, for the named property for all atoms in the simulation.
+The data is ordered by count and then by atom ID.  See the extract()
+method in the src/atom.cpp file for a list of valid names.  Again, new
+names could easily be added if the property you want is missing.  The
+vector can be used via normal Python subscripting.  If atom IDs are
+not consecutively ordered within LAMMPS, a None is returned as
+indication of an error.
 
 Note that the data structure gather_atoms("x") returns is different
 from the data structure returned by extract_atom("x") in four ways.
@@ -711,14 +715,18 @@ assigning a new values to the extract_atom() array.  To do this with
 the gather_atoms() vector, you need to change values in the vector,
 then invoke the scatter_atoms() method.
 
-The scatter_atoms() method takes a vector of ints or doubles as
-specified by type, of length count*natoms, for the property of all the
-atoms in the simulation specified by name, ordered by bount and then
-by atom ID.  It uses the vector of data to overwrite the corresponding
-properties for each atom inside LAMMPS.  This requires LAMMPS to have
-its "map" option enabled; see the "atom_modify"_atom_modify.html
-command for details.  If it is not, or if atom IDs are not
-consecutively ordered, no coordinates are reset.
+The scatter_atoms() method allows any per-atom property (coordinates,
+velocities, etc) to be inserted into LAMMPS, overwriting the current
+property.  It takes a vector of ints or doubles as specified by type,
+of length count*natoms, for the named property for all atoms in the
+simulation.  The data should be ordered by count and then by atom ID.
+See the extract() method in the src/atom.cpp file for a list of valid
+names.  Again, new names could easily be added if the property you
+want is missing.  It uses the vector of data to overwrite the
+corresponding properties for each atom inside LAMMPS.  This requires
+LAMMPS to have its "map" option enabled; see the
+"atom_modify"_atom_modify.html command for details.  If it is not, or
+if atom IDs are not consecutively ordered, no coordinates are reset.
 
 The array of coordinates passed to scatter_atoms() must be a ctypes
 vector of ints or doubles, allocated and initialized something like
@@ -734,7 +742,7 @@ x\[2\] = z coord of atom with ID 1
 x\[3\] = x coord of atom with ID 2
 ...
 x\[n3-1\] = z coord of atom with ID natoms
-lmp.scatter_coords("x",1,3,x) :pre
+lmp.scatter_atoms("x",1,3,x) :pre
 
 Alternatively, you can just change values in the vector returned by
 gather_atoms("x",1,3), since it is a ctypes vector of doubles.
diff --git a/doc/src/compute_modify.txt b/doc/src/compute_modify.txt
index acf14526a7..637f9b5e41 100644
--- a/doc/src/compute_modify.txt
+++ b/doc/src/compute_modify.txt
@@ -14,27 +14,29 @@ compute_modify compute-ID keyword value ... :pre
 
 compute-ID = ID of the compute to modify :ulb,l
 one or more keyword/value pairs may be listed :l
-keyword = {extra} or {dynamic} :l
-  {extra} value = N
+keyword = {extra/dof} or {extra} or {dynamic/dof} or {dynamic} :l
+  {extra/dof} value = N
     N = # of extra degrees of freedom to subtract
-  {dynamic} value = {yes} or {no}
-    yes/no = do or do not recompute the number of atoms contributing to the temperature :pre
+  {extra} syntax is identical to {extra/dof}, will be disabled at some point
+  {dynamic/dof} value = {yes} or {no}
+    yes/no = do or do not recompute the number of atoms contributing to the temperature
+  {dynamic} syntax is identical to {dynamic/dof}, will be disabled at some point :pre
 :ule
 
 [Examples:]
 
-compute_modify myTemp extra 0
-compute_modify newtemp dynamic yes extra 600 :pre
+compute_modify myTemp extra/dof 0
+compute_modify newtemp dynamic/dof yes extra/dof 600 :pre
 
 [Description:]
 
 Modify one or more parameters of a previously defined compute.  Not
 all compute styles support all parameters.
 
-The {extra} keyword refers to how many degrees-of-freedom are
-subtracted (typically from 3N) as a normalizing factor in a
-temperature computation.  Only computes that compute a temperature use
-this option.  The default is 2 or 3 for "2d or 3d
+The {extra/dof} or {extra} keyword refers to how many
+degrees-of-freedom are subtracted (typically from 3N) as a normalizing
+factor in a temperature computation.  Only computes that compute a
+temperature use this option.  The default is 2 or 3 for "2d or 3d
 systems"_dimension.html which is a correction factor for an ensemble
 of velocities with zero total linear momentum. For compute
 temp/partial, if one or more velocity components are excluded, the
@@ -43,14 +45,18 @@ number for the {extra} parameter if you need to add
 degrees-of-freedom.  See the "compute
 temp/asphere"_compute_temp_asphere.html command for an example.
 
-The {dynamic} keyword determines whether the number of atoms N in the
-compute group is re-computed each time a temperature is computed.
-Only compute styles that calculate a temperature use this option.  By
-default, N is assumed to be constant.  If you are adding atoms to the
-system (see the "fix pour"_fix_pour.html or "fix
-deposit"_fix_deposit.html commands) or expect atoms to be lost
-(e.g. due to evaporation), then this option should be used to insure
-the temperature is correctly normalized.
+The {dynamic/dof} or {dynamic} keyword determines whether the number
+of atoms N in the compute group is re-computed each time a temperature
+is computed.  Only compute styles that calculate a temperature use
+this option.  By default, N is assumed to be constant.  If you are
+adding atoms to the system (see the "fix pour"_fix_pour.html, "fix
+deposit"_fix_deposit.html and "fix gcmc"_fix_gcmc.html commands) or
+expect atoms to be lost (e.g. due to evaporation), then this option
+should be used to insure the temperature is correctly normalized.
+
+NOTE: The {extra} and {dynamic} keywords should not be used as they
+are deprecated (March 2017) and will eventually be disabled.  Instead,
+use the equivalent {extra/dof} and {dynamic/dof} keywords.
 
 [Restrictions:] none
 
@@ -60,5 +66,5 @@ the temperature is correctly normalized.
 
 [Default:]
 
-The option defaults are extra = 2 or 3 for 2d or 3d systems and
-dynamic = no.
+The option defaults are extra/dof = 2 or 3 for 2d or 3d systems and
+dynamic/dof = no.
diff --git a/doc/src/compute_sna_atom.txt b/doc/src/compute_sna_atom.txt
index 0332ab9c1c..e2df706473 100644
--- a/doc/src/compute_sna_atom.txt
+++ b/doc/src/compute_sna_atom.txt
@@ -24,7 +24,7 @@ twojmax = band limit for bispectrum components (non-negative integer) :l
 R_1, R_2,... = list of cutoff radii, one for each type (distance units) :l
 w_1, w_2,... = list of neighbor weights, one for each type  :l
 zero or more keyword/value pairs may be appended :l
-keyword = {diagonal} or {rmin0} or {switchflag} :l
+keyword = {diagonal} or {rmin0} or {switchflag} or {bzeroflag} :l
   {diagonal} value = {0} or {1} or {2} or {3}
      {0} = all j1, j2, j <= twojmax, j2 <= j1
      {1} = subset satisfying j1 == j2
@@ -33,7 +33,10 @@ keyword = {diagonal} or {rmin0} or {switchflag} :l
   {rmin0} value = parameter in distance to angle conversion (distance units)
   {switchflag} value = {0} or {1}
      {0} = do not use switching function
-     {1} = use switching function :pre
+     {1} = use switching function
+  {bzeroflag} value = {0} or {1}
+     {0} = do not subtract B0
+     {1} = subtract B0 :pre
 :ule
 
 [Examples:]
@@ -153,6 +156,12 @@ ordered in which they are listed
 The keyword {switchflag} can be used to turn off the switching
 function.
 
+The keyword {bzeroflag} determines whether or not {B0}, the bispectrum
+components of an atom with no neighbors, are subtracted from
+the calculated bispectrum components. This optional keyword is only
+available for compute {sna/atom}, as {snad/atom} and {snav/atom}
+are unaffected by the removal of constant terms.
+
 NOTE: If you have a bonded system, then the settings of
 "special_bonds"_special_bonds.html command can remove pairwise
 interactions between atoms in the same bond, angle, or dihedral.  This
@@ -222,7 +231,7 @@ LAMMPS"_Section_start.html#start_3 section for more info.
 [Default:]
 
 The optional keyword defaults are {diagonal} = 0, {rmin0} = 0,
-{switchflag} = 1.
+{switchflag} = 1, {bzeroflag} = 0.
 
 :line
 
diff --git a/doc/src/fix_deposit.txt b/doc/src/fix_deposit.txt
index a1dd5f6434..477c14ea89 100644
--- a/doc/src/fix_deposit.txt
+++ b/doc/src/fix_deposit.txt
@@ -150,6 +150,12 @@ treated as rigid bodies, use the {rigid} keyword, specifying as its
 value the ID of a separate "fix rigid/small"_fix_rigid.html
 command which also appears in your input script.
 
+NOTE: If you wish the new rigid molecules (and other rigid molecules)
+to be thermostatted correctly via "fix rigid/small/nvt"_fix_rigid.html
+or "fix rigid/small/npt"_fix_rigid.html, then you need to use the
+"fix_modify dynamic/dof yes" command for the rigid fix.  This is to
+inform that fix that the molecule count will vary dynamically.
+
 If you wish to insert molecules via the {mol} keyword, that will have
 their bonds or angles constrained via SHAKE, use the {shake} keyword,
 specifying as its value the ID of a separate "fix
diff --git a/doc/src/fix_gcmc.txt b/doc/src/fix_gcmc.txt
index 4eab458a65..1a419628e9 100644
--- a/doc/src/fix_gcmc.txt
+++ b/doc/src/fix_gcmc.txt
@@ -26,6 +26,8 @@ zero or more keyword/value pairs may be appended to args :l
 keyword = {mol}, {region}, {maxangle}, {pressure}, {fugacity_coeff}, {full_energy}, {charge}, {group}, {grouptype}, {intra_energy}, {tfac_insert}, or {overlap_cutoff}
   {mol} value = template-ID
     template-ID = ID of molecule template specified in a separate "molecule"_molecule.html command
+  {rigid} value = fix-ID
+    fix-ID = ID of "fix rigid/small"_fix_rigid.html command
   {shake} value = fix-ID
     fix-ID = ID of "fix shake"_fix_shake.html command
   {region} value = region-ID
@@ -161,6 +163,17 @@ soon generate an error when it tries to find bonded neighbors.  LAMMPS will
 warn you if any of the atoms eligible for deletion have a non-zero
 molecule ID, but does not check for this at the time of deletion.
 
+If you wish to insert molecules via the {mol} keyword, that will be
+treated as rigid bodies, use the {rigid} keyword, specifying as its
+value the ID of a separate "fix rigid/small"_fix_rigid.html
+command which also appears in your input script.
+
+NOTE: If you wish the new rigid molecules (and other rigid molecules)
+to be thermostatted correctly via "fix rigid/small/nvt"_fix_rigid.html
+or "fix rigid/small/npt"_fix_rigid.html, then you need to use the
+"fix_modify dynamic/dof yes" command for the rigid fix.  This is to
+inform that fix that the molecule count will vary dynamically.
+
 If you wish to insert molecules via the {mol} keyword, that will have
 their bonds or angles constrained via SHAKE, use the {shake} keyword,
 specifying as its value the ID of a separate "fix
@@ -343,10 +356,6 @@ may no longer exist since it might have been deleted by the first
 fix gcmc command. An existing template molecule will need to be
 referenced by the user for each subsequent fix gcmc command.
 
-Because molecule insertion does not work in combination with
-fix rigid, simulataneous use of fix rigid or fix rigid/small
-with this fix is not allowed.
-
 [Related commands:]
 
 "fix atom/swap"_fix_atom_swap.html,
diff --git a/doc/src/fix_modify.txt b/doc/src/fix_modify.txt
index 9c95cdc452..5e097bb34f 100644
--- a/doc/src/fix_modify.txt
+++ b/doc/src/fix_modify.txt
@@ -14,11 +14,13 @@ fix_modify fix-ID keyword value ... :pre
 
 fix-ID = ID of the fix to modify :ulb,l
 one or more keyword/value pairs may be appended :l
-keyword = {temp} or {press} or {energy} or {respa} :l
+keyword = {temp} or {press} or {energy} or {respa} or {dynamic/dof} :l
   {temp} value = compute ID that calculates a temperature
   {press} value = compute ID that calculates a pressure
   {energy} value = {yes} or {no}
-  {respa} value = {1} to {max respa level} or {0} (for outermost level) :pre
+  {respa} value = {1} to {max respa level} or {0} (for outermost level)
+  {dynamic/dof} value = {yes} or {no}
+    yes/no = do or do not recompute the number of degrees of freedom (DOF) contributing to the temperature :pre
 :ule
 
 [Examples:]
@@ -78,6 +80,27 @@ enabled to support this feature; if not, {fix_modify} will report an
 error. Active fixes with a custom RESPA level setting are reported
 with their specified level at the beginning of a r-RESPA run.
 
+The {dynamic/dof} keyword determines whether the number of atoms N in
+the fix group and their associated degrees of freedom are re-computed
+each time a temperature is computed.  Only fix styles that calculate
+their own internal temperature use this option.  Currently this is
+only the "fix rigid/nvt/small"_fix_rigid.html and "fix
+rigid/npt/small"_fix_rigid.html commands for the purpose of
+thermostatting rigid body translation and rotation.  By default, N and
+their DOF are assumed to be constant.  If you are adding atoms or
+molecules to the system (see the "fix pour"_fix_pour.html, "fix
+deposit"_fix_deposit.html, and "fix gcmc"_fix_gcmc.html commands) or
+expect atoms or molecules to be lost (e.g. due to exiting the
+simulation box or via "fix evaporation"_fix_evaporation.html), then
+this option should be used to insure the temperature is correctly
+normalized.
+
+NOTE: Other thermostatting fixes, such as "fix nvt"_fix_nh.html, do
+not use the {dynamic/dof} keyword because they use a temperature
+compute to calculate temperature.  See the "compute_modify
+dynamic/dof"_compute_modify.html command for a similar way to insure
+correct temperature normalization for those thermostats.
+
 [Restrictions:] none
 
 [Related commands:]
diff --git a/doc/src/fix_pour.txt b/doc/src/fix_pour.txt
index f2c10184dd..54f78287e0 100644
--- a/doc/src/fix_pour.txt
+++ b/doc/src/fix_pour.txt
@@ -117,6 +117,12 @@ treated as rigid bodies, use the {rigid} keyword, specifying as its
 value the ID of a separate "fix rigid/small"_fix_rigid.html
 command which also appears in your input script.
 
+NOTE: If you wish the new rigid molecules (and other rigid molecules)
+to be thermostatted correctly via "fix rigid/small/nvt"_fix_rigid.html
+or "fix rigid/small/npt"_fix_rigid.html, then you need to use the
+"fix_modify dynamic/dof yes" command for the rigid fix.  This is to
+inform that fix that the molecule count will vary dynamically.
+
 If you wish to insert molecules via the {mol} keyword, that will have
 their bonds or angles constrained via SHAKE, use the {shake} keyword,
 specifying as its value the ID of a separate "fix
diff --git a/doc/src/pair_snap.txt b/doc/src/pair_snap.txt
index a24c6c3160..ab7313832c 100644
--- a/doc/src/pair_snap.txt
+++ b/doc/src/pair_snap.txt
@@ -131,14 +131,15 @@ The SNAP parameter file can contain blank and comment lines (start
 with #) anywhere. Each non-blank non-comment line must contain one
 keyword/value pair. The required keywords are {rcutfac} and
 {twojmax}. Optional keywords are {rfac0}, {rmin0}, {diagonalstyle},
-and {switchflag}.
+{switchflag}, and {bzeroflag}.
 
 The default values for these keywords are
 
 {rfac0} = 0.99363
 {rmin0} = 0.0
 {diagonalstyle} = 3
-{switchflag} = 0 :ul
+{switchflag} = 0
+{bzeroflag} = 1 :ul
 
 Detailed definitions of these keywords are given on the "compute
 sna/atom"_compute_sna_atom.html doc page.
diff --git a/examples/accelerate/in.lc b/examples/accelerate/in.lc
index 66e3916fa4..4dbdbb8554 100644
--- a/examples/accelerate/in.lc
+++ b/examples/accelerate/in.lc
@@ -30,7 +30,7 @@ set	      group all quat/random 982381
 compute	       rot all temp/asphere
 group	       spheroid type 1
 variable       dof equal count(spheroid)+3
-compute_modify rot extra ${dof}
+compute_modify rot extra/dof ${dof}
 
 velocity      all create 2.4 41787 loop geom
 
@@ -45,7 +45,7 @@ thermo	      100
 # equilibration run
 
 fix	       1 all npt/asphere temp 2.4 2.4 0.1 iso 5.0 8.0 0.1
-compute_modify 1_temp extra ${dof}
+compute_modify 1_temp extra/dof ${dof}
 run	       200
 
 # dynamics run
diff --git a/examples/deposit/in.deposit.atom b/examples/deposit/in.deposit.atom
index 8f29040399..3e48276bb2 100644
--- a/examples/deposit/in.deposit.atom
+++ b/examples/deposit/in.deposit.atom
@@ -23,7 +23,7 @@ region          mobile block 0 5 0 5 2 INF
 group		mobile region mobile
 
 compute		add addatoms temp
-compute_modify	add dynamic yes extra 0
+compute_modify	add dynamic/dof yes extra/dof 0
 
 fix		1 addatoms nve
 fix		2 mobile langevin 1.0 1.0 0.1 587283
diff --git a/examples/deposit/in.deposit.molecule b/examples/deposit/in.deposit.molecule
index 870d740724..16b0b9a72a 100644
--- a/examples/deposit/in.deposit.molecule
+++ b/examples/deposit/in.deposit.molecule
@@ -26,7 +26,7 @@ region          mobile block 0 5 0 5 2 INF
 group		mobile region mobile
 
 compute		add addatoms temp
-compute_modify	add dynamic yes extra 0
+compute_modify	add dynamic/dof yes extra/dof 0
 
 fix		1 addatoms nve
 fix		2 mobile langevin 0.1 0.1 0.1 587283
diff --git a/examples/deposit/in.deposit.molecule.shake b/examples/deposit/in.deposit.molecule.shake
index a1c141088b..7bd701c924 100644
--- a/examples/deposit/in.deposit.molecule.shake
+++ b/examples/deposit/in.deposit.molecule.shake
@@ -26,7 +26,7 @@ region          mobile block 0 5 0 5 2 INF
 group		mobile region mobile
 
 compute		add addatoms temp
-compute_modify	add dynamic yes extra 0
+compute_modify	add dynamic/dof yes extra/dof 0
 
 fix		1 addatoms nve
 fix		2 mobile langevin 0.1 0.1 0.1 587283
diff --git a/examples/ellipse/in.ellipse.gayberne b/examples/ellipse/in.ellipse.gayberne
index 19e4e2414a..fe783ac6de 100644
--- a/examples/ellipse/in.ellipse.gayberne
+++ b/examples/ellipse/in.ellipse.gayberne
@@ -19,7 +19,7 @@ set	     group all quat/random 18238
 compute	     rot all temp/asphere
 group	     spheroid type 1
 variable     dof equal count(spheroid)+2
-compute_modify rot extra ${dof}
+compute_modify rot extra/dof ${dof}
 
 velocity     all create 2.4 87287 loop geom
 
@@ -52,7 +52,7 @@ fix	     1 all npt/asphere temp 2.0 2.0 0.1 iso 0.0 1.0 1.0 &
 	       mtk no pchain 0 tchain 1
 fix	     2 all enforce2d
 
-compute_modify 1_temp extra ${dof}
+compute_modify 1_temp extra/dof ${dof}
 
 # equilibrate to shrink box around dilute system
 
diff --git a/examples/ellipse/in.ellipse.resquared b/examples/ellipse/in.ellipse.resquared
index df52aef66c..82398987f8 100644
--- a/examples/ellipse/in.ellipse.resquared
+++ b/examples/ellipse/in.ellipse.resquared
@@ -19,7 +19,7 @@ set	     group all quat/random 18238
 compute	     rot all temp/asphere
 group	     spheroid type 1
 variable     dof equal count(spheroid)+2
-compute_modify rot extra ${dof}
+compute_modify rot extra/dof ${dof}
 
 velocity     all create 2.4 87287 loop geom
 
@@ -52,7 +52,7 @@ fix	     1 all npt/asphere temp 2.0 2.0 0.1 iso 0.0 1.0 1.0 &
 	       mtk no pchain 0 tchain 1
 fix	     2 all enforce2d
 
-compute_modify 1_temp extra ${dof}
+compute_modify 1_temp extra/dof ${dof}
 
 # equilibrate to shrink box around dilute system
 
diff --git a/examples/granregion/in.granregion.box b/examples/granregion/in.granregion.box
index 91f06744d3..f1f20ad79c 100644
--- a/examples/granregion/in.granregion.box
+++ b/examples/granregion/in.granregion.box
@@ -29,15 +29,15 @@ fix             ins all pour 100 2 4767548 vol 0.4 10 &
 timestep	0.005
 
 compute         1 all temp
-compute_modify  1 dynamic yes
+compute_modify  1 dynamic/dof yes
 
 compute         2 all temp/sphere
-compute_modify  2 dynamic yes
+compute_modify  2 dynamic/dof yes
 
 thermo		100
 thermo_style	custom step atoms temp c_1 c_2 press
 thermo_modify   lost ignore
-compute_modify  thermo_temp dynamic yes
+compute_modify  thermo_temp dynamic/dof yes
 
 #dump		2 all image 100 image.*.jpg type type &
 #		zoom 1.4 adiam 1.0 box no 0.0 axes yes 0.9 0.03
diff --git a/examples/pour/in.pour b/examples/pour/in.pour
index 9b24b2906c..332ccf8b62 100644
--- a/examples/pour/in.pour
+++ b/examples/pour/in.pour
@@ -33,7 +33,7 @@ compute		1 all erotate/sphere
 thermo_style	custom step atoms ke c_1 vol
 thermo		1000
 thermo_modify	lost ignore norm no
-compute_modify	thermo_temp dynamic yes
+compute_modify	thermo_temp dynamic/dof yes
 
 #dump		id all atom 1000 dump.pour
 
diff --git a/examples/pour/in.pour.2d b/examples/pour/in.pour.2d
index 02fc794ffd..5fa057c746 100644
--- a/examples/pour/in.pour.2d
+++ b/examples/pour/in.pour.2d
@@ -39,7 +39,7 @@ compute		1 all erotate/sphere
 thermo_style	custom step atoms ke c_1 vol
 thermo		1000
 thermo_modify	lost ignore norm no
-compute_modify	thermo_temp dynamic yes
+compute_modify	thermo_temp dynamic/dof yes
 
 #dump		id all atom 250 dump.pour
 
diff --git a/examples/pour/in.pour.2d.molecule b/examples/pour/in.pour.2d.molecule
index 84869eb233..a65c72875f 100644
--- a/examples/pour/in.pour.2d.molecule
+++ b/examples/pour/in.pour.2d.molecule
@@ -46,7 +46,7 @@ compute		1 all erotate/sphere
 compute         Tsphere all temp/sphere
 thermo_style	custom step atoms ke c_1 vol
 thermo_modify	lost ignore norm no temp Tsphere
-compute_modify	Tsphere dynamic yes
+compute_modify	Tsphere dynamic/dof yes
 
 thermo		1000
 
diff --git a/examples/snap/Ta06A.snapparam b/examples/snap/Ta06A.snapparam
index 0627253341..7b30312f59 100644
--- a/examples/snap/Ta06A.snapparam
+++ b/examples/snap/Ta06A.snapparam
@@ -12,3 +12,4 @@ gamma 1
 rfac0 0.99363
 rmin0 0
 diagonalstyle 3
+bzeroflag 0
diff --git a/examples/snap/W_2940_2017_2.snapparam b/examples/snap/W_2940_2017_2.snapparam
index e0b20005e6..f17961bdd7 100644
--- a/examples/snap/W_2940_2017_2.snapparam
+++ b/examples/snap/W_2940_2017_2.snapparam
@@ -10,3 +10,4 @@ gamma 1
 rfac0 0.99363
 rmin0 0
 diagonalstyle 3
+bzeroflag 0
diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp
index e638cd8873..780df25bdf 100644
--- a/src/MC/fix_gcmc.cpp
+++ b/src/MC/fix_gcmc.cpp
@@ -60,7 +60,7 @@ using namespace MathConst;
 // this must be lower than MAXENERGYSIGNAL
 // by a large amount, so that it is still
 // less than total energy when negative
-// energy changes are adddd to MAXENERGYSIGNAL
+// energy changes are added to MAXENERGYSIGNAL
 
 #define MAXENERGYTEST 1.0e50
 
@@ -72,7 +72,7 @@ FixGCMC::FixGCMC(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg),
   idregion(NULL), full_flag(0), ngroups(0), groupstrings(NULL), ngrouptypes(0), grouptypestrings(NULL),
   grouptypebits(NULL), grouptypes(NULL), local_gas_list(NULL), atom_coord(NULL), random_equal(NULL), random_unequal(NULL), 
-  coords(NULL), imageflags(NULL), idshake(NULL)
+  coords(NULL), imageflags(NULL), idrigid(NULL), idshake(NULL), fixrigid(NULL), fixshake(NULL)
 {
   if (narg < 11) error->all(FLERR,"Illegal fix gcmc command");
 
@@ -182,8 +182,12 @@ FixGCMC::FixGCMC(LAMMPS *lmp, int narg, char **arg) :
   if (charge_flag && atom->q == NULL)
     error->all(FLERR,"Fix gcmc atom has charge, but atom style does not");
 
+  if (rigidflag && mode == ATOM)
+    error->all(FLERR,"Cannot use fix gcmc rigid and not molecule");
   if (shakeflag && mode == ATOM)
     error->all(FLERR,"Cannot use fix gcmc shake and not molecule");
+  if (rigidflag && shakeflag)
+    error->all(FLERR,"Cannot use fix gcmc rigid and shake");
 
   // setup of coords and imageflags array
 
@@ -241,11 +245,11 @@ void FixGCMC::options(int narg, char **arg)
   pressure_flag = false;
   pressure = 0.0;
   fugacity_coeff = 1.0;
+  rigidflag = 0;
   shakeflag = 0;
   charge = 0.0;
   charge_flag = false;
   full_flag = false;
-  idshake = NULL;
   ngroups = 0;
   int ngroupsmax = 0;
   groupstrings = NULL;
@@ -302,6 +306,14 @@ void FixGCMC::options(int narg, char **arg)
       charge = force->numeric(FLERR,arg[iarg+1]);
       charge_flag = true;
       iarg += 2;
+    } else if (strcmp(arg[iarg],"rigid") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix gcmc command");
+      int n = strlen(arg[iarg+1]) + 1;
+      delete [] idrigid;
+      idrigid = new char[n];
+      strcpy(idrigid,arg[iarg+1]);
+      rigidflag = 1;
+      iarg += 2;
     } else if (strcmp(arg[iarg],"shake") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix gcmc command");
       int n = strlen(arg[iarg+1]) + 1;
@@ -374,6 +386,7 @@ FixGCMC::~FixGCMC()
   memory->destroy(coords);
   memory->destroy(imageflags);
 
+  delete [] idrigid;
   delete [] idshake;
 
   if (ngroups > 0) {
@@ -477,6 +490,21 @@ void FixGCMC::init()
                "Fix gcmc molecule command requires that "
                "atoms have molecule attributes");
 
+  // if rigidflag defined, check for rigid/small fix
+  // its molecule template must be same as this one
+
+  fixrigid = NULL;
+  if (rigidflag) {
+    int ifix = modify->find_fix(idrigid);
+    if (ifix < 0) error->all(FLERR,"Fix gcmc rigid fix does not exist");
+    fixrigid = modify->fix[ifix];
+    int tmp;
+    if (onemols != (Molecule **) fixrigid->extract("onemol",tmp))
+      error->all(FLERR,
+                 "Fix gcmc and fix rigid/small not using "
+                 "same molecule template ID");
+  }
+
   // if shakeflag defined, check for SHAKE fix
   // its molecule template must be same as this one
 
@@ -491,13 +519,6 @@ void FixGCMC::init()
                  "same molecule template ID");
   }
 
-  // check for fix rigid
-
-  for (int irigid = 0; irigid < modify->nfix; irigid++) {
-    if (strncmp(modify->fix[irigid]->style,"rigid",5) == 0)
-      error->all(FLERR,"Fix gcmc can not currently be used with any rigid fix");
-  }
-
   if (domain->dimension == 2)
     error->all(FLERR,"Cannot use fix gcmc in a 2d simulation");
 
@@ -1353,10 +1374,15 @@ void FixGCMC::attempt_molecule_insertion()
         modify->create_attribute(m);
       }
     }
-      
-    if (shakeflag)
+
+    // FixRigidSmall::set_molecule stores rigid body attributes
+    // FixShake::set_molecule stores shake info for molecule
+    
+    if (rigidflag)
+      fixrigid->set_molecule(nlocalprev,maxtag_all,imol,com_coord,vnew,quat);
+    else if (shakeflag)
       fixshake->set_molecule(nlocalprev,maxtag_all,imol,com_coord,vnew,quat);
-      
+
     atom->natoms += natoms_per_molecule;
     if (atom->natoms < 0)
       error->all(FLERR,"Too many total atoms");
@@ -2011,7 +2037,12 @@ void FixGCMC::attempt_molecule_insertion_full()
     }
   }
 
-  if (shakeflag)
+  // FixRigidSmall::set_molecule stores rigid body attributes
+  // FixShake::set_molecule stores shake info for molecule
+
+  if (rigidflag)
+    fixrigid->set_molecule(nlocalprev,maxtag_all,imol,com_coord,vnew,quat);
+  else if (shakeflag)
     fixshake->set_molecule(nlocalprev,maxtag_all,imol,com_coord,vnew,quat);
 
   atom->natoms += natoms_per_molecule;
@@ -2095,7 +2126,7 @@ double FixGCMC::energy(int i, int itype, tagint imolecule, double *coord)
     int jtype = type[j];
 
     // if overlap check requested, if overlap,
-    // return signal value = MAXENERGYSIGNAL 
+    // return signal value for energy 
 
     if (overlap_flag && rsq < overlap_cutoff)
       return MAXENERGYSIGNAL;
@@ -2145,7 +2176,7 @@ double FixGCMC::energy_full()
   int vflag = 0;
 
   // if overlap check requested, if overlap,
-  // return signal value = MAXENERGYSIGNAL 
+  // return signal value for energy 
 
   if (overlap_flag) {
     double delx,dely,delz,rsq;
diff --git a/src/MC/fix_gcmc.h b/src/MC/fix_gcmc.h
index 8261f6e38c..9b2184dda2 100644
--- a/src/MC/fix_gcmc.h
+++ b/src/MC/fix_gcmc.h
@@ -128,9 +128,9 @@ class FixGCMC : public Fix {
   int imol,nmol;
   double **coords;
   imageint *imageflags;
-  class Fix *fixshake;
-  int shakeflag;
-  char *idshake;
+  class Fix *fixrigid, *fixshake;
+  int rigidflag, shakeflag;
+  char *idrigid, *idshake;
   int triclinic;                         // 0 = orthog box, 1 = triclinic
 
   class Compute *c_pe;
diff --git a/src/MISC/fix_deposit.cpp b/src/MISC/fix_deposit.cpp
index d841482f8e..9c1082f816 100644
--- a/src/MISC/fix_deposit.cpp
+++ b/src/MISC/fix_deposit.cpp
@@ -234,7 +234,7 @@ void FixDeposit::init()
   fixrigid = NULL;
   if (rigidflag) {
     int ifix = modify->find_fix(idrigid);
-    if (ifix < 0) error->all(FLERR,"Fix pour rigid fix does not exist");
+    if (ifix < 0) error->all(FLERR,"Fix deposit rigid fix does not exist");
     fixrigid = modify->fix[ifix];
     int tmp;
     if (onemols != (Molecule **) fixrigid->extract("onemol",tmp))
diff --git a/src/RIGID/fix_rigid_nh_small.cpp b/src/RIGID/fix_rigid_nh_small.cpp
index b08c80b626..199c04bd9c 100644
--- a/src/RIGID/fix_rigid_nh_small.cpp
+++ b/src/RIGID/fix_rigid_nh_small.cpp
@@ -321,36 +321,8 @@ void FixRigidNHSmall::init()
 void FixRigidNHSmall::setup(int vflag)
 {
   FixRigidSmall::setup(vflag);
-
-  // total translational and rotational degrees of freedom
-
-  nf_t = dimension * nlocal_body;
-  if (dimension == 3) {
-    nf_r = dimension * nlocal_body;
-    for (int ibody = 0; ibody < nlocal_body; ibody++) {
-      Body *b = &body[ibody];
-      for (int k = 0; k < dimension; k++)
-        if (fabs(b->inertia[k]) < EPSILON) nf_r--;
-    }
-  } else if (dimension == 2) {
-    nf_r = nlocal_body;
-    for (int ibody = 0; ibody < nlocal_body; ibody++) {
-      Body *b = &body[ibody];
-      if (fabs(b->inertia[2]) < EPSILON) nf_r--;
-    }
-  }
-
-  double nf[2], nfall[2];
-  nf[0] = nf_t;
-  nf[1] = nf_r;
-  MPI_Allreduce(nf,nfall,2,MPI_DOUBLE,MPI_SUM,world);
-  nf_t = nfall[0];
-  nf_r = nfall[1];
-
-  g_f = nf_t + nf_r;
-  onednft = 1.0 + (double)(dimension) / (double)g_f;
-  onednfr = (double) (dimension) / (double)g_f;
-
+  compute_dof();
+  
   double mbody[3];
   akin_t = akin_r = 0.0;
   for (int ibody = 0; ibody < nlocal_body; ibody++) {
@@ -608,6 +580,7 @@ void FixRigidNHSmall::initial_integrate(int vflag)
 
   if (tstat_flag) {
     compute_temp_target();
+    if (dynamic) compute_dof();
     nhc_temp_integrate();
   }
 
@@ -1273,6 +1246,40 @@ void FixRigidNHSmall::nh_epsilon_dot()
   mtk_term2 /= g_f;
 }
 
+/* ---------------------------------------------------------------------- */
+
+void FixRigidNHSmall::compute_dof()
+{
+  // total translational and rotational degrees of freedom
+
+  nf_t = dimension * nlocal_body;
+  if (dimension == 3) {
+    nf_r = dimension * nlocal_body;
+    for (int ibody = 0; ibody < nlocal_body; ibody++) {
+      Body *b = &body[ibody];
+      for (int k = 0; k < dimension; k++)
+        if (fabs(b->inertia[k]) < EPSILON) nf_r--;
+    }
+  } else if (dimension == 2) {
+    nf_r = nlocal_body;
+    for (int ibody = 0; ibody < nlocal_body; ibody++) {
+      Body *b = &body[ibody];
+      if (fabs(b->inertia[2]) < EPSILON) nf_r--;
+    }
+  }
+
+  double nf[2], nfall[2];
+  nf[0] = nf_t;
+  nf[1] = nf_r;
+  MPI_Allreduce(nf,nfall,2,MPI_DOUBLE,MPI_SUM,world);
+  nf_t = nfall[0];
+  nf_r = nfall[1];
+
+  g_f = nf_t + nf_r;
+  onednft = 1.0 + (double)(dimension) / (double)g_f;
+  onednfr = (double) (dimension) / (double)g_f;
+}
+
 /* ----------------------------------------------------------------------
    pack entire state of Fix into one write
 ------------------------------------------------------------------------- */
@@ -1508,4 +1515,3 @@ void FixRigidNHSmall::deallocate_order()
   delete [] wdti2;
   delete [] wdti4;
 }
-
diff --git a/src/RIGID/fix_rigid_nh_small.h b/src/RIGID/fix_rigid_nh_small.h
index 07510a59e3..6b51933277 100644
--- a/src/RIGID/fix_rigid_nh_small.h
+++ b/src/RIGID/fix_rigid_nh_small.h
@@ -78,7 +78,8 @@ class FixRigidNHSmall : public FixRigidSmall {
   virtual void compute_temp_target();
   void compute_press_target();
   void nh_epsilon_dot();
-
+  void compute_dof();
+  
   void allocate_chain();
   void allocate_order();
   void deallocate_chain();
diff --git a/src/SNAP/compute_sna_atom.cpp b/src/SNAP/compute_sna_atom.cpp
index 326d2d620d..ad934535ab 100644
--- a/src/SNAP/compute_sna_atom.cpp
+++ b/src/SNAP/compute_sna_atom.cpp
@@ -34,7 +34,7 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) :
   radelem(NULL), wjelem(NULL)
 {
   double rmin0, rfac0;
-  int twojmax, switchflag;
+  int twojmax, switchflag, bzeroflag;
   radelem = NULL;
   wjelem = NULL;
 
@@ -48,6 +48,7 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) :
   diagonalstyle = 0;
   rmin0 = 0.0;
   switchflag = 1;
+  bzeroflag = 0;
 
   // offset by 1 to match up with types
 
@@ -100,19 +101,24 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) :
 	error->all(FLERR,"Illegal compute sna/atom command");
       switchflag = atoi(arg[iarg+1]);
       iarg += 2;
+    } else if (strcmp(arg[iarg],"bzeroflag") == 0) {
+      if (iarg+2 > narg)
+	error->all(FLERR,"Illegal compute sna/atom command");
+      bzeroflag = atoi(arg[iarg+1]);
+      iarg += 2;
     } else error->all(FLERR,"Illegal compute sna/atom command");
   }
 
   snaptr = new SNA*[comm->nthreads];
 #if defined(_OPENMP)
-#pragma omp parallel default(none) shared(lmp,rfac0,twojmax,rmin0,switchflag)
+#pragma omp parallel default(none) shared(lmp,rfac0,twojmax,rmin0,switchflag,bzeroflag)
 #endif
   {
     int tid = omp_get_thread_num();
 
     // always unset use_shared_arrays since it does not work with computes
     snaptr[tid] = new SNA(lmp,rfac0,twojmax,diagonalstyle,
-                          0 /*use_shared_arrays*/, rmin0,switchflag);
+                          0 /*use_shared_arrays*/, rmin0,switchflag,bzeroflag);
   }
 
   ncoeff = snaptr[0]->ncoeff;
diff --git a/src/SNAP/compute_snad_atom.cpp b/src/SNAP/compute_snad_atom.cpp
index efd4cafbca..73452427bd 100644
--- a/src/SNAP/compute_snad_atom.cpp
+++ b/src/SNAP/compute_snad_atom.cpp
@@ -34,7 +34,7 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) :
   radelem(NULL), wjelem(NULL)
 {
   double rfac0, rmin0;
-  int twojmax, switchflag;
+  int twojmax, switchflag, bzeroflag;
   radelem = NULL;
   wjelem = NULL;
 
@@ -48,7 +48,8 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) :
   diagonalstyle = 0;
   rmin0 = 0.0;
   switchflag = 1;
-
+  bzeroflag = 0;
+  
   // process required arguments
   memory->create(radelem,ntypes+1,"sna/atom:radelem"); // offset by 1 to match up with types
   memory->create(wjelem,ntypes+1,"sna/atom:wjelem");
@@ -98,14 +99,14 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) :
 
   snaptr = new SNA*[comm->nthreads];
 #if defined(_OPENMP)
-#pragma omp parallel default(none) shared(lmp,rfac0,twojmax,rmin0,switchflag)
+#pragma omp parallel default(none) shared(lmp,rfac0,twojmax,rmin0,switchflag,bzeroflag)
 #endif
   {
     int tid = omp_get_thread_num();
 
     // always unset use_shared_arrays since it does not work with computes
     snaptr[tid] = new SNA(lmp,rfac0,twojmax,diagonalstyle,
-                          0 /*use_shared_arrays*/, rmin0,switchflag);
+                          0 /*use_shared_arrays*/, rmin0,switchflag,bzeroflag);
   }
 
   ncoeff = snaptr[0]->ncoeff;
diff --git a/src/SNAP/compute_snav_atom.cpp b/src/SNAP/compute_snav_atom.cpp
index c1398864e0..f75b02fba7 100644
--- a/src/SNAP/compute_snav_atom.cpp
+++ b/src/SNAP/compute_snav_atom.cpp
@@ -34,7 +34,7 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) :
   radelem(NULL), wjelem(NULL)
 {
   double rfac0, rmin0;
-  int twojmax, switchflag;
+  int twojmax, switchflag, bzeroflag;
   radelem = NULL;
   wjelem = NULL;
 
@@ -50,6 +50,7 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) :
   diagonalstyle = 0;
   rmin0 = 0.0;
   switchflag = 1;
+  bzeroflag = 0;
 
   // process required arguments
   memory->create(radelem,ntypes+1,"sna/atom:radelem"); // offset by 1 to match up with types
@@ -100,14 +101,14 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) :
 
   snaptr = new SNA*[comm->nthreads];
 #if defined(_OPENMP)
-#pragma omp parallel default(none) shared(lmp,rfac0,twojmax,rmin0,switchflag)
+#pragma omp parallel default(none) shared(lmp,rfac0,twojmax,rmin0,switchflag,bzeroflag)
 #endif
   {
     int tid = omp_get_thread_num();
 
     // always unset use_shared_arrays since it does not work with computes
     snaptr[tid] = new SNA(lmp,rfac0,twojmax,diagonalstyle,
-                          0 /*use_shared_arrays*/, rmin0,switchflag);
+                          0 /*use_shared_arrays*/, rmin0,switchflag,bzeroflag);
   }
 
   ncoeff = snaptr[0]->ncoeff;
diff --git a/src/SNAP/pair_snap.cpp b/src/SNAP/pair_snap.cpp
index dc84b0be05..06c2e48488 100644
--- a/src/SNAP/pair_snap.cpp
+++ b/src/SNAP/pair_snap.cpp
@@ -1411,7 +1411,7 @@ void PairSNAP::coeff(int narg, char **arg)
     int tid = omp_get_thread_num();
     sna[tid] = new SNA(lmp,rfac0,twojmax,
                        diagonalstyle,use_shared_arrays,
-		       rmin0,switchflag);
+		       rmin0,switchflag,bzeroflag);
     if (!use_shared_arrays)
       sna[tid]->grow_rij(nmax);
   }
@@ -1635,6 +1635,7 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
   rmin0 = 0.0;
   diagonalstyle = 3;
   switchflag = 1;
+  bzeroflag = 0;
   // open SNAP parameter file on proc 0
 
   FILE *fpparam;
@@ -1697,6 +1698,8 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
       diagonalstyle = atoi(keyval);
     else if (strcmp(keywd,"switchflag") == 0)
       switchflag = atoi(keyval);
+    else if (strcmp(keywd,"bzeroflag") == 0)
+      bzeroflag = atoi(keyval);
     else
       error->all(FLERR,"Incorrect SNAP parameter file");
   }
diff --git a/src/SNAP/pair_snap.h b/src/SNAP/pair_snap.h
index a6395bfd6a..559d3ef571 100644
--- a/src/SNAP/pair_snap.h
+++ b/src/SNAP/pair_snap.h
@@ -98,7 +98,7 @@ protected:
   double *wjelem;               // elements weights
   double **coeffelem;           // element bispectrum coefficients
   int *map;                     // mapping from atom types to elements
-  int twojmax, diagonalstyle, switchflag;
+  int twojmax, diagonalstyle, switchflag, bzeroflag;
   double rcutfac, rfac0, rmin0, wj1, wj2;
   int rcutfacflag, twojmaxflag; // flags for required parameters
   int gammaoneflag;              // 1 if parameter gamma is 1
diff --git a/src/SNAP/sna.cpp b/src/SNAP/sna.cpp
index 8b16b89336..2c20e78b71 100644
--- a/src/SNAP/sna.cpp
+++ b/src/SNAP/sna.cpp
@@ -115,14 +115,15 @@ using namespace MathConst;
 
 SNA::SNA(LAMMPS* lmp, double rfac0_in,
          int twojmax_in, int diagonalstyle_in, int use_shared_arrays_in,
-         double rmin0_in, int switch_flag_in) : Pointers(lmp)
+         double rmin0_in, int switch_flag_in, int bzero_flag_in) : Pointers(lmp)
 {
   wself = 1.0;
-
+  
   use_shared_arrays = use_shared_arrays_in;
   rfac0 = rfac0_in;
   rmin0 = rmin0_in;
   switch_flag = switch_flag_in;
+  bzero_flag = bzero_flag_in;
 
   twojmax = twojmax_in;
   diagonalstyle = diagonalstyle_in;
@@ -142,6 +143,12 @@ SNA::SNA(LAMMPS* lmp, double rfac0_in,
   nmax = 0;
   idxj = NULL;
 
+  if (bzero_flag) {
+    double www = wself*wself*wself;
+    for(int j = 0; j <= twojmax; j++)
+      bzero[j] = www*(j+1);
+  }
+  
 #ifdef TIMING_INFO
   timers = new double[20];
   for(int i = 0; i < 20; i++) timers[i] = 0;
@@ -151,6 +158,7 @@ SNA::SNA(LAMMPS* lmp, double rfac0_in,
 
   build_indexlist();
 
+  
 }
 
 /* ---------------------------------------------------------------------- */
@@ -539,13 +547,11 @@ void SNA::compute_bi()
           j <= MIN(twojmax, j1 + j2); j += 2) {
         barray[j1][j2][j] = 0.0;
 
-	for(int mb = 0; 2*mb < j; mb++) {
-	  for(int ma = 0; ma <= j; ma++) {
+	for(int mb = 0; 2*mb < j; mb++)
+	  for(int ma = 0; ma <= j; ma++)
             barray[j1][j2][j] +=
               uarraytot_r[j][ma][mb] * zarray_r[j1][j2][j][ma][mb] +
 	      uarraytot_i[j][ma][mb] * zarray_i[j1][j2][j][ma][mb];
-	  }
-	}
 
 	// For j even, special treatment for middle column
 
@@ -562,6 +568,8 @@ void SNA::compute_bi()
 	}
 
         barray[j1][j2][j] *= 2.0;
+	if (bzero_flag)
+	  barray[j1][j2][j] -= bzero[j];
       }
     }
 
@@ -1512,6 +1520,12 @@ void SNA::create_twojmax_arrays()
   memory->create(uarray_i, jdim, jdim, jdim,
                  "sna:uarray");
 
+  if (bzero_flag)
+    memory->create(bzero, jdim,"sna:bzero");
+  else
+    bzero = NULL;
+  
+
   if(!use_shared_arrays) {
     memory->create(uarraytot_r, jdim, jdim, jdim,
                    "sna:uarraytot");
@@ -1541,6 +1555,9 @@ void SNA::destroy_twojmax_arrays()
   memory->destroy(uarray_r);
   memory->destroy(uarray_i);
 
+  if (bzero_flag)
+    memory->destroy(bzero);
+
   if(!use_shared_arrays) {
     memory->destroy(uarraytot_r);
     memory->destroy(zarray_r);
diff --git a/src/SNAP/sna.h b/src/SNAP/sna.h
index c8bce915fb..d05ad0fb84 100644
--- a/src/SNAP/sna.h
+++ b/src/SNAP/sna.h
@@ -31,7 +31,7 @@ struct SNA_LOOPINDICES {
 class SNA : protected Pointers {
 
 public:
-  SNA(LAMMPS*, double, int, int, int, double, int);
+  SNA(LAMMPS*, double, int, int, int, double, int, int);
 
   SNA(LAMMPS* lmp) : Pointers(lmp) {};
   ~SNA();
@@ -139,6 +139,8 @@ private:
   // Self-weight
   double wself;
 
+  int bzero_flag; // 1 if bzero subtracted from barray
+  double *bzero;  // array of B values for isolated atoms
 };
 
 }
diff --git a/src/balance.cpp b/src/balance.cpp
index 52f6072a6c..47e7c0969b 100644
--- a/src/balance.cpp
+++ b/src/balance.cpp
@@ -256,6 +256,7 @@ void Balance::command(int narg, char **arg)
   // insure particles are in current box & update box via shrink-wrap
   // init entire system since comm->setup is done
   // comm::init needs neighbor::init needs pair::init needs kspace::init, etc
+  // must reset atom map after exchange() since it clears it
 
   MPI_Barrier(world);
   double start_time = MPI_Wtime();
@@ -267,6 +268,7 @@ void Balance::command(int narg, char **arg)
   domain->reset_box();
   comm->setup();
   comm->exchange();
+  if (atom->map_style) atom->map_set();
   if (domain->triclinic) domain->lamda2x(atom->nlocal);
 
   // imbinit = initial imbalance
diff --git a/src/comm.cpp b/src/comm.cpp
index b558b3fd8c..871675ca8d 100644
--- a/src/comm.cpp
+++ b/src/comm.cpp
@@ -266,7 +266,8 @@ void Comm::modify_params(int narg, char **arg)
     } else if (strcmp(arg[iarg],"cutoff") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal comm_modify command");
       if (mode == MULTI)
-        error->all(FLERR,"Use cutoff/multi keyword to set cutoff in multi mode");
+        error->all(FLERR,
+                   "Use cutoff/multi keyword to set cutoff in multi mode");
       cutghostuser = force->numeric(FLERR,arg[iarg+1]);
       if (cutghostuser < 0.0)
         error->all(FLERR,"Invalid cutoff in comm_modify command");
diff --git a/src/compute.cpp b/src/compute.cpp
index d1db42a24e..00a3984aab 100644
--- a/src/compute.cpp
+++ b/src/compute.cpp
@@ -127,11 +127,15 @@ void Compute::modify_params(int narg, char **arg)
 
   int iarg = 0;
   while (iarg < narg) {
-    if (strcmp(arg[iarg],"extra") == 0) {
+    // added more specific keywords in Mar17
+    // at some point, remove generic extra and dynamic
+    if (strcmp(arg[iarg],"extra") == 0 || 
+        strcmp(arg[iarg],"extra/dof") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal compute_modify command");
       extra_dof = force->numeric(FLERR,arg[iarg+1]);
       iarg += 2;
-    } else if (strcmp(arg[iarg],"dynamic") == 0) {
+    } else if (strcmp(arg[iarg],"dynamic") == 0 || 
+               strcmp(arg[iarg],"dynamic/dof") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal compute_modify command");
       if (strcmp(arg[iarg+1],"no") == 0) dynamic_user = 0;
       else if (strcmp(arg[iarg+1],"yes") == 0) dynamic_user = 1;
diff --git a/src/fix.cpp b/src/fix.cpp
index 0a95bcc692..ca9a69606c 100644
--- a/src/fix.cpp
+++ b/src/fix.cpp
@@ -71,6 +71,7 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) :
   restart_pbc = 0;
   wd_header = wd_section = 0;
   dynamic_group_allow = 0;
+  dynamic = 0;
   dof_flag = 0;
   special_alter_flag = 0;
   enforce2d_flag = 0;
@@ -130,7 +131,13 @@ void Fix::modify_params(int narg, char **arg)
 
   int iarg = 0;
   while (iarg < narg) {
-    if (strcmp(arg[iarg],"energy") == 0) {
+    if (strcmp(arg[iarg],"dynamic/dof") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command");
+      if (strcmp(arg[iarg+1],"no") == 0) dynamic = 0;
+      else if (strcmp(arg[iarg+1],"yes") == 0) dynamic = 1;
+      else error->all(FLERR,"Illegal fix_modify command");
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"energy") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command");
       if (strcmp(arg[iarg+1],"no") == 0) thermo_energy = 0;
       else if (strcmp(arg[iarg+1],"yes") == 0) thermo_energy = 1;
diff --git a/src/fix.h b/src/fix.h
index 8005da1add..d91937848a 100644
--- a/src/fix.h
+++ b/src/fix.h
@@ -217,6 +217,8 @@ class Fix : protected Pointers {
   int copymode;   // if set, do not deallocate during destruction
                   // required when classes are used as functors by Kokkos
 
+  int dynamic;    // recount atoms for temperature computes
+
   void ev_setup(int, int);
   void ev_tally(int, int *, double, double, double *);
   void v_setup(int);
diff --git a/src/fix_momentum.h b/src/fix_momentum.h
index 05fd7ff7c8..6e17f75a7f 100644
--- a/src/fix_momentum.h
+++ b/src/fix_momentum.h
@@ -34,7 +34,6 @@ class FixMomentum : public Fix {
  protected:
   int linear,angular,rescale;
   int xflag,yflag,zflag;
-  int dynamic;
   double masstotal;
 };
 
diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp
index f4a0a71a4a..0f1f3395dd 100644
--- a/src/fix_nh.cpp
+++ b/src/fix_nh.cpp
@@ -811,7 +811,6 @@ void FixNH::setup(int vflag)
           (etap_mass[ich-1]*etap_dot[ich-1]*etap_dot[ich-1] -
            boltz * t_target) / etap_mass[ich];
     }
-
   }
 }
 
diff --git a/src/library.cpp b/src/library.cpp
index da491f7152..ca3276ee8e 100644
--- a/src/library.cpp
+++ b/src/library.cpp
@@ -862,7 +862,8 @@ void lammps_scatter_atoms(void *ptr, char *name,
     int i,j,m,offset;
     void *vptr = lmp->atom->extract(name);
     if(vptr == NULL) {
-        lmp->error->warning(FLERR,"lammps_scatter_atoms: unknown property name");
+        lmp->error->warning(FLERR,
+                            "lammps_scatter_atoms: unknown property name");
         return;
     }
 
diff --git a/src/special.cpp b/src/special.cpp
index c4de11e09c..3fb5ec8077 100644
--- a/src/special.cpp
+++ b/src/special.cpp
@@ -75,8 +75,8 @@ void Special::build()
     const double * const special_lj   = force->special_lj;
     const double * const special_coul = force->special_coul;
     fprintf(screen,"Finding 1-2 1-3 1-4 neighbors ...\n"
-                   " Special bond factors lj:   %-10g %-10g %-10g\n"
-                   " Special bond factors coul: %-10g %-10g %-10g\n",
+                   "  special bond factors lj:   %-10g %-10g %-10g\n"
+                   "  special bond factors coul: %-10g %-10g %-10g\n",
                    special_lj[1],special_lj[2],special_lj[3],
                    special_coul[1],special_coul[2],special_coul[3]);
   }
@@ -498,6 +498,7 @@ void Special::dedup()
 
   // re-create map
 
+  atom->map_init(0);
   atom->nghost = 0;
   atom->map_set();
 }
@@ -649,6 +650,7 @@ void Special::combine()
 
   // re-create map
 
+  atom->map_init(0);
   atom->nghost = 0;
   atom->map_set();
 }