From e361598ce771347d63f0fc068cf83c207d626336 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Tue, 26 Jun 2007 21:47:46 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@679 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- doc/fix_ave_spatial.html | 48 ++++++++++++++++++++++++------ doc/fix_ave_spatial.txt | 51 +++++++++++++++++++++++++------- src/fix_ave_spatial.cpp | 63 +++++++++++++++++++++++----------------- src/fix_ave_spatial.h | 2 +- 4 files changed, 118 insertions(+), 46 deletions(-) diff --git a/doc/fix_ave_spatial.html b/doc/fix_ave_spatial.html index f34f1bec0f..5e6b9a8698 100644 --- a/doc/fix_ave_spatial.html +++ b/doc/fix_ave_spatial.html @@ -43,7 +43,7 @@ 
keyword = norm or units
   norm value = all or sample
-  units value = box or lattice 
+  units value = box or lattice or reduced
@@ -51,7 +51,7 @@ 

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

Description:

@@ -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.

+

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 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 -value. A box value selects standard distance units as defined by -the units command, e.g. Angstroms for units = real or -metal. A lattice value means the distance units are in lattice -spacings. The lattice command must have been -previously used to define the lattice spacing. +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-orthongal (triclinic) simulation boxes, only the +reduced option may be used. +

+

A box value selects standard distance units as defined by the +units command, e.g. Angstroms for units = real or metal. +A lattice value means the distance units are in lattice spacings. +The lattice 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 @@ -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.

If the density or atom keyword is used, or the compute keyword with a compute that calculates a single quantity per atom, then a diff --git a/doc/fix_ave_spatial.txt b/doc/fix_ave_spatial.txt index 9b15930a0d..f4707aa256 100644 --- a/doc/fix_ave_spatial.txt +++ b/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 -value. 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. +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-orthongal (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 {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. - diff --git a/src/fix_ave_spatial.cpp b/src/fix_ave_spatial.cpp index fe39151532..c527621f32 100644 --- a/src/fix_ave_spatial.cpp +++ b/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; - else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; + if (strcmp(arg[iarg+1],"box") == 0) scaleflag = BOX; + 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 *boxhi = domain->boxhi; + double *boxlo,*boxhi,*prd; + 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) - origin = 0.5 * (boxlo[dim] + boxhi[dim]); + else if (originflag == CENTER) origin = 0.5 * (boxlo[dim] + boxhi[dim]); - if (origin < domain->boxlo[dim]) { - m = static_cast ((domain->boxlo[dim] - origin) * invdelta); + if (origin < boxlo[dim]) { + m = static_cast ((boxlo[dim] - origin) * invdelta); lo = origin + m*delta; } else { - m = static_cast ((origin - domain->boxlo[dim]) * invdelta); + m = static_cast ((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]) { - m = static_cast ((domain->boxhi[dim] - origin) * invdelta); + if (origin < boxhi[dim]) { + m = static_cast ((boxhi[dim] - origin) * invdelta); hi = origin + m*delta; if (hi < boxhi[dim]) hi += delta; } else { - m = static_cast ((origin - domain->boxhi[dim]) * invdelta); + m = static_cast ((origin - boxhi[dim]) * invdelta); hi = origin - m*delta; } offset = lo; nlayers = static_cast ((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) { diff --git a/src/fix_ave_spatial.h b/src/fix_ave_spatial.h index bd2ff9906a..2e23ac7821 100644 --- a/src/fix_ave_spatial.h +++ b/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;