git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@679 f3b2605a-c512-4ea7-a41b-209d697bcdaa

doc/fix_ave_spatial.html

@@ -43,7 +43,7 @@
 
 <PRE>keyword = <I>norm</I> or <I>units</I>
   <I>norm</I> value = <I>all</I> or <I>sample</I>
-  <I>units</I> value = <I>box</I> or <I>lattice</I> 
+  <I>units</I> value = <I>box</I> or <I>lattice</I> or <I>reduced</I> 
 </PRE>
 
 </UL>
@@ -51,7 +51,7 @@
 </P>
 <PRE>fix 1 all ave/spatial 10000 10000 z lower 2.0 centro.profile compute myCentro
 fix 1 flow ave/spatial 100 1000 y 0.0 1.0 vel.profile atom vx norm sample
-fix 1 flow ave/spatial 100 1000 y 0.0 1.0 dens.profile density mass 
+fix 1 flow ave/spatial 100 1000 y 0.0 2.5 dens.profile density mass 
 </PRE>
 <P><B>Description:</B>
 </P>
@@ -94,13 +94,41 @@ 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 align with 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>.
+</P>
 <P>The <I>units</I> keyword determines the meaning of the distance units used
-for the layer thickness <I>delta</I> and <I>origin</I> if it is a coordinate
+for the layer thickness <I>delta</I> and for <I>origin</I> if it is a coordinate
-value.  A <I>box</I> value selects standard distance units as defined by
+value.  For orthogonal simulation boxes, any of the 3 options may be
-the <A HREF = "units.html">units</A> command, e.g. Angstroms for units = real or
+used.  For non-orthongal (triclinic) simulation boxes, only the
-metal.  A <I>lattice</I> value means the distance units are in lattice
+<I>reduced</I> option may be used.
-spacings.  The <A HREF = "lattice.html">lattice</A> command must have been
+</P>
-previously used to define the lattice spacing.
+<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>Nevery</I> and <I>Nfreq</I> arguments specify how the property calculated
 for each layer is time-averaged.  The property is calculated once each
@@ -126,7 +154,9 @@ 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 atoms and the
-value(s) are average quantities.
+value(s) are average quantities.  Note that for non-orthogonal
+(triclinic) simulation boxes, the "coord" of the layer will be a
+normalized value between 0 and 1.
 </P>
 <P>If the <I>density</I> or <I>atom</I> keyword is used, or the <I>compute</I> keyword
 with a compute that calculates a single quantity per atom, then a

doc/fix_ave_spatial.txt

@@ -30,7 +30,7 @@ style = {density} or {atom} or {compute} :l
 zero or more keyword/value pairs may be appended to the args :l
 keyword = {norm} or {units}
   {norm} value = {all} or {sample}
-  {units} value = {box} or {lattice} :pre
+  {units} value = {box} or {lattice} or {reduced} :pre
   
 :ule
 
@@ -38,7 +38,7 @@ keyword = {norm} or {units}
 
 fix 1 all ave/spatial 10000 10000 z lower 2.0 centro.profile compute myCentro
 fix 1 flow ave/spatial 100 1000 y 0.0 1.0 vel.profile atom vx norm sample
-fix 1 flow ave/spatial 100 1000 y 0.0 1.0 dens.profile density mass :pre
+fix 1 flow ave/spatial 100 1000 y 0.0 2.5 dens.profile density mass :pre
 
 [Description:]
 
@@ -81,13 +81,41 @@ 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 align with 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}.
+
 The {units} keyword determines the meaning of the distance units used
-for the layer thickness {delta} and {origin} if it is a coordinate
+for the layer thickness {delta} and for {origin} if it is a coordinate
-value.  A {box} value selects standard distance units as defined by
+value.  For orthogonal simulation boxes, any of the 3 options may be
-the "units"_units.html command, e.g. Angstroms for units = real or
+used.  For non-orthongal (triclinic) simulation boxes, only the
-metal.  A {lattice} value means the distance units are in lattice
+{reduced} option may be used.
-spacings.  The "lattice"_lattice.html command must have been
+
-previously used to define the lattice spacing.
+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 {Nevery} and {Nfreq} arguments specify how the property calculated
 for each layer is time-averaged.  The property is calculated once each
@@ -113,7 +141,11 @@ 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 atoms and the
-value(s) are average quantities.
+value(s) are average quantities.  For orthogonal simulation boxes,
+"coord" is printed in box units, whether the {units} keyword has a
+value of {box} or {lattice}.  For non-orthogonal (triclinic)
+simulation boxes, "coord" is printed in reduced units (0-1) which is
+the required value for the {units} keyword.
 
 If the {density} or {atom} keyword is used, or the {compute} keyword
 with a compute that calculates a single quantity per atom, then a
@@ -149,4 +181,3 @@ minimization"_minimize.html.
 [Default:]
 
 The option defaults are norm = all and units = lattice.
-

src/fix_ave_spatial.cpp

@@ -33,6 +33,7 @@ using namespace LAMMPS_NS;
 enum{LOWER,CENTER,UPPER,COORD};
 enum{DENSITY_MASS,DENSITY_NUM,VX,VY,VZ,FX,FY,FZ,COMPUTE};
 enum{SAMPLE,ALL};
+enum{BOX,LATTICE,REDUCED};
 
 /* ---------------------------------------------------------------------- */
 
@@ -92,7 +93,7 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) :
   // parse optional args
 
   normflag = ALL;
-  int scaleflag = 1;
+  scaleflag = BOX;
 
   int iarg = 11;
   while (iarg < narg) {
@@ -104,8 +105,9 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) :
       iarg += 2;
     } else if (strcmp(arg[iarg],"units") == 0) {
       if (iarg+2 > narg) error->all("Illegal fix ave/spatial command");
-      if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
+      if (strcmp(arg[iarg+1],"box") == 0) scaleflag = BOX;
-      else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
+      else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = LATTICE;
+      else if (strcmp(arg[iarg+1],"reduced") == 0) scaleflag = REDUCED;
       else error->all("Illegal fix ave/spatial command");
       iarg += 2;
     } else error->all("Illegal fix ave/spatial command");
@@ -119,10 +121,14 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) :
 
   // setup scaling
 
-  if (scaleflag && domain->lattice == NULL)
+  int triclinic = domain->triclinic;
+  if (triclinic == 1 && scaleflag != REDUCED)
+    error->all("Fix ave/spatial for triclinic boxes requires units reduced");
+
+  if (scaleflag == LATTICE && domain->lattice == NULL)
     error->all("Use of fix ave/spatial with undefined lattice");
 
-  if (scaleflag) {
+  if (scaleflag == LATTICE) {
     xscale = domain->lattice->xlattice;
     yscale = domain->lattice->ylattice;
     zscale = domain->lattice->zlattice;
@@ -147,15 +153,6 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) :
   if (delta <= 0.0) error->all("Illegal fix ave/spatial command");
   invdelta = 1.0/delta;
 
-  if (domain->triclinic) {
-    if (dim == 0 && (domain->xy != 0.0 || domain->xz != 0.0))
-      error->all("Cannot (yet) use fix ave/spatial with triclinic box");
-    if (dim == 1 && (domain->xy != 0.0 || domain->yz != 0.0))
-      error->all("Cannot (yet) use fix ave/spatial with triclinic box");
-    if (dim == 2 && (domain->xz != 0.0 || domain->yz != 0.0))
-      error->all("Cannot (yet) use fix ave/spatial with triclinic box");
-  }
-
   nvalues = 1;
   if (which == COMPUTE) {
     int icompute = modify->find_compute(id_compute);
@@ -240,34 +237,42 @@ void FixAveSpatial::end_of_step()
   // allocate and initialize arrays based on new layer count
 
   if (nsum == 0) {
-    double *boxlo = domain->boxlo;
+    double *boxlo,*boxhi,*prd;
-    double *boxhi = domain->boxhi;
+    if (scaleflag == REDUCED) {
+      boxlo = domain->boxlo_lamda;
+      boxhi = domain->boxhi_lamda;
+      prd = domain->prd_lamda;
+    } else {
+      boxlo = domain->boxlo;
+      boxhi = domain->boxhi;
+      prd = domain->prd;
+    }
+
     if (originflag == LOWER) origin = boxlo[dim];
     else if (originflag == UPPER) origin = boxhi[dim];
-    else if (originflag == CENTER)
+    else if (originflag == CENTER) origin = 0.5 * (boxlo[dim] + boxhi[dim]);
-      origin = 0.5 * (boxlo[dim] + boxhi[dim]);
 
-    if (origin < domain->boxlo[dim]) {
+    if (origin < boxlo[dim]) {
-      m = static_cast<int> ((domain->boxlo[dim] - origin) * invdelta);
+      m = static_cast<int> ((boxlo[dim] - origin) * invdelta);
       lo = origin + m*delta;
     } else {
-      m = static_cast<int> ((origin - domain->boxlo[dim]) * invdelta);
+      m = static_cast<int> ((origin - boxlo[dim]) * invdelta);
       lo = origin - m*delta;
-      if (lo > domain->boxlo[dim]) lo -= delta;
+      if (lo > boxlo[dim]) lo -= delta;
     }
-    if (origin < domain->boxhi[dim]) {
+    if (origin < boxhi[dim]) {
-      m = static_cast<int> ((domain->boxhi[dim] - origin) * invdelta);
+      m = static_cast<int> ((boxhi[dim] - origin) * invdelta);
       hi = origin + m*delta;
       if (hi < boxhi[dim]) hi += delta;
     } else {
-      m = static_cast<int> ((origin - domain->boxhi[dim]) * invdelta);
+      m = static_cast<int> ((origin - boxhi[dim]) * invdelta);
       hi = origin - m*delta;
     }
 
     offset = lo;
     nlayers = static_cast<int> ((hi-lo) * invdelta + 0.5);
     double volume = domain->xprd * domain->yprd * domain->zprd;
-    layer_volume = delta * volume/domain->prd[dim];
+    layer_volume = delta * volume/prd[dim];
 
     if (nlayers > maxlayer) {
       maxlayer = nlayers;
@@ -307,6 +312,8 @@ void FixAveSpatial::end_of_step()
 
   // perform the computation for one sample
   // sum within each layer, only include atoms in fix group
+  // insure array index is within bounds (since atoms can be outside box)
+  // if scaleflag = REDUCED, convert box coords to lamda coords
   // DENSITY_MASS adds mass to values
   // DENSITY_NUM adds 1 to values
   // ATOM adds atom vector to values
@@ -316,6 +323,8 @@ void FixAveSpatial::end_of_step()
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
+  if (scaleflag == REDUCED) domain->x2lamda(nlocal);
+
   if (which == DENSITY_MASS) {
     int *type = atom->type;
     double *mass = atom->mass;
@@ -386,6 +395,8 @@ void FixAveSpatial::end_of_step()
     }
   }
 
+  if (scaleflag == REDUCED) domain->lamda2x(nlocal);
+
   // average a single sample
 
   if (normflag == ALL) {

src/fix_ave_spatial.h

@@ -35,7 +35,7 @@ class FixAveSpatial : public Fix {
   char *id_compute;
   FILE *fp;
 
-  int nlayers,nvalues,nsum,maxlayer;
+  int nlayers,nvalues,nsum,maxlayer,scaleflag;
   int compute_size_peratom;
   double xscale,yscale,zscale;
   double layer_volume;