diff --git a/doc/src/Section_howto.txt b/doc/src/Section_howto.txt
index aedc2fc64fff9eb445f206f53c549423e8b29b14..e83f5399c0a116c086fd3df5bc09414c7398d4b0 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 2145cfcbce7d04441a723142236adb97edc6c3f5..b994a564027e8b86b3230cdfc67c69f9b754c9c2 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
diff --git a/doc/src/compute_modify.txt b/doc/src/compute_modify.txt
index acf14526a71f64ca59755ac7060c9508fb613250..637f9b5e417da481cb173670cf3eb21c2a7fb929 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 0332ab9c1c51c7142cf5e8990ba0817665ad9e34..e2df706473e178202b1fbb10985cade08dce831f 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 a1dd5f6434d147ed6d1b2c2d855bd100dfb98534..477c14ea8955ed45763bb1dd27d1242898e22ff5 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 4eab458a65484b99a49776f5a8430cc91360c794..1a419628e9dc896fcd7f7b6eea12af8ab805008b 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 9c95cdc452e498ae3a3f2c51ce04ea97375f66b2..5e097bb34f6fde43a238b9ecfbdf9c50e6f3d98d 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 f2c10184ddb4f997a80d765e45d2a6ee1620fde3..54f78287e01d98368f1938122b1c9317e8446659 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 a24c6c3160ff6f3d337ae22a6ae9be591b121f0b..ab7313832cc0aa57ca4539e82466449c397d3754 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 66e3916fa48ed9c527f7321b07b101ef008f7d6a..4dbdbb8554f53942cfebd6bbcc2c583d536b1246 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 8f29040399cd1bed99e6a4f86cfcdc1a443514f0..3e48276bb214015d4c0e5e9b81d68d1142e717a0 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 870d740724d1d45fe56963d9061d6644715ef628..16b0b9a72ae59e4ac394913247e2f3c525068a55 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 a1c141088b93ded8f957046a008663ce26d3c3b7..7bd701c9243c973e38e6b33418906c10aaf49a84 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 19e4e2414af4ec4ce0cd4b2958b8ff065aa09559..fe783ac6de7140fc615a5e9b4458ac7bc583069b 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 df52aef66cb67081b51813dfa63c1b03c77f45d4..82398987f809a7ac4eb96ac1d55c2817a7f2a886 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 91f06744d35ada733c855c1a29a09013542f056c..f1f20ad79cdf636670f08d0f3aa5faa44a23b3d0 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 9b24b2906c1f1ffb509fa7ebc241267ae4ef7735..332ccf8b62dd5d472518639e153211c4034ae26e 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 02fc794ffdcbe3db9829d008d5f77bda95b2d8bf..5fa057c746c8687dd48631a05150eece251d1d9d 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 84869eb233a7073ee0b84094ebe925013822125c..a65c72875f0537716e86b7707bfe843e62cb6774 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 0627253341a4c134b4a13a9749a60abcb7bc3027..7b30312f59abbbb1759580e35f0279fdf1cc1c98 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 e0b20005e6a50d4a264b91995cbb1ffe96d2fcd7..f17961bdd737f379302feca3b87a1f1c60690226 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 e638cd887365100e6b96d462f06759acc3eff6aa..780df25bdf7ab95c3dd74c069012165fd52da4f2 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 8261f6e38c92137f8682ce890dc12490d191a30c..9b2184dda2810f1dea8ddf1de0a183b9773bc46d 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 d841482f8e99a4715ec46d22eb5e839888d6f732..9c1082f8161685714fa8525c6797fef9c40d177b 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 b08c80b626bb59752c7656dd814c3c4e0dc45848..199c04bd9c483938c2d90463f5ee0737dce7f64f 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 07510a59e3b486bd8e993c2711a415e807029acb..6b51933277bb2bcbe0d46f33a69cce0169829897 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 326d2d620d4d9d28379fd6e14ed07ff307e03198..ad934535abd560c0b405109b17d4ce1d53a8c814 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 efd4cafbcaed1f0f1bb34ac20e07d43e994dc4ac..73452427bd48e42810f872eb6d4ba1c5aea4a696 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 c1398864e0e5c12552fb804a24bec53bff0e9752..f75b02fba72ea5f7fef25edc992c60c1cfc3ab4b 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 dc84b0be055a41ba2dfd400e605aadd803256b12..06c2e48488e15d8ae072e5a7b16faafc0ed8a1e1 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 a6395bfd6a3d308a255ca17107149b36bb929697..559d3ef5717f46d0f75bb0518010b46dd30700c8 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 8b16b89336d35aadf4d5452e4e7d6e677457bf29..2c20e78b71b2cba9d26553ff90f29fd46790ac76 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 c8bce915fb5d6476e67a9f6c0787566a80f048e9..d05ad0fb8460b0d517752c1326e594f256ed41d6 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 52f6072a6c014a229ab05973927bf4e6b106052d..47e7c0969b15f6917475af5020879ed0220dbb66 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 b558b3fd8c0e0d4b37e4fb3b0b87a5e332089105..871675ca8dc52ea1259a58a09f0afb869f168a71 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 d1db42a24eb27472cf2b8c2db1b1c9bcb22ee497..00a3984aab3c668b776216e12b8dd7871b1d96de 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 0a95bcc6920d82c929f1c9fca7348e67c95fb4dd..ca9a69606c375e01032a05eaebd602f760ffe77e 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 8005da1addca8956991c283111901016f551d544..d91937848ac8ad4f42b0d86973a4afc3a374518f 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 05fd7ff7c886203be427fcf2dc71007b46b050a0..6e17f75a7f59b86fb90d511cde3356267126ddd6 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 f4a0a71a4a2e8d0b944ee9b05e29dc9a69f9b15f..0f1f3395dd51d8a5bcddb7fe3ccbcb63c454f055 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 da491f71528a319a69a2e28f822ac72e6ad14399..ca3276ee8ed59e486ce71924bf3ced417c147c26 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 c4de11e09c5f9232573444b90261bbb8995fd00c..3fb5ec8077e21130009981b82f69958131072fe7 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();
}