fix ave/chunk fixes, 2d disc option, fix_modify dynamic/dof

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
+the extract_atom() method, see the extract() method in the
-as LAMMPS runs, so you need to re-call the function to assure a
+src/atom.cpp file for a list of valid per-atom properties.  New names
-current pointer or returned value(s).
+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,
+gather function collects peratom info of the requested type (atom
-types, forces, etc) from all processors, orders them by atom ID, and
+coords, types, forces, etc) from all processors, orders them by atom
-returns a full list to each calling processor.  The scatter function
+ID, and returns a full list to each calling processor.  The scatter
-does the inverse.  It distributes the same kinds of values, 
+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

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
+Again, new names could easily be added if the property you want is not
-doubles or integers, or a pointer to an array of doubles (double **)
+listed.  A pointer to a vector of doubles or integers, or a pointer to
-or integers (int **) is returned.  You need to specify the appropriate
+an array of doubles (double **) or integers (int **) is returned.  You
-data type via the type argument.
+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
+The gather_atoms() method allows any per-atom property (coordinates,
-as specified by type, of length count*natoms, for the property of all
+velocities, etc) to be extracted from LAMMPS.  It returns a ctypes
-the atoms in the simulation specified by name, ordered by count and
+vector of ints or doubles as specified by type, of length
-then by atom ID.  The vector can be used via normal Python
+count*natoms, for the named property for all atoms in the simulation.
-subscripting.  If atom IDs are not consecutively ordered within
+The data is ordered by count and then by atom ID.  See the extract()
-LAMMPS, a None is returned as indication of an error.
+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
+The scatter_atoms() method allows any per-atom property (coordinates,
-specified by type, of length count*natoms, for the property of all the
+velocities, etc) to be inserted into LAMMPS, overwriting the current
-atoms in the simulation specified by name, ordered by bount and then
+property.  It takes a vector of ints or doubles as specified by type,
-by atom ID.  It uses the vector of data to overwrite the corresponding
+of length count*natoms, for the named property for all atoms in the
-properties for each atom inside LAMMPS.  This requires LAMMPS to have
+simulation.  The data should be ordered by count and then by atom ID.
-its "map" option enabled; see the "atom_modify"_atom_modify.html
+See the extract() method in the src/atom.cpp file for a list of valid
-command for details.  If it is not, or if atom IDs are not
+names.  Again, new names could easily be added if the property you
-consecutively ordered, no coordinates are reset.
+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.

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
+keyword = {extra/dof} or {extra} or {dynamic/dof} or {dynamic} :l
-  {extra} value = N
+  {extra/dof} value = N
     N = # of extra degrees of freedom to subtract
-  {dynamic} value = {yes} or {no}
+  {extra} syntax is identical to {extra/dof}, will be disabled at some point
-    yes/no = do or do not recompute the number of atoms contributing to the temperature :pre
+  {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 myTemp extra/dof 0
-compute_modify newtemp dynamic yes extra 600 :pre
+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
+The {extra/dof} or {extra} keyword refers to how many
-subtracted (typically from 3N) as a normalizing factor in a
+degrees-of-freedom are subtracted (typically from 3N) as a normalizing
-temperature computation.  Only computes that compute a temperature use
+factor in a temperature computation.  Only computes that compute a
-this option.  The default is 2 or 3 for "2d or 3d
+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
+The {dynamic/dof} or {dynamic} keyword determines whether the number
-compute group is re-computed each time a temperature is computed.
+of atoms N in the compute group is re-computed each time a temperature
-Only compute styles that calculate a temperature use this option.  By
+is computed.  Only compute styles that calculate a temperature use
-default, N is assumed to be constant.  If you are adding atoms to the
+this option.  By default, N is assumed to be constant.  If you are
-system (see the "fix pour"_fix_pour.html or "fix
+adding atoms to the system (see the "fix pour"_fix_pour.html, "fix
-deposit"_fix_deposit.html commands) or expect atoms to be lost
+deposit"_fix_deposit.html and "fix gcmc"_fix_gcmc.html commands) or
-(e.g. due to evaporation), then this option should be used to insure
+expect atoms to be lost (e.g. due to evaporation), then this option
-the temperature is correctly normalized.
+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
+The option defaults are extra/dof = 2 or 3 for 2d or 3d systems and
-dynamic = no.
+dynamic/dof = no.

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
 

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

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,

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:]

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

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.

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

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

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

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

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
 

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
 

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

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
 

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
 

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
 

examples/snap/Ta06A.snapparam

@@ -12,3 +12,4 @@ gamma 1
 rfac0 0.99363
 rmin0 0
 diagonalstyle 3
+bzeroflag 0

examples/snap/W_2940_2017_2.snapparam

@@ -10,3 +10,4 @@ gamma 1
 rfac0 0.99363
 rmin0 0
 diagonalstyle 3
+bzeroflag 0

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");
 
@@ -1354,7 +1375,12 @@ void FixGCMC::attempt_molecule_insertion()
       }
     }
 
-    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;
@@ -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;

src/MC/fix_gcmc.h

@@ -128,9 +128,9 @@ class FixGCMC : public Fix {
   int imol,nmol;
   double **coords;
   imageint *imageflags;
-  class Fix *fixshake;
+  class Fix *fixrigid, *fixshake;
-  int shakeflag;
+  int rigidflag, shakeflag;
-  char *idshake;
+  char *idrigid, *idshake;
   int triclinic;                         // 0 = orthog box, 1 = triclinic
 
   class Compute *c_pe;

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))

src/RIGID/fix_rigid_nh_small.cpp

@@ -321,35 +321,7 @@ void FixRigidNHSmall::init()
 void FixRigidNHSmall::setup(int vflag)
 {
   FixRigidSmall::setup(vflag);
-
+  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;
   
   double mbody[3];
   akin_t = akin_r = 0.0;
@@ -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;
 }
-

src/RIGID/fix_rigid_nh_small.h

@@ -78,6 +78,7 @@ 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();

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;

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,6 +48,7 @@ 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
@@ -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;

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;

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");
   }

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

src/SNAP/sna.cpp

@@ -115,7 +115,7 @@ 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;
   
@@ -123,6 +123,7 @@ SNA::SNA(LAMMPS* lmp, double rfac0_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 mb = 0; 2*mb < j; mb++)
-	  for(int ma = 0; ma <= j; ma++) {
+	  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);

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
 };
 
 }

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

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");

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;

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;

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);

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;
 };
 

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];
     }
-
   }
 }
 

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;
     }
 

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 lj:   %-10g %-10g %-10g\n"
-                   " Special bond factors coul: %-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();
 }