diff --git a/doc/fix_ave_spatial.html b/doc/fix_ave_spatial.html
index 27f1c3fac..1653aa853 100644
--- a/doc/fix_ave_spatial.html
+++ b/doc/fix_ave_spatial.html
@@ -1,332 +1,336 @@
 <HTML>
 <CENTER><A HREF = "http://lammps.sandia.gov">LAMMPS WWW Site</A> - <A HREF = "Manual.html">LAMMPS Documentation</A> - <A HREF = "Section_commands.html#comm">LAMMPS Commands</A> 
 </CENTER>
 
 
 
 
 
 
 <HR>
 
 <H3>fix ave/spatial command 
 </H3>
 <P><B>Syntax:</B>
 </P>
 <PRE>fix ID group-ID ave/spatial Nevery Nrepeat Nfreq dim origin delta value1 value2 ... keyword args ... 
 </PRE>
 <UL><LI>ID, group-ID are documented in <A HREF = "fix.html">fix</A> command 
 
 <LI>ave/spatial = style name of this fix command 
 
 <LI>Nevery = use input values every this many timesteps 
 
 <LI>Nrepeat = # of times to use input values for calculating averages 
 
 <LI>Nfreq = calculate averages every this many timesteps 
 
 <LI>dim = <I>x</I> or <I>y</I> or <I>z</I> 
 
 <LI>origin = <I>lower</I> or <I>center</I> or <I>upper</I> or coordinate value (distance units) 
 
 <LI>delta = thickness of spatial layers in dim (distance units) 
 
 <LI>one or more input values can be listed 
 
 <LI>value = x, y, z, vx, vy, vz, fx, fy, fz, density/mass, density/number, c_ID, c_ID[I], f_ID, f_ID[I], v_name 
 
 <PRE>  x,y,z,vx,vy,vz,fx,fy,fz = atom attribute (velocity, force component)
   density/number, density/mass = number or mass density
   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
   f_ID[I] = Ith column of per-atom array calculated by a fix with ID
   v_name = per-atom vector calculated by an atom-style variable with name 
 </PRE>
 <LI>zero or more keyword/arg pairs may be appended 
 
 <LI>keyword = <I>norm</I> or <I>units</I> or <I>file</I> or <I>ave</I> or <I>title1</I> or <I>title2</I> or <I>title3</I> 
 
 <PRE>  <I>units</I> arg = <I>box</I> or <I>lattice</I> or <I>reduced</I>
   <I>norm</I> arg = <I>all</I> or <I>sample</I>
+  <I>region</I> arg = region-ID
+    region-ID = ID of region atoms must be in to contribute to spatial averaging
   <I>ave</I> args = <I>one</I> or <I>running</I> or <I>window M</I>
     one = output new average value every Nfreq steps
     running = output cumulative average of all previous Nfreq steps
     window M = output average of M most recent Nfreq steps
   <I>file</I> arg = filename
     filename = file to write results to
   <I>title1</I> arg = string
     string = text to print as 1st line of output file
   <I>title2</I> arg = string
     string = text to print as 2nd line of output file
   <I>title3</I> arg = string
     string = text to print as 3rd line of output file 
 </PRE>
 
 </UL>
 <P><B>Examples:</B>
 </P>
 <PRE>fix 1 all ave/spatial 10000 1 10000 z lower 0.02 c_myCentro units reduced &
                       title1 "My output values"
 fix 1 flow ave/spatial 100 10 1000 y 0.0 1.0 vx vz norm sample file vel.profile
 fix 1 flow ave/spatial 100 5 1000 y 0.0 2.5 density/mass ave running 
 </PRE>
 <P><B>Description:</B>
 </P>
 <P>Use one or more per-atom vectors as inputs every few timesteps, bin
 them spatially by layer in a dimension, and average the layer values
 over longer timescales.  The resulting layer averages can be used by
 other <A HREF = "Section_howto.html#4_15">output commands</A> such as <A HREF = "thermo_style.html">thermo_style
 custom</A>, and can also be written to a file.
 </P>
 <P>The group specified with the command means only atoms within the group
-contribute to layer averages.
+contribute to layer averages.  If the <I>region</I> keyword is used, the
+atom must be in both the group and the specified geometric
+<A HREF = "region.html">region</A> in order to contribute to layer averages.
 </P>
 <P>Each listed value can be an atom attribute (position, velocity, force
 component), a mass or number density, or the result of a
 <A HREF = "compute.html">compute</A> or <A HREF = "fix.html">fix</A> or the evaluation of an
 atom-style <A HREF = "variable.html">variable</A>.  In the latter cases, the
 compute, fix, or variable must produce a per-atom quantity, not a
 global quantity.  If you wish to time-average global quantities from a
 compute, fix, or variable, then see the <A HREF = "fix_ave_time.html">fix
 ave/time</A> command.
 </P>
 <P><A HREF = "compute.html">Computes</A> that produce per-atom quantities are those
 which have the word <I>atom</I> in their style name.  See the doc pages for
 individual <A HREF = "fix.html">fixes</A> to determine which ones produce per-atom
 quantities.  <A HREF = "variable.html">Variables</A> of style <I>atom</I> are the only
 ones that can be used with this fix since all other styles of variable
 produce global quantities.
 </P>
 <P>The per-atom values of each input vector are binned and averaged
 independently of the per-atom values in other input vectors.
 </P>
 <HR>
 
 <P>The <I>Nevery</I>, <I>Nrepeat</I>, and <I>Nfreq</I> arguments specify on what
 timesteps the input values will be used to bin them into layers and
 contribute to the average.  The final averaged quantities are
 generated every <I>Nfreq</I> timesteps.  The average is over <I>Nrepeat</I>
 quantities, computed in the preceding portion of the simulation every
 <I>Nevery</I> timesteps.  <I>Nfreq</I> must be a multiple of <I>Nevery</I> and
 <I>Nevery</I> must be non-zero even if <I>Nrepeat</I> is 1.  Also, the timesteps
 contributing to the average value cannot overlap, i.e. Nfreq >
 (Nrepeat-1)*Nevery is required.
 </P>
 <P>For example, if Nevery=2, Nrepeat=6, and Nfreq=100, then values on
 timesteps 90,92,94,96,98,100 will be used to compute the final average
 on timestep 100.  Similarly for timesteps 190,192,194,196,198,200 on
 timestep 200, etc.  If Nrepeat=1 and Nfreq = 100, then no time
 averaging is done; values are simply generated on timesteps
 100,200,etc.
 </P>
 <HR>
 
 <P>Each per-atom property is also averaged over atoms in each layer,
 where the layers are in a particular <I>dim</I> and have a thickness given
 by <I>delta</I>.  Every Nfreq steps, when an averaging is being performed
 and the per-atom property is calculated for the first time, the number
 of layers and the layer boundaries are computed.  Thus if the
 simulation box changes size during a simulation, the number of layers
 and their boundaries may also change.  Layers are defined relative to
 a specified <I>origin</I>, which may be the lower/upper edge of the box (in
 <I>dim</I>) or its center point, or a specified coordinate value.  Starting
 at the origin, sufficient layers are created in both directions to
 completely cover the box.  On subsequent timesteps every atom is
 mapped to one of the layers.  Atoms beyond the lowermost/uppermost
 layer are counted in the first/last layer.
 </P>
 <P>For orthogonal simulation boxes, the layers are "slices" aligned with
 the xyz coordinate axes.  For non-orthogonal (triclinic) simulation
 boxes, the layers are "tilted slices" that are parallel to the tilted
 faces of the box.  See the <A HREF = "region.html">region prism</A> command for a
 discussion of the geometry of tilted boxes in LAMMPS.  As described
 there, a tilted simulation box has edge vectors a,b,c.  In that
 nomenclature, layers in the x dimension have faces with normals in the
 "b" cross "c" direction.  Layers in y have faces normal to the "a"
 cross "c" direction.  And layers in z have faces normal to the "a"
 cross "b" direction.  Note that in order to define the thickness and
 position of these tilted layers in an unambiguous fashion, the <I>units</I>
 option must be set to <I>reduced</I> when using a non-orthogonal simulation
 box, as discussed below.
 </P>
 <HR>
 
 <P>The atom attribute values (x,y,z,vx,vy,vz,fx,fy,fz) are
 self-explanatory.  Note that other atom attributes can be used as
 inputs to this fix by using the <A HREF = "compute_property_atom.html">compute
 property/atom</A> command and then specifying
 an input value from that compute.
 </P>
 <P>The <I>density/number</I> value means the number density is computed in
 each layer, i.e. a weighting of 1 for each atom.  The <I>density/mass</I>
 value means the mass density is computed in each layer, i.e. each atom
 is weighted by its mass.  The resulting density is normalized by the
 volume of the layer so that units of number/volume or density are
 output.  See the <A HREF = "units.html">units</A> command doc page for the
 definition of density for each choice of units, e.g. gram/cm^3.
 </P>
 <P>If a value begins with "c_", a compute ID must follow which has been
 previously defined in the input script.  If no bracketed integer is
 appended, the per-atom vector calculated by the compute is used.  If a
 bracketed interger is appended, the Ith column of the per-atom array
 calculated by the compute is used.  Users can also write code for
 their own compute styles and <A HREF = "Section_modify.html">add them to LAMMPS</A>.
 </P>
 <P>If a value begins with "f_", a fix ID must follow which has been
 previously defined in the input script.  If no bracketed integer is
 appended, the per-atom vector calculated by the fix is used.  If a
 bracketed integer is appended, the Ith column of the per-atom array
 calculated by the fix is used.  Note that some fixes only produce
 their values on certain timesteps, which must be compatible with
 <I>Nevery</I>, else an error results.  Users can also write code for their
 own fix styles and <A HREF = "Section_modify.html">add them to LAMMPS</A>.
 </P>
 <P>If a value begins with "v_", a variable name must follow which has
 been previously defined in the input script.  Variables of style
 <I>atom</I> can reference thermodynamic keywords and various per-atom
 attributes, or invoke other computes, fixes, or variables when they
 are evaluated, so this is a very general means of generating per-atom
 quantities to spatially average.
 </P>
 <HR>
 
 <P>Additional optional keywords also affect the operation of this fix.
 </P>
 <P>The <I>units</I> keyword determines the meaning of the distance units used
 for the layer thickness <I>delta</I> and for <I>origin</I> if it is a coordinate
 value.  For orthogonal simulation boxes, any of the 3 options may be
 used.  For non-orthogonal (triclinic) simulation boxes, only the
 <I>reduced</I> option may be used.
 </P>
 <P>A <I>box</I> value selects standard distance units as defined by the
 <A HREF = "units.html">units</A> command, e.g. Angstroms for units = real or metal.
 A <I>lattice</I> value means the distance units are in lattice spacings.
 The <A HREF = "lattice.html">lattice</A> command must have been previously used to
 define the lattice spacing.  A <I>reduced</I> value means normalized
 unitless values between 0 and 1, which represent the lower and upper
 faces of the simulation box respectively.  Thus an <I>origin</I> value of
 0.5 means the center of the box in any dimension.  A <I>delta</I> value of
 0.1 means 10 layers span the box in any dimension. 
 </P>
 <P>Consider a non-orthogonal box, with layers in the x dimension.  No
 matter how the box is tilted, an <I>origin</I> of 0.0 means start layers at
 the lower "b" cross "c" plane of the simulation box and an <I>origin</I> of
 1.0 means to start layers at the upper "b" cross "c" face of the box.
 A <I>delta</I> value of 0.1 means there will be 10 layers from 0.0 to 1.0,
 regardless of the current size or shape of the simulation box.
 </P>
 <P>The <I>norm</I> keyword affects how averaging is done for the output
 produced every <I>Nfreq</I> timesteps.  For an <I>all</I> setting, a layer
 quantity is summed over all atoms in all <I>Nrepeat</I> samples, as is the
 count of atoms in the layer.  The printed value for the layer is
 Total-quantity / Total-count.  In other words it is an average over
 the entire <I>Nfreq</I> timescale.
 </P>
 <P>For a <I>sample</I> setting, the layer quantity is summed over atoms for
 only a single sample, as is the count, and a "average sample value" is
 computed, i.e. Sample-quantity / Sample-count.  The printed value for
 the layer is the average of the <I>Nrepeat</I> "average sample values", In
 other words it is an average of an average.
 </P>
 <P>The <I>ave</I> keyword determines how the layer values produced every
 <I>Nfreq</I> steps are averaged with layer values produced on previous
 steps that were multiples of <I>Nfreq</I>, before they are accessed by
 another output command or written to a file.
 </P>
 <P>If the <I>ave</I> setting is <I>one</I>, then the layer values produced on
 timesteps that are multiples of <I>Nfreq</I> are independent of each other;
 they are output as-is without further averaging.
 </P>
 <P>If the <I>ave</I> setting is <I>running</I>, then the layer values produced on
 timesteps that are multiples of <I>Nfreq</I> are summed and averaged in a
 cumulative sense before being output.  Each output layer value is
 thus the average of the layer value produced on that timestep with all
 preceding values for the same layer.  This running average begins
 when the fix is defined; it can only be restarted by deleting the fix
 via the <A HREF = "unfix.html">unfix</A> command, or re-defining the fix by
 re-specifying it.
 </P>
 <P>If the <I>ave</I> setting is <I>window</I>, then the layer values produced on
 timesteps that are multiples of <I>Nfreq</I> are summed and averaged within
 a moving "window" of time, so that the last M values for the same
 layer are used to produce the output.  E.g. if M = 3 and Nfreq = 1000,
 then the output on step 10000 will be the average of the individual
 layer values on steps 8000,9000,10000.  Outputs on early steps will
 average over less than M values if they are not available.
 </P>
 <P>The <I>file</I> keyword allows a filename to be specified.  Every <I>Nfreq</I>
 timesteps, a section of layer info will be written to a text file in
 the following format.  A line with the timestep and number of layers
 is written.  Then one line per layer is written, containing the layer
 ID (1-N), the coordinate of the center of the layer, the number of
 atoms in the layer, and one or more calculated values.  The number of
 values in each line corresponds to the number of values specified in
 the fix ave/spatial command.  The number of atoms and the value(s) are
 average quantities.  If the value of the <I>units</I> keyword is <I>box</I> or
 <I>lattice</I>, the "coord" is printed in box units.  If the value of the
 <I>units</I> keyword is <I>reduced</I>, the "coord" is printed in reduced units
 (0-1).
 </P>
 <P>The <I>title1</I> and <I>title2</I> and <I>title3</I> keywords allow specification of
 the strings that will be printed as the first 3 lines of the output
 file, assuming the <I>file</I> keyword was used.  LAMMPS uses default
 values for each of these, so they do not need to be specified.
 </P>
 <P>By default, these header lines are as follows:
 </P>
 <PRE># Spatial-averaged data for fix ID and group name
 # Timestep Number-of-layers
 # Layer Coord Count value1 value2 ... 
 </PRE>
 <P>In the first line, ID and name are replaced with the fix-ID and group
 name.  The second line describes the two values that are printed at
 the first of each section of output.  In the third line the values are
 replaced with the appropriate fields from the fix ave/spatial command.
 </P>
 <HR>
 
 <P><B>Restart, fix_modify, output, run start/stop, minimize info:</B>
 </P>
 <P>No information about this fix is written to <A HREF = "restart.html">binary restart
 files</A>.  None of the <A HREF = "fix_modify.html">fix_modify</A> options
 are relevant to this fix.
 </P>
 <P>This fix computes a global array of values which can be accessed by
 various <A HREF = "Section_howto.html#4_15">output commands</A>.  The values can
 only be accessed on timesteps that are multiples of <I>Nfreq</I> since that
 is when averaging is performed.  The global array has # of rows =
 Nlayers and # of columns = Nvalues+2.  The first column has the layer
 coordinate (center of the layer), the 2nd column has the count of
 atoms in that layer, and the remaining columns are the Nvalue
 quantities.  When the array is accessed with an I that exceeds the
 current number of layers, than a 0.0 is returned by the fix instead of
 an error, since the number of layers can vary as a simulation runs,
 depending on the simulation box size.  The array values calculated by
 this fix are "intensive", since they are already normalized by the
 count of atoms in each layer.
 </P>
 <P>No parameter of this fix can be used with the <I>start/stop</I> keywords of
 the <A HREF = "run.html">run</A> command.  This fix is not invoked during <A HREF = "minimize.html">energy
 minimization</A>.
 </P>
 <P><B>Restrictions:</B>
 </P>
 <P>When the <I>ave</I> keyword is set to <I>running</I> or <I>window</I> then the number
 of layers must remain the same during the simulation, so that the
 appropriate averaging can be done.  This will be the case if the
 simulation box size doesn't change or if the <I>units</I> keyword is set to
 <I>reduced</I>.
 </P>
 <P><B>Related commands:</B>
 </P>
 <P><A HREF = "compute.html">compute</A>, <A HREF = "fix_ave_atom.html">fix ave/atom</A>, <A HREF = "fix_ave_histo.html">fix
 ave/histo</A>, <A HREF = "fix_ave_time.html">fix ave/time</A>,
 <A HREF = "variable.html">variable</A>, <A HREF = "fix_ave_correlate.html">fix ave/correlate</A>,
 </P>
 <P><B>Default:</B>
 </P>
 <P>The option defaults are units = lattice, norm = all, no file output,
 and ave = one, title 1,2,3 = strings as described above.
 </P>
 </HTML>
diff --git a/doc/fix_ave_spatial.txt b/doc/fix_ave_spatial.txt
index 9e6611a96..11880661b 100644
--- a/doc/fix_ave_spatial.txt
+++ b/doc/fix_ave_spatial.txt
@@ -1,314 +1,318 @@
 "LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
 
 :link(lws,http://lammps.sandia.gov)
 :link(ld,Manual.html)
 :link(lc,Section_commands.html#comm)
 
 :line
 
 fix ave/spatial command :h3
 
 [Syntax:]
 
 fix ID group-ID ave/spatial Nevery Nrepeat Nfreq dim origin delta value1 value2 ... keyword args ... :pre
 
 ID, group-ID are documented in "fix"_fix.html command :ulb,l
 ave/spatial = style name of this fix command :l
 Nevery = use input values every this many timesteps :l
 Nrepeat = # of times to use input values for calculating averages :l
 Nfreq = calculate averages every this many timesteps :l
 dim = {x} or {y} or {z} :l
 origin = {lower} or {center} or {upper} or coordinate value (distance units) :l
 delta = thickness of spatial layers in dim (distance units) :l
 one or more input values can be listed :l
 value = x, y, z, vx, vy, vz, fx, fy, fz, density/mass, density/number, c_ID, c_ID\[I\], f_ID, f_ID\[I\], v_name :l
   x,y,z,vx,vy,vz,fx,fy,fz = atom attribute (velocity, force component)
   density/number, density/mass = number or mass density
   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
   f_ID\[I\] = Ith column of per-atom array calculated by a fix with ID
   v_name = per-atom vector calculated by an atom-style variable with name :pre
 
 zero or more keyword/arg pairs may be appended :l
 keyword = {norm} or {units} or {file} or {ave} or {title1} or {title2} or {title3} :l
   {units} arg = {box} or {lattice} or {reduced}
   {norm} arg = {all} or {sample}
+  {region} arg = region-ID
+    region-ID = ID of region atoms must be in to contribute to spatial averaging
   {ave} args = {one} or {running} or {window M}
     one = output new average value every Nfreq steps
     running = output cumulative average of all previous Nfreq steps
     window M = output average of M most recent Nfreq steps
   {file} arg = filename
     filename = file to write results to
   {title1} arg = string
     string = text to print as 1st line of output file
   {title2} arg = string
     string = text to print as 2nd line of output file
   {title3} arg = string
     string = text to print as 3rd line of output file :pre
 :ule
 
 [Examples:]
 
 fix 1 all ave/spatial 10000 1 10000 z lower 0.02 c_myCentro units reduced &
                       title1 "My output values"
 fix 1 flow ave/spatial 100 10 1000 y 0.0 1.0 vx vz norm sample file vel.profile
 fix 1 flow ave/spatial 100 5 1000 y 0.0 2.5 density/mass ave running :pre
 
 [Description:]
 
 Use one or more per-atom vectors as inputs every few timesteps, bin
 them spatially by layer in a dimension, and average the layer values
 over longer timescales.  The resulting layer averages can be used by
 other "output commands"_Section_howto.html#4_15 such as "thermo_style
 custom"_thermo_style.html, and can also be written to a file.
 
 The group specified with the command means only atoms within the group
-contribute to layer averages.
+contribute to layer averages.  If the {region} keyword is used, the
+atom must be in both the group and the specified geometric
+"region"_region.html in order to contribute to layer averages.
 
 Each listed value can be an atom attribute (position, velocity, force
 component), a mass or number density, or the result of a
 "compute"_compute.html or "fix"_fix.html or the evaluation of an
 atom-style "variable"_variable.html.  In the latter cases, the
 compute, fix, or variable must produce a per-atom quantity, not a
 global quantity.  If you wish to time-average global quantities from a
 compute, fix, or variable, then see the "fix
 ave/time"_fix_ave_time.html command.
 
 "Computes"_compute.html that produce per-atom quantities are those
 which have the word {atom} in their style name.  See the doc pages for
 individual "fixes"_fix.html to determine which ones produce per-atom
 quantities.  "Variables"_variable.html of style {atom} are the only
 ones that can be used with this fix since all other styles of variable
 produce global quantities.
 
 The per-atom values of each input vector are binned and averaged
 independently of the per-atom values in other input vectors.
 
 :line
 
 The {Nevery}, {Nrepeat}, and {Nfreq} arguments specify on what
 timesteps the input values will be used to bin them into layers and
 contribute to the average.  The final averaged quantities are
 generated every {Nfreq} timesteps.  The average is over {Nrepeat}
 quantities, computed in the preceding portion of the simulation every
 {Nevery} timesteps.  {Nfreq} must be a multiple of {Nevery} and
 {Nevery} must be non-zero even if {Nrepeat} is 1.  Also, the timesteps
 contributing to the average value cannot overlap, i.e. Nfreq >
 (Nrepeat-1)*Nevery is required.
 
 For example, if Nevery=2, Nrepeat=6, and Nfreq=100, then values on
 timesteps 90,92,94,96,98,100 will be used to compute the final average
 on timestep 100.  Similarly for timesteps 190,192,194,196,198,200 on
 timestep 200, etc.  If Nrepeat=1 and Nfreq = 100, then no time
 averaging is done; values are simply generated on timesteps
 100,200,etc.
 
 :line
 
 Each per-atom property is also averaged over atoms in each layer,
 where the layers are in a particular {dim} and have a thickness given
 by {delta}.  Every Nfreq steps, when an averaging is being performed
 and the per-atom property is calculated for the first time, the number
 of layers and the layer boundaries are computed.  Thus if the
 simulation box changes size during a simulation, the number of layers
 and their boundaries may also change.  Layers are defined relative to
 a specified {origin}, which may be the lower/upper edge of the box (in
 {dim}) or its center point, or a specified coordinate value.  Starting
 at the origin, sufficient layers are created in both directions to
 completely cover the box.  On subsequent timesteps every atom is
 mapped to one of the layers.  Atoms beyond the lowermost/uppermost
 layer are counted in the first/last layer.
 
 For orthogonal simulation boxes, the layers are "slices" aligned with
 the xyz coordinate axes.  For non-orthogonal (triclinic) simulation
 boxes, the layers are "tilted slices" that are parallel to the tilted
 faces of the box.  See the "region prism"_region.html command for a
 discussion of the geometry of tilted boxes in LAMMPS.  As described
 there, a tilted simulation box has edge vectors a,b,c.  In that
 nomenclature, layers in the x dimension have faces with normals in the
 "b" cross "c" direction.  Layers in y have faces normal to the "a"
 cross "c" direction.  And layers in z have faces normal to the "a"
 cross "b" direction.  Note that in order to define the thickness and
 position of these tilted layers in an unambiguous fashion, the {units}
 option must be set to {reduced} when using a non-orthogonal simulation
 box, as discussed below.
 
 :line
 
 The atom attribute values (x,y,z,vx,vy,vz,fx,fy,fz) are
 self-explanatory.  Note that other atom attributes can be used as
 inputs to this fix by using the "compute
 property/atom"_compute_property_atom.html command and then specifying
 an input value from that compute.
 
 The {density/number} value means the number density is computed in
 each layer, i.e. a weighting of 1 for each atom.  The {density/mass}
 value means the mass density is computed in each layer, i.e. each atom
 is weighted by its mass.  The resulting density is normalized by the
 volume of the layer so that units of number/volume or density are
 output.  See the "units"_units.html command doc page for the
 definition of density for each choice of units, e.g. gram/cm^3.
 
 If a value begins with "c_", a compute ID must follow which has been
 previously defined in the input script.  If no bracketed integer is
 appended, the per-atom vector calculated by the compute is used.  If a
 bracketed interger is appended, the Ith column of the per-atom array
 calculated by the compute is used.  Users can also write code for
 their own compute styles and "add them to LAMMPS"_Section_modify.html.
 
 If a value begins with "f_", a fix ID must follow which has been
 previously defined in the input script.  If no bracketed integer is
 appended, the per-atom vector calculated by the fix is used.  If a
 bracketed integer is appended, the Ith column of the per-atom array
 calculated by the fix is used.  Note that some fixes only produce
 their values on certain timesteps, which must be compatible with
 {Nevery}, else an error results.  Users can also write code for their
 own fix styles and "add them to LAMMPS"_Section_modify.html.
 
 If a value begins with "v_", a variable name must follow which has
 been previously defined in the input script.  Variables of style
 {atom} can reference thermodynamic keywords and various per-atom
 attributes, or invoke other computes, fixes, or variables when they
 are evaluated, so this is a very general means of generating per-atom
 quantities to spatially average.
 
 :line
 
 Additional optional keywords also affect the operation of this fix.
 
 The {units} keyword determines the meaning of the distance units used
 for the layer thickness {delta} and for {origin} if it is a coordinate
 value.  For orthogonal simulation boxes, any of the 3 options may be
 used.  For non-orthogonal (triclinic) simulation boxes, only the
 {reduced} option may be used.
 
 A {box} value selects standard distance units as defined by the
 "units"_units.html command, e.g. Angstroms for units = real or metal.
 A {lattice} value means the distance units are in lattice spacings.
 The "lattice"_lattice.html command must have been previously used to
 define the lattice spacing.  A {reduced} value means normalized
 unitless values between 0 and 1, which represent the lower and upper
 faces of the simulation box respectively.  Thus an {origin} value of
 0.5 means the center of the box in any dimension.  A {delta} value of
 0.1 means 10 layers span the box in any dimension. 
 
 Consider a non-orthogonal box, with layers in the x dimension.  No
 matter how the box is tilted, an {origin} of 0.0 means start layers at
 the lower "b" cross "c" plane of the simulation box and an {origin} of
 1.0 means to start layers at the upper "b" cross "c" face of the box.
 A {delta} value of 0.1 means there will be 10 layers from 0.0 to 1.0,
 regardless of the current size or shape of the simulation box.
 
 The {norm} keyword affects how averaging is done for the output
 produced every {Nfreq} timesteps.  For an {all} setting, a layer
 quantity is summed over all atoms in all {Nrepeat} samples, as is the
 count of atoms in the layer.  The printed value for the layer is
 Total-quantity / Total-count.  In other words it is an average over
 the entire {Nfreq} timescale.
 
 For a {sample} setting, the layer quantity is summed over atoms for
 only a single sample, as is the count, and a "average sample value" is
 computed, i.e. Sample-quantity / Sample-count.  The printed value for
 the layer is the average of the {Nrepeat} "average sample values", In
 other words it is an average of an average.
 
 The {ave} keyword determines how the layer values produced every
 {Nfreq} steps are averaged with layer values produced on previous
 steps that were multiples of {Nfreq}, before they are accessed by
 another output command or written to a file.
 
 If the {ave} setting is {one}, then the layer values produced on
 timesteps that are multiples of {Nfreq} are independent of each other;
 they are output as-is without further averaging.
 
 If the {ave} setting is {running}, then the layer values produced on
 timesteps that are multiples of {Nfreq} are summed and averaged in a
 cumulative sense before being output.  Each output layer value is
 thus the average of the layer value produced on that timestep with all
 preceding values for the same layer.  This running average begins
 when the fix is defined; it can only be restarted by deleting the fix
 via the "unfix"_unfix.html command, or re-defining the fix by
 re-specifying it.
 
 If the {ave} setting is {window}, then the layer values produced on
 timesteps that are multiples of {Nfreq} are summed and averaged within
 a moving "window" of time, so that the last M values for the same
 layer are used to produce the output.  E.g. if M = 3 and Nfreq = 1000,
 then the output on step 10000 will be the average of the individual
 layer values on steps 8000,9000,10000.  Outputs on early steps will
 average over less than M values if they are not available.
 
 The {file} keyword allows a filename to be specified.  Every {Nfreq}
 timesteps, a section of layer info will be written to a text file in
 the following format.  A line with the timestep and number of layers
 is written.  Then one line per layer is written, containing the layer
 ID (1-N), the coordinate of the center of the layer, the number of
 atoms in the layer, and one or more calculated values.  The number of
 values in each line corresponds to the number of values specified in
 the fix ave/spatial command.  The number of atoms and the value(s) are
 average quantities.  If the value of the {units} keyword is {box} or
 {lattice}, the "coord" is printed in box units.  If the value of the
 {units} keyword is {reduced}, the "coord" is printed in reduced units
 (0-1).
 
 The {title1} and {title2} and {title3} keywords allow specification of
 the strings that will be printed as the first 3 lines of the output
 file, assuming the {file} keyword was used.  LAMMPS uses default
 values for each of these, so they do not need to be specified.
 
 By default, these header lines are as follows:
 
 # Spatial-averaged data for fix ID and group name
 # Timestep Number-of-layers
 # Layer Coord Count value1 value2 ... :pre
 
 In the first line, ID and name are replaced with the fix-ID and group
 name.  The second line describes the two values that are printed at
 the first of each section of output.  In the third line the values are
 replaced with the appropriate fields from the fix ave/spatial command.
 
 :line
 
 [Restart, fix_modify, output, run start/stop, minimize info:]
 
 No information about this fix is written to "binary restart
 files"_restart.html.  None of the "fix_modify"_fix_modify.html options
 are relevant to this fix.
 
 This fix computes a global array of values which can be accessed by
 various "output commands"_Section_howto.html#4_15.  The values can
 only be accessed on timesteps that are multiples of {Nfreq} since that
 is when averaging is performed.  The global array has # of rows =
 Nlayers and # of columns = Nvalues+2.  The first column has the layer
 coordinate (center of the layer), the 2nd column has the count of
 atoms in that layer, and the remaining columns are the Nvalue
 quantities.  When the array is accessed with an I that exceeds the
 current number of layers, than a 0.0 is returned by the fix instead of
 an error, since the number of layers can vary as a simulation runs,
 depending on the simulation box size.  The array values calculated by
 this fix are "intensive", since they are already normalized by the
 count of atoms in each layer.
 
 No parameter of this fix can be used with the {start/stop} keywords of
 the "run"_run.html command.  This fix is not invoked during "energy
 minimization"_minimize.html.
 
 [Restrictions:]
 
 When the {ave} keyword is set to {running} or {window} then the number
 of layers must remain the same during the simulation, so that the
 appropriate averaging can be done.  This will be the case if the
 simulation box size doesn't change or if the {units} keyword is set to
 {reduced}.
 
 [Related commands:]
 
 "compute"_compute.html, "fix ave/atom"_fix_ave_atom.html, "fix
 ave/histo"_fix_ave_histo.html, "fix ave/time"_fix_ave_time.html,
 "variable"_variable.html, "fix ave/correlate"_fix_ave_correlate.html,
 
 [Default:]
 
 The option defaults are units = lattice, norm = all, no file output,
 and ave = one, title 1,2,3 = strings as described above.