2 new computes: chunk/spread/atom and reduce/chunk

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,

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.

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

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

doc/src/compute_chunk_spread_atom.txt

+"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

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
 

doc/src/compute_reduce_chunk.txt

+"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

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

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

src/compute_chunk_spread_atom.cpp

+/* ----------------------------------------------------------------------
+   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;
+}

src/compute_chunk_spread_atom.h

+/* -*- 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.
+
+*/

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"

src/compute_reduce_chunk.cpp

+/* ----------------------------------------------------------------------
+   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;
+}

src/compute_reduce_chunk.h

+/* -*- 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.
+
+*/