From ac2cf8c4ed91c1bcb39a2160432e34e9fcdde724 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Tue, 30 Aug 2022 12:49:09 -0600 Subject: [PATCH] add support for discard yes/no of out-of-bin atoms --- doc/src/fix_ave_grid.rst | 16 +++++++++- src/fix_ave_grid.cpp | 63 ++++++++++++++++++++++++++++++---------- src/fix_ave_grid.h | 4 +-- src/grid2d.cpp | 8 ++--- src/grid3d.cpp | 12 ++++---- 5 files changed, 74 insertions(+), 29 deletions(-) diff --git a/doc/src/fix_ave_grid.rst b/doc/src/fix_ave_grid.rst index c3a1169acd..7e98f3d619 100644 --- a/doc/src/fix_ave_grid.rst +++ b/doc/src/fix_ave_grid.rst @@ -42,10 +42,13 @@ Syntax f_ID[I] = Ith column of per-grid array calculated by a fix with ID, I can include wildcard (see below) * zero or more keyword/arg pairs may be appended -* keyword = *norm* or *ave* or *bias* or *adof* or *cdof* +* keyword = *discard* or *norm* or *ave* or *bias* or *adof* or *cdof* .. parsed-literal:: + *discard* arg = *yes* or *no* + yes = discard an atom outside grid in a non-periodic dimension + no = remap an atom outside grid in a non-periodic dimension to first or last grid cell *norm* arg = *all* or *sample* or *none* = how output on *Nfreq* steps is normalized all = output is sum of atoms across all *Nrepeat* samples, divided by atom count sample = output is sum of *Nrepeat* sample averages, divided by *Nrepeat* @@ -296,6 +299,17 @@ Additional optional keywords also affect the operation of this fix and its outputs. Some are only applicable to per-atom mode. Some are applicable to both per-atom and per-grid mode. +The *discard* keyword is only applicable to per-atom mode. If a +dimension of the system is non-periodic, then grid cells will only +span the box dimension (fixed or shrink-wrap boundaries as set by the +:doc:`boundary` command). An atom may thus be slighlty outside the +range of grid cells on a particular timestep. If *discard* is set to +*yes* (the default), then the atom will be assigned to the closest +grid cell (lowest or highest) in that dimension. If *discard* is set +to *no* the atom will be ignored. + +---------- + The *norm* keyword is only applicable to per-atom mode. In per-grid mode, the *norm* keyword setting is ignored. The output grid value on an *Nfreq* timestep is the sum of the grid values in each of the diff --git a/src/fix_ave_grid.cpp b/src/fix_ave_grid.cpp index 3482ccea63..391df34b8a 100644 --- a/src/fix_ave_grid.cpp +++ b/src/fix_ave_grid.cpp @@ -36,6 +36,7 @@ using namespace FixConst; enum{ALL,SAMPLE,NONORM}; enum{ONE,RUNNING,WINDOW}; +enum{DISCARD,KEEP}; // OFFSET avoids outside-of-box atoms being rounded to grid pts incorrectly // SHIFT = 0.0 assigns atoms to lower-left grid pt @@ -166,6 +167,7 @@ FixAveGrid::FixAveGrid(LAMMPS *lmp, int narg, char **arg) : // optional args + discardflag = DISCARD; normflag = ALL; aveflag = ONE; nwindow = 0; @@ -175,7 +177,14 @@ FixAveGrid::FixAveGrid(LAMMPS *lmp, int narg, char **arg) : cdof = 0.0; while (iarg < nargnew) { - if (strcmp(arg[iarg],"norm") == 0) { + if (strcmp(arg[iarg],"discard") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/grid command"); + if (strcmp(arg[iarg+1],"yes") == 0) discardflag = DISCARD; + else if (strcmp(arg[iarg+1],"no") == 0) discardflag = KEEP; + else error->all(FLERR,"Illegal fix ave/grid command"); + iarg += 2; + + } else if (strcmp(arg[iarg],"norm") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/grid command"); if (strcmp(arg[iarg+1],"all") == 0) normflag = ALL; else if (strcmp(arg[iarg+1],"sample") == 0) normflag = SAMPLE; @@ -794,10 +803,18 @@ void FixAveGrid::atom2grid() { int i,j,m,n,ix,iy,iz; + double **count2d = grid_sample->count2d; + double **vec2d = grid_sample->vec2d; + double ***array2d = grid_sample->array2d; + double ***count3d = grid_sample->count3d; + double ***vec3d = grid_sample->vec3d; + double ****array3d = grid_sample->array3d; + // bin[i][dim] = indices of bin each atom is in - // not set if group mask does not match - // also count atoms contributing to each bin + // skip atom if group mask does not match // check if any atom is out of bounds for my local grid + // for nonperiodic dim, remap atom to first/last bin if out of bounds + // count atoms contributing to each bin double *boxlo,*prd; int *periodicity = domain->periodicity; @@ -826,13 +843,6 @@ void FixAveGrid::atom2grid() memory->create(skip,maxatom,"ave/grid:skip"); } - double **count2d = grid_sample->count2d; - double **vec2d = grid_sample->vec2d; - double ***array2d = grid_sample->array2d; - double ***count3d = grid_sample->count3d; - double ***vec3d = grid_sample->vec3d; - double ****array3d = grid_sample->array3d; - if (triclinic) domain->x2lamda(nlocal); int flag = 0; @@ -850,12 +860,20 @@ void FixAveGrid::atom2grid() if (ix < nxlo_out || ix > nxhi_out) { if (periodicity[0]) flag = 1; - else skip[i] = 1; + else if (discardflag == KEEP) { + if (ix < nxlo_out && nxlo_out == 0) ix = 0; + else if (ix > nxhi_out && nxhi_out == nxgrid-1) ix = nxgrid-1; + else flag = 1; + } else skip[i] = 1; continue; } if (iy < nylo_out || iy > nyhi_out) { if (periodicity[1]) flag = 1; - else skip[i] = 1; + else if (discardflag == KEEP) { + if (iy < nylo_out && nylo_out == 0) iy = 0; + else if (iy > nyhi_out && nyhi_out == nygrid-1) iy = nygrid-1; + else flag = 1; + } else skip[i] = 1; continue; } @@ -878,17 +896,30 @@ void FixAveGrid::atom2grid() if (ix < nxlo_out || ix > nxhi_out) { if (periodicity[0]) flag = 1; - else skip[i] = 1; - continue; + else if (discardflag == KEEP) { + if (ix < nxlo_out && nxlo_out == 0) ix = 0; + else if (ix > nxhi_out && nxhi_out == nxgrid-1) ix = nxgrid-1; + else flag = 1; + } else skip[i] = 1; } + if (iy < nylo_out || iy > nyhi_out) { if (periodicity[1]) flag = 1; - else skip[i] = 1; + else if (discardflag == KEEP) { + if (iy < nylo_out && nylo_out == 0) iy = 0; + else if (iy > nyhi_out && nyhi_out == nygrid-1) iy = nygrid-1; + else flag = 1; + } else skip[i] = 1; continue; } + if (iz < nzlo_out || iz > nzhi_out) { if (periodicity[2]) flag = 1; - else skip[i] = 1; + else if (discardflag == KEEP) { + if (iz < nzlo_out && nzlo_out == 0) iz = 0; + else if (iz > nzhi_out && nzhi_out == nzgrid-1) iz = nzgrid-1; + else flag = 1; + } else skip[i] = 1; continue; } diff --git a/src/fix_ave_grid.h b/src/fix_ave_grid.h index cd162ae9a5..fe8d2cd98d 100644 --- a/src/fix_ave_grid.h +++ b/src/fix_ave_grid.h @@ -50,8 +50,8 @@ class FixAveGrid : public Fix { int nvalues; int nrepeat, irepeat; bigint nvalid, nvalid_last; - int modeatom,modegrid; - int normflag,aveflag,nwindow; + int modeatom, modegrid; + int discardflag, normflag, aveflag, nwindow; int running_count; int window_count,window_oldest,window_newest; diff --git a/src/grid2d.cpp b/src/grid2d.cpp index 27296d43ff..fe17e84137 100644 --- a/src/grid2d.cpp +++ b/src/grid2d.cpp @@ -129,13 +129,13 @@ Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm, int *periodicity = domain->periodicity; if (!periodicity[0]) { - oxlo = MAX(1,oxlo); - oxhi = MIN(gnx-1,oxhi); + oxlo = MAX(0,oxlo); + oxhi = MIN(nx-1,oxhi); } if (!periodicity[1]) { - oylo = MAX(1,oylo); - oyhi = MIN(gnx-1,oyhi); + oylo = MAX(0,oylo); + oyhi = MIN(ny-1,oyhi); } // error check on size of grid stored by this proc diff --git a/src/grid3d.cpp b/src/grid3d.cpp index a8786b28d8..d16c901cf6 100644 --- a/src/grid3d.cpp +++ b/src/grid3d.cpp @@ -135,18 +135,18 @@ Grid3d::Grid3d(LAMMPS *lmp, MPI_Comm gcomm, int *periodicity = domain->periodicity; if (!periodicity[0]) { - oxlo = MAX(1,oxlo); - oxhi = MIN(gnx-1,oxhi); + oxlo = MAX(0,oxlo); + oxhi = MIN(nx-1,oxhi); } if (!periodicity[1]) { - oylo = MAX(1,oylo); - oyhi = MIN(gnx-1,oyhi); + oylo = MAX(0,oylo); + oyhi = MIN(ny-1,oyhi); } if (!periodicity[2]) { - ozlo = MAX(1,ozlo); - ozhi = MIN(gnx-1,ozhi); + ozlo = MAX(0,ozlo); + ozhi = MIN(nz-1,ozhi); } // error check on size of grid stored by this proc