diff --git a/doc/src/Commands_compute.txt b/doc/src/Commands_compute.txt
index 028e274c9be4a80555dfcab5a417d738b763bfda..6185634187fb0e5089bf9759b221aae5ee53a53a 100644
--- a/doc/src/Commands_compute.txt
+++ b/doc/src/Commands_compute.txt
@@ -35,6 +35,7 @@ KOKKOS, o = USER-OMP, t = OPT.
"bond/local"_compute_bond_local.html,
"centro/atom"_compute_centro_atom.html,
"chunk/atom"_compute_chunk_atom.html,
+"chunk/spread/atom"_compute_chunk_spread_atom.html,
"cluster/atom"_compute_cluster_atom.html,
"cna/atom"_compute_cna_atom.html,
"cnp/atom"_compute_cnp_atom.html,
@@ -97,6 +98,7 @@ KOKKOS, o = USER-OMP, t = OPT.
"property/local"_compute_property_local.html,
"rdf"_compute_rdf.html,
"reduce"_compute_reduce.html,
+"reduce/chunk"_compute_reduce_chunk.html,
"reduce/region"_compute_reduce.html,
"rigid/local"_compute_rigid_local.html,
"saed"_compute_saed.html,
diff --git a/doc/src/Howto_chunk.txt b/doc/src/Howto_chunk.txt
index 8e52acf4b8693569bac54304bb5f2c0ffd56d0fb..611b71dac7d4546ad85241ac3d0b1786a05aff2b 100644
--- a/doc/src/Howto_chunk.txt
+++ b/doc/src/Howto_chunk.txt
@@ -22,7 +22,7 @@ commands, to calculate various properties of a system:
"fix ave/chunk"_fix_ave_chunk.html
any of the "compute */chunk"_compute.html commands :ul
-Here, each of the 3 kinds of chunk-related commands is briefly
+Here, each of the 4 kinds of chunk-related commands is briefly
overviewed. Then some examples are given of how to compute different
properties with chunk commands.
@@ -83,8 +83,9 @@ chunk.
Compute */chunk commands: :h4
-Currently the following computes operate on chunks of atoms to produce
-per-chunk values.
+The following computes operate on chunks of atoms to produce per-chunk
+values. Any compute whose style name ends in "/chunk" is in this
+category:
"compute com/chunk"_compute_com_chunk.html
"compute gyration/chunk"_compute_gyration_chunk.html
@@ -111,8 +112,8 @@ of a center of mass, which requires summing mass*position over the
atoms and then dividing by summed mass.
All of these computes produce a global vector or global array as
-output, wih one or more values per chunk. They can be used
-in various ways:
+output, wih one or more values per chunk. The output can be used in
+various ways:
As input to the "fix ave/time"_fix_ave_time.html command, which can
write the values to a file and optionally time average them. :ulb,l
@@ -122,9 +123,27 @@ histogram values across chunks. E.g. a histogram of cluster sizes or
molecule diffusion rates. :l
As input to special functions of "equal-style
-variables"_variable.html, like sum() and max(). E.g. to find the
-largest cluster or fastest diffusing molecule. :l
-:ule
+variables"_variable.html, like sum() and max() and ave(). E.g. to
+find the largest cluster or fastest diffusing molecule or average
+radius-of-gyration of a set of molecules (chunks). :l :ule
+
+Other chunk commands: :h4
+
+"compute chunk/spread/atom"_compute_chunk_spread_atom.html
+"compute reduce/chunk"_compute_reduce_chunk.html :ul
+
+The "compute chunk/spread/atom"_compute_chunk_spread_atom.html command
+spreads per-chunk values to each atom in the chunk, producing per-atom
+values as its output. This can be useful for outputting per-chunk
+values to a per-atom "dump file"_dump.html. Or for using an atom's
+associated chunk value in an "atom-style variable"_variable.html.
+
+The "compute reduce/chunk"_compute_reduce_chunk.html command reduces a
+peratom value across the atoms in each chunk to produce a value per
+chunk. When used with the "compute
+chunk/spread/atom"_compute_chunk_spread_atom.html command it can
+create peratom values that induce a new set of chunks with a second
+"compute chunk/atom"_compute_chunk_atom.html command.
Example calculations with chunks :h4
@@ -164,3 +183,13 @@ compute cluster all cluster/atom 1.0
compute cc1 all chunk/atom c_cluster compress yes
compute size all property/chunk cc1 count
fix 1 all ave/histo 100 1 100 0 20 20 c_size mode vector ave running beyond ignore file tmp.histo :pre
+
+(6) An example of using a per-chunk value to apply per-atom forces to
+compress individual polymer chains (molecules) in a mixture, is
+explained on the "compute
+chunk/spread/atom"_compute_chunk_spread_atom.html command doc page.
+
+(7) An example of using one set of per-chunk values for molecule
+chunks, to create a 2nd set of micelle-scale chunks (clustered
+molecules, due to hydrophobicity), is explained on the "compute
+chunk/reduce"_compute_reduce_chunk.html command doc page.
diff --git a/doc/src/compute.txt b/doc/src/compute.txt
index 8facb4de63ebfb527e11108ba02e381db6a1425f..72ea1c193004b132a32af49797e12a58478bdaad 100644
--- a/doc/src/compute.txt
+++ b/doc/src/compute.txt
@@ -183,6 +183,7 @@ compute"_Commands_compute.html doc page are followed by one or more of
"bond/local"_compute_bond_local.html - distance and energy of each bond
"centro/atom"_compute_centro_atom.html - centro-symmetry parameter for each atom
"chunk/atom"_compute_chunk_atom.html - assign chunk IDs to each atom
+"chunk/spread/atom"_compute_chunk_spread_atom.html - spreads chunk values to each atom in chunk
"cluster/atom"_compute_cluster_atom.html - cluster ID for each atom
"cna/atom"_compute_cna_atom.html - common neighbor analysis (CNA) for each atom
"com"_compute_com.html - center-of-mass of group of atoms
@@ -225,6 +226,7 @@ compute"_Commands_compute.html doc page are followed by one or more of
"property/chunk"_compute_property_chunk.html - extract various per-chunk attributes
"rdf"_compute_rdf.html - radial distribution function g(r) histogram of group of atoms
"reduce"_compute_reduce.html - combine per-atom quantities into a single global value
+"reduce/chunk"_compute_reduce_chunk.html - reduce per-atom quantities within each chunk
"reduce/region"_compute_reduce.html - same as compute reduce, within a region
"rigid/local"_compute_rigid_local.html - extract rigid body attributes
"slice"_compute_slice.html - extract values from global vector or array
diff --git a/doc/src/compute_chunk_atom.txt b/doc/src/compute_chunk_atom.txt
index e76b51e6ec9df377f4e1354f0090c588f6dd4b53..c29a5600a982659325cd6af6948eb8ed16de22d9 100644
--- a/doc/src/compute_chunk_atom.txt
+++ b/doc/src/compute_chunk_atom.txt
@@ -14,7 +14,7 @@ compute ID group-ID chunk/atom style args keyword values ... :pre
ID, group-ID are documented in "compute"_compute.html command :ulb,l
chunk/atom = style name of this compute command :l
-style = {bin/1d} or {bin/2d} or {bin/3d} or {bin/sphere} or {type} or {molecule} or {compute/fix/variable}
+style = {bin/1d} or {bin/2d} or {bin/3d} or {bin/sphere} or {type} or {molecule} or c_ID, c_ID\[I\], f_ID, f_ID\[I\], v_name
{bin/1d} args = dim origin delta
dim = {x} or {y} or {z}
origin = {lower} or {center} or {upper} or coordinate value (distance units)
@@ -40,7 +40,7 @@ style = {bin/1d} or {bin/2d} or {bin/3d} or {bin/sphere} or {type} or {molecule}
ncbin = # of concentric circle bins between rmin and rmax
{type} args = none
{molecule} args = none
- {compute/fix/variable} = c_ID, c_ID\[I\], f_ID, f_ID\[I\], v_name with no args
+ c_ID, c_ID\[I\], f_ID, f_ID\[I\], v_name args = none
c_ID = per-atom vector calculated by a compute with ID
c_ID\[I\] = Ith column of per-atom array calculated by a compute with ID
f_ID = per-atom vector calculated by a fix with ID
@@ -85,7 +85,8 @@ compute 1 all chunk/atom bin/1d z lower 0.02 units reduced
compute 1 all chunk/atom bin/2d z lower 1.0 y 0.0 2.5
compute 1 all chunk/atom molecule region sphere nchunk once ids once compress yes
compute 1 all chunk/atom bin/sphere 5 5 5 2.0 5.0 5 discard yes
-compute 1 all chunk/atom bin/cylinder z lower 2 10 10 2.0 5.0 3 discard yes :pre
+compute 1 all chunk/atom bin/cylinder z lower 2 10 10 2.0 5.0 3 discard yes
+compute 1 all chunk/atom c_cluster :pre
[Description:]
@@ -386,8 +387,8 @@ described below, which resets {Nchunk}. The {limit} keyword is then
applied to the new {Nchunk} value, exactly as described in the
preceding paragraph. Note that in this case, all atoms will end up
with chunk IDs <= {Nc}, but their original values (e.g. molecule ID or
-compute/fix/variable value) may have been > {Nc}, because of the
-compression operation.
+compute/fix/variable) may have been > {Nc}, because of the compression
+operation.
If {compress yes} is set, and the {compress} keyword comes after the
{limit} keyword, then the {limit} value of {Nc} is applied first to
diff --git a/doc/src/compute_chunk_spread_atom.txt b/doc/src/compute_chunk_spread_atom.txt
new file mode 100644
index 0000000000000000000000000000000000000000..59662035e531479a36ac73e116abc31a6e9924ff
--- /dev/null
+++ b/doc/src/compute_chunk_spread_atom.txt
@@ -0,0 +1,173 @@
+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Commands_all.html)
+
+:line
+
+compute chunk/spread/atom command :h3
+
+[Syntax:]
+
+compute ID group-ID chunk/spread/atom chunkID input1 input2 ... :pre
+
+ID, group-ID are documented in "compute"_compute.html command :ulb,l
+chunk/spread/atom = style name of this compute command :l
+chunkID = ID of "compute chunk/atom"_compute_chunk_atom.html command :l
+one or more inputs can be listed :l
+input = c_ID, c_ID\[N\], f_ID, f_ID\[N\] :l
+ c_ID = global vector calculated by a compute with ID
+ c_ID\[I\] = Ith column of global array calculated by a compute with ID, I can include wildcard (see below)
+ f_ID = global vector calculated by a fix with ID
+ f_ID\[I\] = Ith column of global array calculated by a fix with ID, I can include wildcard (see below) :pre
+:ule
+
+[Examples:]
+
+compute 1 all chunk/spread/atom mychunk c_com[*] c_gyration :pre
+
+[Description:]
+
+Define a calculation that "spreads" one or more per-chunk values to
+each atom in the chunk. This can be useful for creating a "dump
+file"_dump.html where each atom lists info about the chunk it is in,
+e.g. for post-processing purposes. It can also be used in "atom-style
+variables"_variable.html that need info about the chunk each atom is
+in. Examples are given below.
+
+In LAMMPS, chunks are collections of atoms defined by a "compute
+chunk/atom"_compute_chunk_atom.html command, which assigns each atom
+to a single chunk (or no chunk). The ID for this command is specified
+as chunkID. For example, a single chunk could be the atoms in a
+molecule or atoms in a spatial bin. See the "compute
+chunk/atom"_compute_chunk_atom.html and "Howto chunk"_Howto_chunk.html
+doc pages for details of how chunks can be defined and examples of how
+they can be used to measure properties of a system.
+
+For inputs that are computes, they must be a compute that calculates
+per-chunk values. These are computes whose style names end in
+"/chunk".
+
+For inputs that are fixes, they should be a a fix that calculates
+per-chunk values. For example, "fix ave/chunk"_fix_ave_chunk.html or
+"fix ave/time"_fix_ave_time.html (assuming it is time-averaging
+per-chunk data).
+
+For each atom, this compute accesses its chunk ID from the specified
+{chunkID} compute, then accesses the per-chunk value in each input.
+Those values are copied to this compute to become the output for that
+atom.
+
+The values generated by this compute will be 0.0 for atoms not in the
+specified compute group {group-ID}. They will also be 0.0 if the atom
+is not in a chunk, as assigned by the {chunkID} compute. They will
+also be 0.0 if the current chunk ID for the atom is out-of-bounds with
+respect to the number of chunks stored by a particular input compute
+or fix.
+
+NOTE: LAMMPS does not check that a compute or fix which calculates
+per-chunk values uses the same definition of chunks as this compute.
+It's up to you to be consistent. Likewise, for a fix input, LAMMPS
+does not check that it is per-chunk data. It only checks that the fix
+produces a global vector or array.
+
+:line
+
+Each listed input is operated on independently.
+
+If a bracketed index I is used, it can be specified using a wildcard
+asterisk with the index to effectively specify multiple values. This
+takes the form "*" or "*n" or "n*" or "m*n". If N = the number of
+columns in the array, then an asterisk with no numeric values means
+all indices from 1 to N. A leading asterisk means all indices from 1
+to n (inclusive). A trailing asterisk means all indices from n to N
+(inclusive). A middle asterisk means all indices from m to n
+(inclusive).
+
+Using a wildcard is the same as if the individual columns of the array
+had been listed one by one. E.g. these 2 compute chunk/spread/atom
+commands are equivalent, since the "compute
+com/chunk"_compute_com_chunk.html command creates a per-atom array
+with 3 columns:
+
+compute com all com/chunk mychunk
+compute 10 all chunk/spread/atom mychunk c_com\[*\]
+compute 10 all chunk/spread/atom mychunk c_com\[1\] c_com\[2\] c_com\[3\] :pre
+
+:line
+
+Here is an example of writing a dump file the with the center-of-mass
+(COM) for the chunk each atom is in:
+
+compute cmol all chunk/atom molecule
+compute com all com/chunk cmol
+compute comchunk all chunk/spread/atom cmol c_com[*]
+dump 1 all custom 50 tmp.dump id mol type x y z c_comchunk[*]
+dump_modify 1 sort id :pre
+
+The same per-chunk data for each atom could be used to define per-atom
+forces for the "fix addforce"_fix_addforce.html command. In this
+example the forces act to pull atoms of an extended polymer chain
+towards its COM in an attractive manner.
+
+compute prop all property/atom xu yu zu
+variable k equal 0.1
+variable fx atom v_k*(c_comchunk\[1\]-c_prop\[1\])
+variable fy atom v_k*(c_comchunk\[2\]-c_prop\[2\])
+variable fz atom v_k*(c_comchunk\[3\]-c_prop\[3\])
+fix 3 all addforce v_fx v_fy v_fz :pre
+
+Note that "compute property/atom"_compute_property_atom.html is used
+to generate unwrapped coordinates for use in the per-atom force
+calculation, so that the effect of periodic boundaries is accounted
+for properly.
+
+Over time this applied force could shrink each polymer chain's radius
+of gyration in a polymer mixture simulation. Here is output after
+adding the above lines to the bench/in.chain script. The thermo
+output is shown for 1000 steps, where the last column is the average
+radius of gyration over all 320 chains in the 32000 atom system:
+
+compute gyr all gyration/chunk cmol
+variable ave equal ave(c_gyr)
+thermo_style custom step etotal press v_ave :pre
+
+ 0 22.394765 4.6721833 5.128278
+ 100 22.445002 4.8166709 5.0348372
+ 200 22.500128 4.8790392 4.9364875
+ 300 22.534686 4.9183766 4.8590693
+ 400 22.557196 4.9492211 4.7937849
+ 500 22.571017 4.9161853 4.7412008
+ 600 22.573944 5.0229708 4.6931243
+ 700 22.581804 5.0541301 4.6440647
+ 800 22.584683 4.9691734 4.6000016
+ 900 22.59128 5.0247538 4.5611513
+ 1000 22.586832 4.94697 4.5238362 :pre
+
+:line
+
+[Output info:]
+
+This compute calculates a per-atom vector or array, which can be
+accessed by any command that uses per-atom values from a compute as
+input. See the "Howto output"_Howto_output.html doc page for an
+overview of LAMMPS output options.
+
+The output is a per-atom vector if a single input value is specified,
+otherwise a per-atom array is output. The number of columns in the
+array is the number of inputs provided. The per-atom values for the
+vector or each column of the array will be in whatever
+"units"_units.html the corresponding input value is in.
+
+The vector or array values are "intensive".
+
+[Restrictions:] none
+
+[Related commands:]
+
+"compute chunk/atom"_compute_chunk_atom.html, "fix
+ave/chunk"_fix_ave_chunk.html, "compute
+reduce/chunk"_compute_reduce_chunk.html
+
+[Default:] none
diff --git a/doc/src/compute_reduce.txt b/doc/src/compute_reduce.txt
index ef3c7c6489e03ef3c1056edefd8d4dc419766cdf..0bd2accf3c08ba9d7f50cb008c8b8151fe8b4012 100644
--- a/doc/src/compute_reduce.txt
+++ b/doc/src/compute_reduce.txt
@@ -97,9 +97,9 @@ equivalent, since the "compute stress/atom"_compute_stress_atom.html
command creates a per-atom array with 6 columns:
compute myPress all stress/atom NULL
-compute 2 all reduce min myPress\[*\]
-compute 2 all reduce min myPress\[1\] myPress\[2\] myPress\[3\] &
- myPress\[4\] myPress\[5\] myPress\[6\] :pre
+compute 2 all reduce min c_myPress\[*\]
+compute 2 all reduce min c_myPress\[1\] c_myPress\[2\] c_myPress\[3\] &
+ c_myPress\[4\] c_myPress\[5\] c_myPress\[6\] :pre
:line
diff --git a/doc/src/compute_reduce_chunk.txt b/doc/src/compute_reduce_chunk.txt
new file mode 100644
index 0000000000000000000000000000000000000000..3628c0f1ad7a0a0e1add2773d24ecf97d86b89ee
--- /dev/null
+++ b/doc/src/compute_reduce_chunk.txt
@@ -0,0 +1,177 @@
+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Commands_all.html)
+
+:line
+
+compute reduce/chunk command :h3
+
+[Syntax:]
+
+compute ID group-ID reduce/chunk chunkID mode input1 input2 ... :pre
+
+ID, group-ID are documented in "compute"_compute.html command :ulb,l
+reduce/chunk = style name of this compute command :l
+chunkID = ID of "compute chunk/atom"_compute_chunk_atom.html command :l
+mode = {sum} or {min} or {max} :l
+one or more inputs can be listed :l
+input = c_ID, c_ID\[N\], f_ID, f_ID\[N\], v_ID :l
+ c_ID = per-atom vector calculated by a compute with ID
+ c_ID\[I\] = Ith column of per-atom array calculated by a compute with ID, I can include wildcard (see below)
+ f_ID = per-atom vector calculated by a fix with ID
+ f_ID\[I\] = Ith column of per-atom array calculated by a fix with ID, I can include wildcard (see below)
+ v_name = per-atom vector calculated by an atom-style variable with name :pre
+:ule
+
+[Examples:]
+
+compute 1 all reduce/chunk/atom mychunk min c_cluster :pre
+
+[Description:]
+
+Define a calculation that reduces one or more per-atom vectors into
+per-chunk values. This can be useful for diagnostic output. Or when
+used in conjunction with the "compute
+chunk/spread/atom"_compute_chunk_spread_atom.html command it can be
+used ot create per-atom values that induce a new set of chunks with a
+second "compute chunk/atom"_compute_chunk_atom.html command. An
+example is given below.
+
+In LAMMPS, chunks are collections of atoms defined by a "compute
+chunk/atom"_compute_chunk_atom.html command, which assigns each atom
+to a single chunk (or no chunk). The ID for this command is specified
+as chunkID. For example, a single chunk could be the atoms in a
+molecule or atoms in a spatial bin. See the "compute
+chunk/atom"_compute_chunk_atom.html and "Howto chunk"_Howto_chunk.html
+doc pages for details of how chunks can be defined and examples of how
+they can be used to measure properties of a system.
+
+For each atom, this compute accesses its chunk ID from the specified
+{chunkID} compute. The per-atom value from an input contributes
+to a per-chunk value corresponding the the chunk ID.
+
+The reduction operation is specified by the {mode} setting and is
+performed over all the per-atom values from the atoms in each chunk.
+The {sum} option adds the pre-atom values to a per-chunk total. The
+{min} or {max} options find the minimum or maximum value of the
+per-atom values for each chunk.
+
+Note that only atoms in the specified group contribute to the
+reduction operation. If the {chunkID} compute returns a 0 for the
+chunk ID of an atom (i.e. the atom is not in a chunk defined by the
+"compute chunk/atom"_compute_chunk_atom.html command), that atom will
+also not contribute to the reduction operation. An input that is a
+compute or fix may define its own group which affects the quantities
+it returns. For example, a compute with return a zero value for atoms
+that are not in the group specified for that compute.
+
+Each listed input is operated on independently. Each input can be the
+result of a "compute"_compute.html or "fix"_fix.html or the evaluation
+of an atom-style "variable"_variable.html.
+
+Note that for values from a compute or fix, the bracketed index I can
+be specified using a wildcard asterisk with the index to effectively
+specify multiple values. This takes the form "*" or "*n" or "n*" or
+"m*n". If N = the size of the vector (for {mode} = scalar) or the
+number of columns in the array (for {mode} = vector), then an asterisk
+with no numeric values means all indices from 1 to N. A leading
+asterisk means all indices from 1 to n (inclusive). A trailing
+asterisk means all indices from n to N (inclusive). A middle asterisk
+means all indices from m to n (inclusive).
+
+Using a wildcard is the same as if the individual columns of the array
+had been listed one by one. E.g. these 2 compute reduce/chunk
+commands are equivalent, since the "compute
+property/chunk"_compute_property_chunk.html command creates a per-atom
+array with 3 columns:
+
+compute prop all property/atom vx vy vz
+compute 10 all reduce/chunk mychunk max c_prop\[*\]
+compute 10 all reduce/chunk mychunk max c_prop\[1\] c_prop\[2\] c_prop\[3\] :pre
+
+:line
+
+Here is an example of using this compute, in conjunction with the
+compute chunk/spread/atom command to identify self-assembled micelles.
+The commands below can be added to the examples/in.micelle script.
+
+Imagine a collection of polymer chains or small molecules with
+hydrophobic end groups. All the hydrophobic (HP) atoms are assigned
+to a group called "phobic".
+
+These commands will assign a unique cluster ID to all HP atoms within
+a specified distance of each other. A cluster will contain all HP
+atoms in a single molecule, but also the HP atoms in nearby molecules,
+e.g. molecules that have clumped to form a micelle due to the
+attraction induced by the hydrophobicity. The output of the
+chunk/reduce command will be a cluster ID per chunk (molecule).
+Molecules with the same cluster ID are in the same micelle.
+
+group phobic type 4 # specific to in.micelle model
+compute cluster phobic cluster/atom 2.0
+compute cmol all chunk/atom molecule
+compute reduce phobic reduce/chunk cmol min c_cluster :pre
+
+This per-chunk info could be output in at least two ways:
+
+fix 10 all ave/time 1000 1 1000 c_reduce file tmp.phobic mode vector :pre
+
+compute spread all chunk/spread/atom cmol c_reduce
+dump 1 all custom 1000 tmp.dump id type mol x y z c_cluster c_spread
+dump_modify 1 sort id :pre
+
+In the first case, each snapshot in the tmp.phobic file will contain
+one line per molecule. Molecules with the same value are in the same
+micelle. In the second case each dump snapshot contains all atoms,
+each with a final field with the cluster ID of the micelle that the HP
+atoms of that atom's molecule belong to.
+
+The result from compute chunk/spread/atom can be used to define a new
+set of chunks, where all the atoms in all the molecules in the same
+micelle are assigned to the same chunk, i.e. one chunk per micelle.
+
+compute micelle all chunk/atom c_spread compress yes :pre
+
+Further analysis on a per-micelle basis can now be performed using any
+of the per-chunk computes listed on the "Howto chunk"_Howto_chunk.html
+doc page. E.g. count the number of atoms in each micelle, calculate
+its center or mass, shape (moments of intertia), radius of gyration,
+etc.
+
+compute prop all property/chunk micelle count
+fix 20 all ave/time 1000 1 1000 c_prop file tmp.micelle mode vector :pre
+
+Each snapshot in the tmp.micelle file will have one line per micelle
+with its count of atoms, plus a first line for a chunk with all the
+solvent atoms. By the time 50000 steps have elapsed there are a
+handful of large micelles.
+
+:line
+
+[Output info:]
+
+This compute calculates a global vector if a single input value is
+specified, otherwise a global array is output. The number of columns
+in the array is the number of inputs provided. The length of the
+vector or the number of vector elements or array rows = the number of
+chunks {Nchunk} as calculated by the specified "compute
+chunk/atom"_compute_chunk_atom.html command. The vector or array can
+be accessed by any command that uses global values from a compute as
+input. See the "Howto output"_Howto_output.html doc page for an
+overview of LAMMPS output options.
+
+The per-atom values for the vector or each column of the array will be
+in whatever "units"_units.html the corresponding input value is in.
+The vector or array values are "intensive".
+
+[Restrictions:] none
+
+[Related commands:]
+
+"compute chunk/atom"_compute_chunk_atom.html, "compute
+reduce"_compute_reduce.html, "compute
+chunk/spread/atom"_compute_chunk_spread_atom.html
+
+[Default:] none
diff --git a/doc/src/computes.txt b/doc/src/computes.txt
index 46dd30f7570f7a87d9874367f67f1d4af85a2fca..6528f78e40e38424828bc0349ec5689c8987d593 100644
--- a/doc/src/computes.txt
+++ b/doc/src/computes.txt
@@ -15,6 +15,7 @@ Computes :h1
compute_bond_local
compute_centro_atom
compute_chunk_atom
+ compute_chunk_spread_atom
compute_cluster_atom
compute_cna_atom
compute_cnp_atom
@@ -72,6 +73,7 @@ Computes :h1
compute_property_local
compute_rdf
compute_reduce
+ compute_reduce_chunk
compute_rigid_local
compute_saed
compute_slice
diff --git a/doc/src/lammps.book b/doc/src/lammps.book
index c296ff4039e0827464ba616081d6bb9630943056..c7f0b2fddd60ae395e084a1f5b2a735c2227fdff 100644
--- a/doc/src/lammps.book
+++ b/doc/src/lammps.book
@@ -406,6 +406,7 @@ compute_bond.html
compute_bond_local.html
compute_centro_atom.html
compute_chunk_atom.html
+compute_chunk_spread_atom.html
compute_cluster_atom.html
compute_cna_atom.html
compute_cnp_atom.html
@@ -463,6 +464,7 @@ compute_property_chunk.html
compute_property_local.html
compute_rdf.html
compute_reduce.html
+compute_reduce_chunk.html
compute_rigid_local.html
compute_saed.html
compute_slice.html
diff --git a/src/compute_chunk_spread_atom.cpp b/src/compute_chunk_spread_atom.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..93a4d30a2e313c305f1abc7875280e2ddaa7745e
--- /dev/null
+++ b/src/compute_chunk_spread_atom.cpp
@@ -0,0 +1,352 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include <cstring>
+#include <cstdlib>
+#include "compute_chunk_spread_atom.h"
+#include "atom.h"
+#include "update.h"
+#include "modify.h"
+#include "fix.h"
+#include "compute.h"
+#include "compute_chunk_atom.h"
+#include "input.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+enum{COMPUTE,FIX};
+
+#define INVOKED_VECTOR 2
+#define INVOKED_ARRAY 4
+#define INVOKED_PERATOM 8
+
+/* ---------------------------------------------------------------------- */
+
+ComputeChunkSpreadAtom::
+ComputeChunkSpreadAtom(LAMMPS *lmp, int narg, char **arg) :
+ Compute(lmp, narg, arg),
+ which(NULL), argindex(NULL), ids(NULL), value2index(NULL), idchunk(NULL)
+{
+ if (narg < 5) error->all(FLERR,"Illegal compute chunk/spread/atom command");
+
+ // ID of compute chunk/atom
+
+ int n = strlen(arg[3]) + 1;
+ idchunk = new char[n];
+ strcpy(idchunk,arg[3]);
+ init_chunk();
+
+ // expand args if any have wildcard character "*"
+
+ int iarg = 4;
+ int expand = 0;
+ char **earg;
+ int nargnew = input->expand_args(narg-iarg,&arg[iarg],1,earg);
+
+ if (earg != &arg[iarg]) expand = 1;
+ arg = earg;
+
+ // parse values
+
+ which = new int[nargnew];
+ argindex = new int[nargnew];
+ ids = new char*[nargnew];
+ value2index = new int[nargnew];
+ nvalues = 0;
+
+ iarg = 0;
+ while (iarg < nargnew) {
+ ids[nvalues] = NULL;
+
+ if (strncmp(arg[iarg],"c_",2) == 0 ||
+ strncmp(arg[iarg],"f_",2) == 0) {
+ if (arg[iarg][0] == 'c') which[nvalues] = COMPUTE;
+ else if (arg[iarg][0] == 'f') which[nvalues] = FIX;
+
+ int n = strlen(arg[iarg]);
+ char *suffix = new char[n];
+ strcpy(suffix,&arg[iarg][2]);
+
+ char *ptr = strchr(suffix,'[');
+ if (ptr) {
+ if (suffix[strlen(suffix)-1] != ']')
+ error->all(FLERR,"Illegal compute chunk/spread/atom command");
+ argindex[nvalues] = atoi(ptr+1);
+ *ptr = '\0';
+ } else argindex[nvalues] = 0;
+
+ n = strlen(suffix) + 1;
+ ids[nvalues] = new char[n];
+ strcpy(ids[nvalues],suffix);
+ nvalues++;
+ delete [] suffix;
+
+ } else error->all(FLERR,"Illegal compute chunk/spread/atom command");
+
+ iarg++;
+ }
+
+ // if wildcard expansion occurred, free earg memory from expand_args()
+
+ if (expand) {
+ for (int i = 0; i < nargnew; i++) delete [] earg[i];
+ memory->sfree(earg);
+ }
+
+ // setup and error check
+ // for compute, must calculate per-chunk values, i.e. style ends in "/chunk"
+ // for fix, assume a global vector or array is per-chunk data
+
+ for (int i = 0; i < nvalues; i++) {
+ if (which[i] == COMPUTE) {
+ int icompute = modify->find_compute(ids[i]);
+ if (icompute < 0)
+ error->all(FLERR,"Compute ID for compute chunk/spread/atom "
+ "does not exist");
+
+ char *ptr = strstr(modify->compute[icompute]->style,"/chunk");
+ if (!ptr || (ptr != modify->compute[icompute]->style +
+ strlen(modify->compute[icompute]->style) - strlen("/chunk")))
+ error->all(FLERR,"Compute for compute chunk/spread/atom "
+ "does not calculate per-chunk values");
+
+ if (argindex[i] == 0) {
+ if (!modify->compute[icompute]->vector_flag)
+ error->all(FLERR,"Compute chunk/spread/atom compute "
+ "does not calculate global vector");
+ } else {
+ if (!modify->compute[icompute]->array_flag)
+ error->all(FLERR,"Compute chunk/spread/atom compute "
+ "does not calculate global array");
+ if (argindex[i] > modify->compute[icompute]->size_array_cols)
+ error->all(FLERR,"Compute chunk/spread/atom compute array "
+ "is accessed out-of-range");
+ }
+
+ } else if (which[i] == FIX) {
+ int ifix = modify->find_fix(ids[i]);
+ if (ifix < 0)
+ error->all(FLERR,"Fix ID for compute chunk/spread/atom does not exist");
+ if (argindex[i] == 0) {
+ if (!modify->fix[ifix]->vector_flag)
+ error->all(FLERR,"Compute chunk/spread/atom fix "
+ "does not calculate global vector");
+ } else {
+ if (!modify->fix[ifix]->array_flag)
+ error->all(FLERR,"Compute chunk/spread/atom fix "
+ "does not calculate global array");
+ if (argindex[i] > modify->fix[ifix]->size_array_cols)
+ error->all(FLERR,"Compute chunk/spread/atom fix array "
+ "is accessed out-of-range");
+ }
+ }
+ }
+
+ // this compute produces a peratom vector or array
+
+ peratom_flag = 1;
+ if (nvalues == 1) size_peratom_cols = 0;
+ else size_peratom_cols = nvalues;
+
+ // per-atom vector or array
+
+ nmax = 0;
+ vector_atom = NULL;
+ array_atom = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeChunkSpreadAtom::~ComputeChunkSpreadAtom()
+{
+ delete [] idchunk;
+
+ delete [] which;
+ delete [] argindex;
+ for (int i = 0; i < nvalues; i++) delete [] ids[i];
+ delete [] ids;
+ delete [] value2index;
+
+ memory->destroy(vector_atom);
+ memory->destroy(array_atom);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeChunkSpreadAtom::init()
+{
+ init_chunk();
+
+ // set indices of all computes,fixes,variables
+
+ for (int m = 0; m < nvalues; m++) {
+ if (which[m] == COMPUTE) {
+ int icompute = modify->find_compute(ids[m]);
+ if (icompute < 0)
+ error->all(FLERR,"Compute ID for compute chunk/spread/atom "
+ "does not exist");
+ value2index[m] = icompute;
+
+ } else if (which[m] == FIX) {
+ int ifix = modify->find_fix(ids[m]);
+ if (ifix < 0)
+ error->all(FLERR,"Fix ID for compute chunk/spread/atom does not exist");
+ value2index[m] = ifix;
+ }
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeChunkSpreadAtom::init_chunk()
+{
+ int icompute = modify->find_compute(idchunk);
+ if (icompute < 0)
+ error->all(FLERR,"Chunk/atom compute does not exist for compute chunk/spread/atom");
+ cchunk = (ComputeChunkAtom *) modify->compute[icompute];
+ if (strcmp(cchunk->style,"chunk/atom") != 0)
+ error->all(FLERR,"Compute chunk/spread/atom does not use chunk/atom compute");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeChunkSpreadAtom::compute_peratom()
+{
+ invoked_peratom = update->ntimestep;
+
+ // grow local vector_atom or array_atom if necessary
+
+ if (atom->nmax > nmax) {
+ if (nvalues == 1) {
+ memory->destroy(vector_atom);
+ nmax = atom->nmax;
+ memory->create(vector_atom,nmax,"chunk/spread/atom:vector_atom");
+ } else {
+ memory->destroy(array_atom);
+ nmax = atom->nmax;
+ memory->create(array_atom,nmax,nvalues,"chunk/spread/atom:array_atom");
+ }
+ }
+
+ // compute chunk/atom assigns atoms to chunk IDs
+ // extract ichunk index vector from compute
+ // ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms
+
+ int nchunk = cchunk->setup_chunks();
+ cchunk->compute_ichunk();
+ int *ichunk = cchunk->ichunk;
+
+ // loop over values, access compute or fix
+ // loop over atoms, use chunk ID of each atom to store value from compute/fix
+
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ int i,m,n,index,nstride;
+ double *ptr;
+
+ for (m = 0; m < nvalues; m++) {
+ n = value2index[m];
+
+ // copy compute/fix values into vector_atom or array_atom
+ // nstride between values for each atom
+
+ if (nvalues == 1) {
+ ptr = vector_atom;
+ nstride = 1;
+ } else {
+ ptr = &array_atom[0][m];
+ nstride = nvalues;
+ }
+
+ // invoke compute if not previously invoked
+
+ if (which[m] == COMPUTE) {
+ Compute *compute = modify->compute[n];
+
+ if (argindex[m] == 0) {
+ if (!(compute->invoked_flag & INVOKED_VECTOR)) {
+ compute->compute_vector();
+ compute->invoked_flag |= INVOKED_VECTOR;
+ }
+ double *cvector = compute->vector;
+ for (i = 0; i < nlocal; i++, ptr += nstride) {
+ *ptr = 0.0;
+ if (!(mask[i] & groupbit)) continue;
+ index = ichunk[i]-1;
+ if (index < 0 || index >= nchunk) continue;
+ *ptr = cvector[index];
+ }
+
+ } else {
+ if (!(compute->invoked_flag & INVOKED_ARRAY)) {
+ compute->compute_array();
+ compute->invoked_flag |= INVOKED_ARRAY;
+ }
+ int icol = argindex[m]-1;
+ double **carray = compute->array;
+ for (i = 0; i < nlocal; i++, ptr += nstride) {
+ *ptr = 0.0;
+ if (!(mask[i] & groupbit)) continue;
+ index = ichunk[i]-1;
+ if (index < 0 || index >= nchunk) continue;
+ *ptr = carray[index][icol];
+ }
+ }
+
+ // access fix data, check if fix frequency is a match
+ // are assuming the fix global vector/array is per-chunk data
+ // check if index exceeds fix output length/rows
+
+ } else if (which[m] == FIX) {
+ Fix *fix = modify->fix[n];
+ if (update->ntimestep % fix->global_freq)
+ error->all(FLERR,"Fix used in compute chunk/spread/atom not "
+ "computed at compatible time");
+
+ if (argindex[m] == 0) {
+ int nfix = fix->size_vector;
+ for (i = 0; i < nlocal; i++, ptr += nstride) {
+ *ptr = 0.0;
+ if (!(mask[i] & groupbit)) continue;
+ index = ichunk[i]-1;
+ if (index < 0 || index >= nchunk || index >= nfix) continue;
+ *ptr = fix->compute_vector(index);
+ }
+
+ } else {
+ int icol = argindex[m]-1;
+ int nfix = fix->size_array_rows;
+ for (i = 0; i < nlocal; i++, ptr += nstride) {
+ *ptr = 0.0;
+ if (!(mask[i] & groupbit)) continue;
+ index = ichunk[i]-1;
+ if (index < 0 || index >= nchunk || index >= nfix) continue;
+ *ptr = fix->compute_array(index,icol);
+ }
+ }
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ memory usage of local atom-based array
+------------------------------------------------------------------------- */
+
+double ComputeChunkSpreadAtom::memory_usage()
+{
+ double bytes = nmax*nvalues * sizeof(double);
+ return bytes;
+}
diff --git a/src/compute_chunk_spread_atom.h b/src/compute_chunk_spread_atom.h
new file mode 100644
index 0000000000000000000000000000000000000000..9a4df080ca2d4de652f357f0a84c46f15fb260a7
--- /dev/null
+++ b/src/compute_chunk_spread_atom.h
@@ -0,0 +1,61 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef COMPUTE_CLASS
+
+ComputeStyle(chunk/spread/atom,ComputeChunkSpreadAtom)
+
+#else
+
+#ifndef LMP_COMPUTE_CHUNK_SPREAD_ATOM_H
+#define LMP_COMPUTE_CHUNK_SPREAD_ATOM_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeChunkSpreadAtom : public Compute {
+ public:
+ ComputeChunkSpreadAtom(class LAMMPS *, int, char **);
+ ~ComputeChunkSpreadAtom();
+ void init();
+ void compute_peratom();
+ double memory_usage();
+
+ protected:
+ int mode,nvalues;
+ char *idchunk;
+
+ int *which,*argindex,*value2index;
+ char **ids;
+
+ int nmax;
+ class ComputeChunkAtom *cchunk;
+
+ void init_chunk();
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory. Check the input script syntax and compare to the
+documentation for the command. You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+*/
diff --git a/src/compute_reduce.cpp b/src/compute_reduce.cpp
index 7c61d666b2d9db9c9625b69bb08dfdafff49fa29..e3c3c5b70a9c3adbad74163e5b0ac9149c24c813 100644
--- a/src/compute_reduce.cpp
+++ b/src/compute_reduce.cpp
@@ -11,6 +11,7 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
+#include <mpi.h>
#include <cstring>
#include <cstdlib>
#include "compute_reduce.h"
diff --git a/src/compute_reduce_chunk.cpp b/src/compute_reduce_chunk.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..40e9dad8dd4f9a7fcb2177551b2489066f06bc61
--- /dev/null
+++ b/src/compute_reduce_chunk.cpp
@@ -0,0 +1,512 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include <mpi.h>
+#include <string.h>
+#include "compute_reduce_chunk.h"
+#include "atom.h"
+#include "update.h"
+#include "modify.h"
+#include "fix.h"
+#include "compute.h"
+#include "compute_chunk_atom.h"
+#include "input.h"
+#include "variable.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+enum{SUM,MINN,MAXX};
+enum{COMPUTE,FIX,VARIABLE};
+
+#define INVOKED_PERATOM 8
+
+#define BIG 1.0e20
+
+/* ---------------------------------------------------------------------- */
+
+ComputeReduceChunk::ComputeReduceChunk(LAMMPS *lmp, int narg, char **arg) :
+ Compute(lmp, narg, arg),
+ vlocal(NULL), vglobal(NULL), alocal(NULL), aglobal(NULL), varatom(NULL),
+ which(NULL), argindex(NULL), value2index(NULL), idchunk(NULL), ids(NULL)
+{
+ if (narg < 6) error->all(FLERR,"Illegal compute reduce/chunk command");
+
+ // ID of compute chunk/atom
+
+ int n = strlen(arg[3]) + 1;
+ idchunk = new char[n];
+ strcpy(idchunk,arg[3]);
+ init_chunk();
+
+ // mode
+
+ if (strcmp(arg[4],"sum") == 0) mode = SUM;
+ else if (strcmp(arg[4],"min") == 0) mode = MINN;
+ else if (strcmp(arg[4],"max") == 0) mode = MAXX;
+ else error->all(FLERR,"Illegal compute reduce/chunk command");
+
+ int iarg = 5;
+
+ // expand args if any have wildcard character "*"
+
+ int expand = 0;
+ char **earg;
+ int nargnew = input->expand_args(narg-iarg,&arg[iarg],1,earg);
+
+ if (earg != &arg[iarg]) expand = 1;
+ arg = earg;
+
+ // parse values until one isn't recognized
+
+ which = new int[nargnew];
+ argindex = new int[nargnew];
+ ids = new char*[nargnew];
+ value2index = new int[nargnew];
+ nvalues = 0;
+
+ iarg = 0;
+ while (iarg < nargnew) {
+ ids[nvalues] = NULL;
+
+ if (strncmp(arg[iarg],"c_",2) == 0 ||
+ strncmp(arg[iarg],"f_",2) == 0 ||
+ strncmp(arg[iarg],"v_",2) == 0) {
+ if (arg[iarg][0] == 'c') which[nvalues] = COMPUTE;
+ else if (arg[iarg][0] == 'f') which[nvalues] = FIX;
+ else if (arg[iarg][0] == 'v') which[nvalues] = VARIABLE;
+
+ int n = strlen(arg[iarg]);
+ char *suffix = new char[n];
+ strcpy(suffix,&arg[iarg][2]);
+
+ char *ptr = strchr(suffix,'[');
+ if (ptr) {
+ if (suffix[strlen(suffix)-1] != ']')
+ error->all(FLERR,"Illegal compute reduce/chunk command");
+ argindex[nvalues] = atoi(ptr+1);
+ *ptr = '\0';
+ } else argindex[nvalues] = 0;
+
+ n = strlen(suffix) + 1;
+ ids[nvalues] = new char[n];
+ strcpy(ids[nvalues],suffix);
+ nvalues++;
+ delete [] suffix;
+
+ } else error->all(FLERR,"Illegal compute reduce/chunk command");
+
+ iarg++;
+ }
+
+ // if wildcard expansion occurred, free earg memory from expand_args()
+
+ if (expand) {
+ for (int i = 0; i < nargnew; i++) delete [] earg[i];
+ memory->sfree(earg);
+ }
+
+ // error check
+
+ for (int i = 0; i < nvalues; i++) {
+ if (which[i] == COMPUTE) {
+ int icompute = modify->find_compute(ids[i]);
+ if (icompute < 0)
+ error->all(FLERR,"Compute ID for compute reduce/chunk does not exist");
+ if (!modify->compute[icompute]->peratom_flag)
+ error->all(FLERR,"Compute reduce/chunk compute does not "
+ "calculate per-atom values");
+ if (argindex[i] == 0 &&
+ modify->compute[icompute]->size_peratom_cols != 0)
+ error->all(FLERR,"Compute reduce/chunk compute does not "
+ "calculate a per-atom vector");
+ if (argindex[i] && modify->compute[icompute]->size_peratom_cols == 0)
+ error->all(FLERR,"Compute reduce/chunk compute does not "
+ "calculate a per-atom array");
+ if (argindex[i] &&
+ argindex[i] > modify->compute[icompute]->size_peratom_cols)
+ error->all(FLERR,
+ "Compute reduce/chunk compute array is accessed out-of-range");
+
+ } else if (which[i] == FIX) {
+ int ifix = modify->find_fix(ids[i]);
+ if (ifix < 0)
+ error->all(FLERR,"Fix ID for compute reduce/chunk does not exist");
+ if (!modify->fix[ifix]->peratom_flag)
+ error->all(FLERR,"Compute reduce/chunk fix does not "
+ "calculate per-atom values");
+ if (argindex[i] == 0 &&
+ modify->fix[ifix]->size_peratom_cols != 0)
+ error->all(FLERR,"Compute reduce/chunk fix does not "
+ "calculate a per-atom vector");
+ if (argindex[i] && modify->fix[ifix]->size_peratom_cols == 0)
+ error->all(FLERR,"Compute reduce/chunk fix does not "
+ "calculate a per-atom array");
+ if (argindex[i] && argindex[i] > modify->fix[ifix]->size_peratom_cols)
+ error->all(FLERR,"Compute reduce/chunk fix array is "
+ "accessed out-of-range");
+
+ } else if (which[i] == VARIABLE) {
+ int ivariable = input->variable->find(ids[i]);
+ if (ivariable < 0)
+ error->all(FLERR,"Variable name for compute reduce/chunk does not exist");
+ if (input->variable->atomstyle(ivariable) == 0)
+ error->all(FLERR,"Compute reduce/chunk variable is "
+ "not atom-style variable");
+ }
+ }
+
+ // this compute produces either a vector or array
+
+ if (nvalues == 1) {
+ vector_flag = 1;
+ size_vector_variable = 1;
+ extvector = 0;
+ } else {
+ array_flag = 1;
+ size_array_rows_variable = 1;
+ size_array_cols = nvalues;
+ extarray = 0;
+ }
+
+ // setup
+
+ if (mode == SUM) initvalue = 0.0;
+ else if (mode == MINN) initvalue = BIG;
+ else if (mode == MAXX) initvalue = -BIG;
+
+ maxchunk = 0;
+ vlocal = vglobal = NULL;
+ alocal = aglobal = NULL;
+
+ maxatom = 0;
+ varatom = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeReduceChunk::~ComputeReduceChunk()
+{
+ delete [] idchunk;
+
+ delete [] which;
+ delete [] argindex;
+ for (int m = 0; m < nvalues; m++) delete [] ids[m];
+ delete [] ids;
+ delete [] value2index;
+
+ memory->destroy(vlocal);
+ memory->destroy(vglobal);
+ memory->destroy(alocal);
+ memory->destroy(aglobal);
+
+ memory->destroy(varatom);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeReduceChunk::init()
+{
+ init_chunk();
+
+ // set indices of all computes,fixes,variables
+
+ for (int m = 0; m < nvalues; m++) {
+ if (which[m] == COMPUTE) {
+ int icompute = modify->find_compute(ids[m]);
+ if (icompute < 0)
+ error->all(FLERR,"Compute ID for compute reduce/chunk does not exist");
+ value2index[m] = icompute;
+
+ } else if (which[m] == FIX) {
+ int ifix = modify->find_fix(ids[m]);
+ if (ifix < 0)
+ error->all(FLERR,"Fix ID for compute reduce/chunk does not exist");
+ value2index[m] = ifix;
+
+ } else if (which[m] == VARIABLE) {
+ int ivariable = input->variable->find(ids[m]);
+ if (ivariable < 0)
+ error->all(FLERR,"Variable name for compute reduce/chunk does not exist");
+ value2index[m] = ivariable;
+ }
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeReduceChunk::init_chunk()
+{
+ int icompute = modify->find_compute(idchunk);
+ if (icompute < 0)
+ error->all(FLERR,"Chunk/atom compute does not exist for "
+ "compute reduce/chunk");
+ cchunk = (ComputeChunkAtom *) modify->compute[icompute];
+ if (strcmp(cchunk->style,"chunk/atom") != 0)
+ error->all(FLERR,"Compute reduce/chunk does not use chunk/atom compute");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeReduceChunk::compute_vector()
+{
+ invoked_vector = update->ntimestep;
+
+ // compute chunk/atom assigns atoms to chunk IDs
+ // extract ichunk index vector from compute
+ // ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms
+
+ nchunk = cchunk->setup_chunks();
+ cchunk->compute_ichunk();
+ ichunk = cchunk->ichunk;
+ if (!nchunk) return;
+
+ size_vector = nchunk;
+
+ if (nchunk > maxchunk) {
+ memory->destroy(vlocal);
+ memory->destroy(vglobal);
+ maxchunk = nchunk;
+ memory->create(vlocal,maxchunk,"reduce/chunk:vlocal");
+ memory->create(vglobal,maxchunk,"reduce/chunk:vglobal");
+ vector = vglobal;
+ }
+
+ // perform local reduction of single peratom value
+
+ compute_one(0,vlocal,1);
+
+ // reduce the per-chunk values across all procs
+
+ if (mode == SUM)
+ MPI_Allreduce(vlocal,vglobal,nchunk,MPI_DOUBLE,MPI_SUM,world);
+ else if (mode == MINN)
+ MPI_Allreduce(vlocal,vglobal,nchunk,MPI_DOUBLE,MPI_MIN,world);
+ else if (mode == MAXX)
+ MPI_Allreduce(vlocal,vglobal,nchunk,MPI_DOUBLE,MPI_MAX,world);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeReduceChunk::compute_array()
+{
+ invoked_array = update->ntimestep;
+
+ // compute chunk/atom assigns atoms to chunk IDs
+ // extract ichunk index vector from compute
+ // ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms
+
+ nchunk = cchunk->setup_chunks();
+ cchunk->compute_ichunk();
+ ichunk = cchunk->ichunk;
+ if (!nchunk) return;
+
+ size_array_rows = nchunk;
+
+ if (nchunk > maxchunk) {
+ memory->destroy(alocal);
+ memory->destroy(aglobal);
+ maxchunk = nchunk;
+ memory->create(alocal,maxchunk,nvalues,"reduce/chunk:alocal");
+ memory->create(aglobal,maxchunk,nvalues,"reduce/chunk:aglobal");
+ array = aglobal;
+ }
+
+ // perform local reduction of all peratom values
+
+ for (int m = 0; m < nvalues; m++) compute_one(m,&alocal[0][m],nvalues);
+
+ // reduce the per-chunk values across all procs
+
+ if (mode == SUM)
+ MPI_Allreduce(&alocal[0][0],&aglobal[0][0],nchunk*nvalues,
+ MPI_DOUBLE,MPI_SUM,world);
+ else if (mode == MINN)
+ MPI_Allreduce(&alocal[0][0],&aglobal[0][0],nchunk*nvalues,
+ MPI_DOUBLE,MPI_MIN,world);
+ else if (mode == MAXX)
+ MPI_Allreduce(&alocal[0][0],&aglobal[0][0],nchunk*nvalues,
+ MPI_DOUBLE,MPI_MAX,world);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeReduceChunk::compute_one(int m, double *vchunk, int nstride)
+{
+ // initialize per-chunk values in accumulation vector
+
+ for (int i = 0; i < nchunk; i += nstride) vchunk[i] = initvalue;
+
+ // loop over my atoms
+ // use peratom input and chunk ID of each atom to update vector
+
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ int index;
+
+ if (which[m] == COMPUTE) {
+ Compute *compute = modify->compute[value2index[m]];
+
+ if (!(compute->invoked_flag & INVOKED_PERATOM)) {
+ compute->compute_peratom();
+ compute->invoked_flag |= INVOKED_PERATOM;
+ }
+
+ if (argindex[m] == 0) {
+ double *vcompute = compute->vector_atom;
+ for (int i = 0; i < nlocal; i++) {
+ if (!(mask[i] & groupbit)) continue;
+ index = ichunk[i]-1;
+ if (index < 0) continue;
+ combine(vchunk[index*nstride],vcompute[i]);
+ }
+ } else {
+ double **acompute = compute->array_atom;
+ int argindexm1 = argindex[m] - 1;
+ for (int i = 0; i < nlocal; i++) {
+ if (!(mask[i] & groupbit)) continue;
+ index = ichunk[i]-1;
+ if (index < 0) continue;
+ combine(vchunk[index*nstride],acompute[i][argindexm1]);
+ }
+ }
+
+ // access fix fields, check if fix frequency is a match
+
+ } else if (which[m] == FIX) {
+ Fix *fix = modify->fix[value2index[m]];
+ if (update->ntimestep % fix->peratom_freq)
+ error->all(FLERR,"Fix used in compute reduce/chunk not "
+ "computed at compatible time");
+
+ if (argindex[m] == 0) {
+ double *vfix = fix->vector_atom;
+ for (int i = 0; i < nlocal; i++) {
+ if (!(mask[i] & groupbit)) continue;
+ index = ichunk[i]-1;
+ if (index < 0) continue;
+ combine(vchunk[index*nstride],vfix[i]);
+ }
+ } else {
+ double **afix = fix->array_atom;
+ int argindexm1 = argindex[m] - 1;
+ for (int i = 0; i < nlocal; i++) {
+ if (!(mask[i] & groupbit)) continue;
+ index = ichunk[i]-1;
+ if (index < 0) continue;
+ combine(vchunk[index*nstride],afix[i][argindexm1]);
+ }
+ }
+
+ // evaluate atom-style variable
+
+ } else if (which[m] == VARIABLE) {
+ if (atom->nmax > maxatom) {
+ memory->destroy(varatom);
+ maxatom = atom->nmax;
+ memory->create(varatom,maxatom,"reduce/chunk:varatom");
+ }
+
+ input->variable->compute_atom(value2index[m],igroup,varatom,1,0);
+ for (int i = 0; i < nlocal; i++) {
+ if (!(mask[i] & groupbit)) continue;
+ index = ichunk[i]-1;
+ if (index < 0) continue;
+ combine(vchunk[index*nstride],varatom[i]);
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ combine two values according to reduction mode
+------------------------------------------------------------------------- */
+
+void ComputeReduceChunk::combine(double &one, double two)
+{
+ if (mode == SUM) one += two;
+ else if (mode == MINN) {
+ if (two < one) one = two;
+ } else if (mode == MAXX) {
+ if (two > one) one = two;
+ }
+}
+
+/* ----------------------------------------------------------------------
+ lock methods: called by fix ave/time
+ these methods insure vector/array size is locked for Nfreq epoch
+ by passing lock info along to compute chunk/atom
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ increment lock counter
+------------------------------------------------------------------------- */
+
+void ComputeReduceChunk::lock_enable()
+{
+ cchunk->lockcount++;
+}
+
+/* ----------------------------------------------------------------------
+ decrement lock counter in compute chunk/atom, if it still exists
+------------------------------------------------------------------------- */
+
+void ComputeReduceChunk::lock_disable()
+{
+ int icompute = modify->find_compute(idchunk);
+ if (icompute >= 0) {
+ cchunk = (ComputeChunkAtom *) modify->compute[icompute];
+ cchunk->lockcount--;
+ }
+}
+
+/* ----------------------------------------------------------------------
+ calculate and return # of chunks = length of vector/array
+------------------------------------------------------------------------- */
+
+int ComputeReduceChunk::lock_length()
+{
+ nchunk = cchunk->setup_chunks();
+ return nchunk;
+}
+
+/* ----------------------------------------------------------------------
+ set the lock from startstep to stopstep
+------------------------------------------------------------------------- */
+
+void ComputeReduceChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep)
+{
+ cchunk->lock(fixptr,startstep,stopstep);
+}
+
+/* ----------------------------------------------------------------------
+ unset the lock
+------------------------------------------------------------------------- */
+
+void ComputeReduceChunk::unlock(Fix *fixptr)
+{
+ cchunk->unlock(fixptr);
+}
+
+/* ----------------------------------------------------------------------
+ memory usage of local data
+------------------------------------------------------------------------- */
+
+double ComputeReduceChunk::memory_usage()
+{
+ double bytes = (bigint) maxatom * sizeof(double);
+ if (nvalues == 1) bytes += (bigint) maxchunk * 2 * sizeof(double);
+ else bytes += (bigint) maxchunk * nvalues * 2 * sizeof(double);
+ return bytes;
+}
diff --git a/src/compute_reduce_chunk.h b/src/compute_reduce_chunk.h
new file mode 100644
index 0000000000000000000000000000000000000000..ac556907b7352cbcf7fb47d44e13a537286bdbae
--- /dev/null
+++ b/src/compute_reduce_chunk.h
@@ -0,0 +1,77 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef COMPUTE_CLASS
+
+ComputeStyle(reduce/chunk,ComputeReduceChunk)
+
+#else
+
+#ifndef LMP_COMPUTE_REDUCE_CHUNK_H
+#define LMP_COMPUTE_REDUCE_CHUNK_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeReduceChunk : public Compute {
+ public:
+ ComputeReduceChunk(class LAMMPS *, int, char **);
+ ~ComputeReduceChunk();
+ void init();
+ void compute_vector();
+ void compute_array();
+
+ void lock_enable();
+ void lock_disable();
+ int lock_length();
+ void lock(class Fix *, bigint, bigint);
+ void unlock(class Fix *);
+
+ double memory_usage();
+
+ private:
+ int mode,nvalues;
+ int *which,*argindex,*value2index;
+ char *idchunk;
+ char **ids;
+
+ int nchunk;
+ int maxchunk,maxatom;
+ double initvalue;
+ double *vlocal,*vglobal;
+ double **alocal,**aglobal;
+ double *varatom;
+
+ class ComputeChunkAtom *cchunk;
+ int *ichunk;
+
+ void init_chunk();
+ void compute_one(int, double *, int);
+ void combine(double &, double);
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory. Check the input script syntax and compare to the
+documentation for the command. You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+*/