Homec4science
Diffusion lammps (master) Browse

lammps/src/compute_chunk_atom.cpp
da0dcbe0bb2amaster

Last Change
Raw File

compute_chunk_atom.cpp
View Options

/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
// NOTE: allow for bin center to be variables for sphere/cylinder
#include <mpi.h>
#include <string.h>
#include <stdlib.h>
#include "compute_chunk_atom.h"
#include "atom.h"
#include "update.h"
#include "force.h"
#include "domain.h"
#include "region.h"
#include "lattice.h"
#include "modify.h"
#include "fix_store.h"
#include "comm.h"
#include "group.h"
#include "input.h"
#include "variable.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include <map>
using namespace LAMMPS_NS;
using namespace MathConst;
enum{BIN1D,BIN2D,BIN3D,BINSPHERE,BINCYLINDER,
TYPE,MOLECULE,COMPUTE,FIX,VARIABLE};
enum{LOWER,CENTER,UPPER,COORD};
enum{BOX,LATTICE,REDUCED};
enum{NODISCARD,MIXED,YESDISCARD};
enum{ONCE,NFREQ,EVERY}; // used in several files
enum{LIMITMAX,LIMITEXACT};
#define IDMAX 1024*1024
#define INVOKED_PERATOM 8
/* ---------------------------------------------------------------------- */
ComputeChunkAtom::ComputeChunkAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
chunk_volume_vec(NULL), coord(NULL), ichunk(NULL), chunkID(NULL),
cfvid(NULL), idregion(NULL), region(NULL), cchunk(NULL), fchunk(NULL),
varatom(NULL), id_fix(NULL), fixstore(NULL), lockfix(NULL), chunk(NULL),
exclude(NULL), hash(NULL)
{
if (narg < 4) error->all(FLERR,"Illegal compute chunk/atom command");
peratom_flag = 1;
size_peratom_cols = 0;
create_attribute = 1;
// chunk style and its args
int iarg;
binflag = 0;
ncoord = 0;
cfvid = NULL;
if (strcmp(arg[3],"bin/1d") == 0) {
binflag = 1;
which = BIN1D;
ncoord = 1;
iarg = 4;
readdim(narg,arg,iarg,0);
iarg += 3;
} else if (strcmp(arg[3],"bin/2d") == 0) {
binflag = 1;
which = BIN2D;
ncoord = 2;
iarg = 4;
readdim(narg,arg,iarg,0);
readdim(narg,arg,iarg+3,1);
iarg += 6;
} else if (strcmp(arg[3],"bin/3d") == 0) {
binflag = 1;
which = BIN3D;
ncoord = 3;
iarg = 4;
readdim(narg,arg,iarg,0);
readdim(narg,arg,iarg+3,1);
readdim(narg,arg,iarg+6,2);
iarg += 9;
} else if (strcmp(arg[3],"bin/sphere") == 0) {
binflag = 1;
which = BINSPHERE;
ncoord = 1;
iarg = 4;
if (iarg+6 > narg) error->all(FLERR,"Illegal compute chunk/atom command");
sorigin_user[0] = force->numeric(FLERR,arg[iarg]);
sorigin_user[1] = force->numeric(FLERR,arg[iarg+1]);
sorigin_user[2] = force->numeric(FLERR,arg[iarg+2]);
sradmin_user = force->numeric(FLERR,arg[iarg+3]);
sradmax_user = force->numeric(FLERR,arg[iarg+4]);
nsbin = force->inumeric(FLERR,arg[iarg+5]);
iarg += 6;
} else if (strcmp(arg[3],"bin/cylinder") == 0) {
binflag = 1;
which = BINCYLINDER;
ncoord = 2;
iarg = 4;
readdim(narg,arg,iarg,0);
iarg += 3;
if (dim[0] == 0) {
cdim1 = 1;
cdim2 = 2;
} else if (dim[0] == 1) {
cdim1 = 0;
cdim2 = 2;
} else {
cdim1 = 0;
cdim2 = 1;
}
if (iarg+5 > narg) error->all(FLERR,"Illegal compute chunk/atom command");
corigin_user[dim[0]] = 0.0;
corigin_user[cdim1] = force->numeric(FLERR,arg[iarg]);
corigin_user[cdim2] = force->numeric(FLERR,arg[iarg+1]);
cradmin_user = force->numeric(FLERR,arg[iarg+2]);
cradmax_user = force->numeric(FLERR,arg[iarg+3]);
ncbin = force->inumeric(FLERR,arg[iarg+4]);
iarg += 5;
} else if (strcmp(arg[3],"type") == 0) {
which = TYPE;
iarg = 4;
} else if (strcmp(arg[3],"molecule") == 0) {
which = MOLECULE;
iarg = 4;
} else if (strstr(arg[3],"c_") == arg[3] ||
strstr(arg[3],"f_") == arg[3] ||
strstr(arg[3],"v_") == arg[3]) {
if (arg[3][0] == 'c') which = COMPUTE;
else if (arg[3][0] == 'f') which = FIX;
else if (arg[3][0] == 'v') which = VARIABLE;
iarg = 4;
int n = strlen(arg[3]);
char *suffix = new char[n];
strcpy(suffix,&arg[3][2]);
char *ptr = strchr(suffix,'[');
if (ptr) {
if (suffix[strlen(suffix)-1] != ']')
error->all(FLERR,"Illegal compute chunk/atom command");
argindex = atoi(ptr+1);
*ptr = '\0';
} else argindex = 0;
n = strlen(suffix) + 1;
cfvid = new char[n];
strcpy(cfvid,suffix);
delete [] suffix;
} else error->all(FLERR,"Illegal compute chunk/atom command");
// optional args
regionflag = 0;
idregion = NULL;
nchunksetflag = 0;
nchunkflag = EVERY;
limit = 0;
limitstyle = LIMITMAX;
limitfirst = 0;
idsflag = EVERY;
compress = 0;
int discardsetflag = 0;
discard = MIXED;
minflag[0] = LOWER;
minflag[1] = LOWER;
minflag[2] = LOWER;
maxflag[0] = UPPER;
maxflag[1] = UPPER;
maxflag[2] = UPPER;
scaleflag = LATTICE;
pbcflag = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"region") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal compute chunk/atom command");
int iregion = domain->find_region(arg[iarg+1]);
if (iregion == -1)
error->all(FLERR,"Region ID for compute chunk/atom does not exist");
int n = strlen(arg[iarg+1]) + 1;
idregion = new char[n];
strcpy(idregion,arg[iarg+1]);
regionflag = 1;
iarg += 2;
} else if (strcmp(arg[iarg],"nchunk") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal compute chunk/atom command");
if (strcmp(arg[iarg+1],"once") == 0) nchunkflag = ONCE;
else if (strcmp(arg[iarg+1],"every") == 0) nchunkflag = EVERY;
else error->all(FLERR,"Illegal compute chunk/atom command");
nchunksetflag = 1;
iarg += 2;
} else if (strcmp(arg[iarg],"limit") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal compute chunk/atom command");
limit = force->inumeric(FLERR,arg[iarg+1]);
if (limit < 0) error->all(FLERR,"Illegal compute chunk/atom command");
if (limit && !compress) limitfirst = 1;
iarg += 2;
if (limit) {
if (iarg+1 > narg)
error->all(FLERR,"Illegal compute chunk/atom command");
if (strcmp(arg[iarg+1],"max") == 0) limitstyle = LIMITMAX;
else if (strcmp(arg[iarg+1],"exact") == 0) limitstyle = LIMITEXACT;
else error->all(FLERR,"Illegal compute chunk/atom command");
iarg++;
}
} else if (strcmp(arg[iarg],"ids") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal compute chunk/atom command");
if (strcmp(arg[iarg+1],"once") == 0) idsflag = ONCE;
else if (strcmp(arg[iarg+1],"nfreq") == 0) idsflag = NFREQ;
else if (strcmp(arg[iarg+1],"every") == 0) idsflag = EVERY;
else error->all(FLERR,"Illegal compute chunk/atom command");
iarg += 2;
} else if (strcmp(arg[iarg],"compress") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal compute chunk/atom command");
else if (strcmp(arg[iarg+1],"no") == 0) compress = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) compress = 1;
else error->all(FLERR,"Illegal compute chunk/atom command");
iarg += 2;
} else if (strcmp(arg[iarg],"discard") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal compute chunk/atom command");
if (strcmp(arg[iarg+1],"mixed") == 0) discard = MIXED;
else if (strcmp(arg[iarg+1],"no") == 0) discard = NODISCARD;
else if (strcmp(arg[iarg+1],"yes") == 0) discard = YESDISCARD;
else error->all(FLERR,"Illegal compute chunk/atom command");
discardsetflag = 1;
iarg += 2;
} else if (strcmp(arg[iarg],"bound") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal compute chunk/atom command");
int idim;
if (strcmp(arg[iarg+1],"x") == 0) idim = 0;
else if (strcmp(arg[iarg+1],"y") == 0) idim = 1;
else if (strcmp(arg[iarg+1],"z") == 0) idim = 2;
else error->all(FLERR,"Illegal compute chunk/atom command");
if (strcmp(arg[iarg+2],"lower") == 0) minflag[idim] = LOWER;
else minflag[idim] = COORD;
if (minflag[idim] == COORD)
minvalue[idim] = force->numeric(FLERR,arg[iarg+2]);
if (strcmp(arg[iarg+3],"upper") == 0) maxflag[idim] = UPPER;
else maxflag[idim] = COORD;
if (maxflag[idim] == COORD)
maxvalue[idim] = force->numeric(FLERR,arg[iarg+3]);
else error->all(FLERR,"Illegal compute chunk/atom command");
iarg += 4;
} else if (strcmp(arg[iarg],"units") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal compute chunk/atom command");
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(FLERR,"Illegal compute chunk/atom command");
iarg += 2;
} else if (strcmp(arg[iarg],"pbc") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal compute chunk/atom command");
if (strcmp(arg[iarg+1],"no") == 0) pbcflag = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) pbcflag = 1;
else error->all(FLERR,"Illegal compute chunk/atom command");
iarg += 2;
} else error->all(FLERR,"Illegal compute chunk/atom command");
}
// set nchunkflag and discard to default values if not explicitly set
// for binning style, also check in init() if simulation box is static,
// which sets nchunkflag = ONCE
if (!nchunksetflag) {
if (binflag) {
if (scaleflag == REDUCED) nchunkflag = ONCE;
else nchunkflag = EVERY;
}
if (which == TYPE) nchunkflag = ONCE;
if (which == MOLECULE) {
if (regionflag) nchunkflag = EVERY;
else nchunkflag = ONCE;
}
if (compress) nchunkflag = EVERY;
}
if (!discardsetflag) {
if (binflag) discard = MIXED;
else discard = YESDISCARD;
}
// error checks
if (which == MOLECULE && !atom->molecule_flag)
error->all(FLERR,"Compute chunk/atom molecule for non-molecular system");
if (!binflag && discard == MIXED)
error->all(FLERR,"Compute chunk/atom without bins "
"cannot use discard mixed");
if (which == BIN1D && delta[0] <= 0.0)
error->all(FLERR,"Illegal compute chunk/atom command");
if (which == BIN2D && (delta[0] <= 0.0 || delta[1] <= 0.0))
error->all(FLERR,"Illegal compute chunk/atom command");
if (which == BIN2D && (dim[0] == dim[1]))
error->all(FLERR,"Illegal compute chunk/atom command");
if (which == BIN3D &&
(delta[0] <= 0.0 || delta[1] <= 0.0 || delta[2] <= 0.0))
error->all(FLERR,"Illegal compute chunk/atom command");
if (which == BIN3D &&
(dim[0] == dim[1] || dim[1] == dim[2] || dim[0] == dim[2]))
error->all(FLERR,"Illegal compute chunk/atom command");
if (which == BINSPHERE) {
if (domain->dimension == 2 && sorigin_user[2] != 0.0)
error->all(FLERR,"Compute chunk/atom sphere z origin must be 0.0 for 2d");
if (sradmin_user < 0.0 || sradmin_user >= sradmax_user || nsbin < 1)
error->all(FLERR,"Illegal compute chunk/atom command");
}
if (which == BINCYLINDER) {
if (delta[0] <= 0.0)
error->all(FLERR,"Illegal compute chunk/atom command");
if (domain->dimension == 2 && dim[0] != 2)
error->all(FLERR,"Compute chunk/atom cylinder axis must be z for 2d");
if (cradmin_user < 0.0 || cradmin_user >= cradmax_user || ncbin < 1)
error->all(FLERR,"Illegal compute chunk/atom command");
}
if (which == COMPUTE) {
int icompute = modify->find_compute(cfvid);
if (icompute < 0)
error->all(FLERR,"Compute ID for compute chunk /atom does not exist");
if (modify->compute[icompute]->peratom_flag == 0)
error->all(FLERR,
"Compute chunk/atom compute does not calculate "
"per-atom values");
if (argindex == 0 &&
modify->compute[icompute]->size_peratom_cols != 0)
error->all(FLERR,"Compute chunk/atom compute does not "
"calculate a per-atom vector");
if (argindex && modify->compute[icompute]->size_peratom_cols == 0)
error->all(FLERR,"Compute chunk/atom compute does not "
"calculate a per-atom array");
if (argindex &&
argindex > modify->compute[icompute]->size_peratom_cols)
error->all(FLERR,"Compute chunk/atom compute array is "
"accessed out-of-range");
}
if (which == FIX) {
int ifix = modify->find_fix(cfvid);
if (ifix < 0)
error->all(FLERR,"Fix ID for compute chunk/atom does not exist");
if (modify->fix[ifix]->peratom_flag == 0)
error->all(FLERR,"Compute chunk/atom fix does not calculate "
"per-atom values");
if (argindex == 0 && modify->fix[ifix]->size_peratom_cols != 0)
error->all(FLERR,
"Compute chunk/atom fix does not calculate a per-atom vector");
if (argindex && modify->fix[ifix]->size_peratom_cols == 0)
error->all(FLERR,
"Compute chunk/atom fix does not calculate a per-atom array");
if (argindex && argindex > modify->fix[ifix]->size_peratom_cols)
error->all(FLERR,"Compute chunk/atom fix array is accessed out-of-range");
}
if (which == VARIABLE) {
int ivariable = input->variable->find(cfvid);
if (ivariable < 0)
error->all(FLERR,"Variable name for compute chunk/atom does not exist");
if (input->variable->atomstyle(ivariable) == 0)
error->all(FLERR,"Compute chunk/atom variable is not "
"atom-style variable");
}
// setup scaling
if (binflag) {
if (domain->triclinic == 1 && scaleflag != REDUCED)
error->all(FLERR,"Compute chunk/atom for triclinic boxes "
"requires units reduced");
}
if (scaleflag == LATTICE) {
xscale = domain->lattice->xlattice;
yscale = domain->lattice->ylattice;
zscale = domain->lattice->zlattice;
} else xscale = yscale = zscale = 1.0;
// apply scaling factors and cylinder dims orthogonal to axis
if (binflag) {
double scale;
if (which == BIN1D || which == BIN2D || which == BIN3D ||
which == BINCYLINDER) {
if (which == BIN1D || which == BINCYLINDER) ndim = 1;
if (which == BIN2D) ndim = 2;
if (which == BIN3D) ndim = 3;
for (int idim = 0; idim < ndim; idim++) {
if (dim[idim] == 0) scale = xscale;
else if (dim[idim] == 1) scale = yscale;
else if (dim[idim] == 2) scale = zscale;
delta[idim] *= scale;
invdelta[idim] = 1.0/delta[idim];
if (originflag[idim] == COORD) origin[idim] *= scale;
if (minflag[idim] == COORD) minvalue[idim] *= scale;
if (maxflag[idim] == COORD) maxvalue[idim] *= scale;
}
} else if (which == BINSPHERE) {
sorigin_user[0] *= xscale;
sorigin_user[1] *= yscale;
sorigin_user[2] *= zscale;
sradmin_user *= xscale; // radii are scaled by xscale
sradmax_user *= xscale;
} else if (which == BINCYLINDER) {
if (dim[0] == 0) {
corigin_user[cdim1] *= yscale;
corigin_user[cdim2] *= zscale;
cradmin_user *= yscale; // radii are scaled by first non-axis dim
cradmax_user *= yscale;
} else if (dim[0] == 1) {
corigin_user[cdim1] *= xscale;
corigin_user[cdim2] *= zscale;
cradmin_user *= xscale;
cradmax_user *= xscale;
} else {
corigin_user[cdim1] *= xscale;
corigin_user[cdim2] *= yscale;
cradmin_user *= xscale;
cradmax_user *= xscale;
}
}
}
// initialize chunk vector and per-chunk info
nmax = 0;
chunk = NULL;
nmaxint = -1;
ichunk = NULL;
exclude = NULL;
nchunk = 0;
chunk_volume_scalar = 1.0;
chunk_volume_vec = NULL;
coord = NULL;
chunkID = NULL;
// computeflag = 1 if this compute might invoke another compute
// during assign_chunk_ids()
if (which == COMPUTE || which == FIX || which == VARIABLE) computeflag = 1;
else computeflag = 0;
// other initializations
invoked_setup = -1;
invoked_ichunk = -1;
id_fix = NULL;
fixstore = NULL;
if (compress) hash = new std::map<tagint,int>();
else hash = NULL;
maxvar = 0;
varatom = NULL;
lockcount = 0;
lockfix = NULL;
if (which == MOLECULE) molcheck = 1;
else molcheck = 0;
}
/* ---------------------------------------------------------------------- */
ComputeChunkAtom::~ComputeChunkAtom()
{
// check nfix in case all fixes have already been deleted
if (id_fix && modify->nfix) modify->delete_fix(id_fix);
delete [] id_fix;
memory->destroy(chunk);
memory->destroy(ichunk);
memory->destroy(exclude);
memory->destroy(chunk_volume_vec);
memory->destroy(coord);
memory->destroy(chunkID);
delete [] idregion;
delete [] cfvid;
delete hash;
memory->destroy(varatom);
}
/* ---------------------------------------------------------------------- */
void ComputeChunkAtom::init()
{
// set and check validity of region
if (regionflag) {
int iregion = domain->find_region(idregion);
if (iregion == -1)
error->all(FLERR,"Region ID for compute chunk/atom does not exist");
region = domain->regions[iregion];
}
// set compute,fix,variable
if (which == COMPUTE) {
int icompute = modify->find_compute(cfvid);
if (icompute < 0)
error->all(FLERR,"Compute ID for compute chunk/atom does not exist");
cchunk = modify->compute[icompute];
} else if (which == FIX) {
int ifix = modify->find_fix(cfvid);
if (ifix < 0)
error->all(FLERR,"Fix ID for compute chunk/atom does not exist");
fchunk = modify->fix[ifix];
} else if (which == VARIABLE) {
int ivariable = input->variable->find(cfvid);
if (ivariable < 0)
error->all(FLERR,"Variable name for compute chunk/atom does not exist");
vchunk = ivariable;
}
// for style MOLECULE, check that no mol IDs exceed MAXSMALLINT
// don't worry about group or optional region
if (which == MOLECULE) {
tagint *molecule = atom->molecule;
int nlocal = atom->nlocal;
tagint maxone = -1;
for (int i = 0; i < nlocal; i++)
if (molecule[i] > maxone) maxone = molecule[i];
tagint maxall;
MPI_Allreduce(&maxone,&maxall,1,MPI_LMP_TAGINT,MPI_MAX,world);
if (maxall > MAXSMALLINT)
error->all(FLERR,"Molecule IDs too large for compute chunk/atom");
}
// for binning, if nchunkflag not already set, set it to ONCE or EVERY
// depends on whether simulation box size is static or dynamic
// reset invoked_setup if this is not first run and box just became static
if (binflag && !nchunksetflag && !compress && scaleflag != REDUCED) {
if (domain->box_change_size == 0) {
if (nchunkflag == EVERY && invoked_setup >= 0) invoked_setup = -1;
nchunkflag = ONCE;
} else nchunkflag = EVERY;
}
// require nchunkflag = ONCE if idsflag = ONCE
// b/c nchunk cannot change if chunk IDs are frozen
// can't check until now since nchunkflag may have been adjusted in init()
if (idsflag == ONCE && nchunkflag != ONCE)
error->all(FLERR,"Compute chunk/atom ids once but nchunk is not once");
// create/destroy fix STORE for persistent chunk IDs as needed
// need to do this if idsflag = ONCE or locks will be used by other commands
// need to wait until init() so that fix command(s) are in place
// they increment lockcount if they lock this compute
// fixstore ID = compute-ID + COMPUTE_STORE, fix group = compute group
// fixstore initializes all values to 0.0
if ((idsflag == ONCE || lockcount) && !fixstore) {
int n = strlen(id) + strlen("_COMPUTE_STORE") + 1;
id_fix = new char[n];
strcpy(id_fix,id);
strcat(id_fix,"_COMPUTE_STORE");
char **newarg = new char*[6];
newarg[0] = id_fix;
newarg[1] = group->names[igroup];
newarg[2] = (char *) "STORE";
newarg[3] = (char *) "peratom";
newarg[4] = (char *) "1";
newarg[5] = (char *) "1";
modify->add_fix(6,newarg);
fixstore = (FixStore *) modify->fix[modify->nfix-1];
delete [] newarg;
}
if ((idsflag != ONCE && !lockcount) && fixstore) {
modify->delete_fix(id_fix);
fixstore = NULL;
}
}
/* ----------------------------------------------------------------------
invoke setup_chunks and/or compute_ichunk if only done ONCE
so that nchunks or chunk IDs are assigned when this compute was specified
as opposed to first times compute_peratom() or compute_ichunk() is called
------------------------------------------------------------------------- */
void ComputeChunkAtom::setup()
{
if (nchunkflag == ONCE) setup_chunks();
if (idsflag == ONCE) compute_ichunk();
else invoked_ichunk = -1;
}
/* ----------------------------------------------------------------------
only called by classes that use per-atom computes in standard way
dump, variable, thermo output, other computes, etc
not called by fix chunk or compute chunk commands
they invoke setup_chunks() and compute_ichunk() directly
------------------------------------------------------------------------- */
void ComputeChunkAtom::compute_peratom()
{
invoked_peratom = update->ntimestep;
// grow floating point chunk vector if necessary
if (atom->nmax > nmax) {
memory->destroy(chunk);
nmax = atom->nmax;
memory->create(chunk,nmax,"chunk/atom:chunk");
vector_atom = chunk;
}
setup_chunks();
compute_ichunk();
// copy integer indices into floating-point chunk vector
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) chunk[i] = ichunk[i];
}
/* ----------------------------------------------------------------------
set lock, so that nchunk will not change from startstep to stopstep
called by fix for duration of time it requires lock
OK if called by multiple fix commands
error if all callers do not have same duration
last caller holds the lock, so it can also unlock
stopstep can be positive for final step of finite-size time window
or can be -1 for infinite-size time window
------------------------------------------------------------------------- */
void ComputeChunkAtom::lock(Fix *fixptr, bigint startstep, bigint stopstep)
{
if (lockfix == NULL) {
lockfix = fixptr;
lockstart = startstep;
lockstop = stopstep;
return;
}
if (startstep != lockstart || stopstep != lockstop)
error->all(FLERR,"Two fix commands using "
"same compute chunk/atom command in incompatible ways");
// set lock to last calling Fix, since it will be last to unlock()
lockfix = fixptr;
}
/* ----------------------------------------------------------------------
unset lock
can only be done by fix command that holds the lock
------------------------------------------------------------------------- */
void ComputeChunkAtom::unlock(Fix *fixptr)
{
if (fixptr != lockfix) return;
lockfix = NULL;
}
/* ----------------------------------------------------------------------
assign chunk IDs from 1 to Nchunk to every atom, or 0 if not in chunk
------------------------------------------------------------------------- */
void ComputeChunkAtom::compute_ichunk()
{
int i;
// skip if already done on this step
if (invoked_ichunk == update->ntimestep) return;
// if old IDs persist via storage in fixstore, then just retrieve them
// yes if idsflag = ONCE, and already done once
// or if idsflag = NFREQ and lock is in place and are on later timestep
// else proceed to recalculate per-atom chunk assignments
int restore = 0;
if (idsflag == ONCE && invoked_ichunk >= 0) restore = 1;
if (idsflag == NFREQ && lockfix && update->ntimestep > lockstart) restore = 1;
if (restore) {
invoked_ichunk = update->ntimestep;
double *vstore = fixstore->vstore;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) ichunk[i] = static_cast<int> (vstore[i]);
return;
}
invoked_ichunk = update->ntimestep;
// assign chunk IDs to atoms
// will exclude atoms not in group or in optional region
// already invoked if this is same timestep as last setup_chunks()
if (update->ntimestep > invoked_setup) assign_chunk_ids();
// compress chunk IDs via hash of the original uncompressed IDs
// also apply discard rule except for binning styles which already did
int nlocal = atom->nlocal;
if (compress) {
if (binflag) {
for (i = 0; i < nlocal; i++) {
if (exclude[i]) continue;
if (hash->find(ichunk[i]) == hash->end()) exclude[i] = 1;
else ichunk[i] = hash->find(ichunk[i])->second;
}
} else if (discard == NODISCARD) {
for (i = 0; i < nlocal; i++) {
if (exclude[i]) continue;
if (hash->find(ichunk[i]) == hash->end()) ichunk[i] = nchunk;
else ichunk[i] = hash->find(ichunk[i])->second;
}
} else {
for (i = 0; i < nlocal; i++) {
if (exclude[i]) continue;
if (hash->find(ichunk[i]) == hash->end()) exclude[i] = 1;
else ichunk[i] = hash->find(ichunk[i])->second;
}
}
// else if no compression apply discard rule by itself
} else {
if (discard == NODISCARD) {
for (i = 0; i < nlocal; i++) {
if (exclude[i]) continue;
if (ichunk[i] < 1 || ichunk[i] > nchunk) ichunk[i] = nchunk;;
}
} else {
for (i = 0; i < nlocal; i++) {
if (exclude[i]) continue;
if (ichunk[i] < 1 || ichunk[i] > nchunk) exclude[i] = 1;
}
}
}
// set ichunk = 0 for excluded atoms
// this should set any ichunk values which have not yet been set
for (i = 0; i < nlocal; i++)
if (exclude[i]) ichunk[i] = 0;
// if newly calculated IDs need to persist, store them in fixstore
// yes if idsflag = ONCE or idsflag = NFREQ and lock is in place
if (idsflag == ONCE || (idsflag == NFREQ && lockfix)) {
double *vstore = fixstore->vstore;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) vstore[i] = ichunk[i];
}
// one-time check if which = MOLECULE and
// any chunks do not contain all atoms in the molecule
if (molcheck) {
check_molecules();
molcheck = 0;
}
}
/* ----------------------------------------------------------------------
setup chunks
return nchunk = # of chunks
all atoms will be assigned a chunk ID from 1 to Nchunk, or 0
also setup any internal state needed to quickly assign atoms to chunks
called from compute_peratom() and also directly from
fix chunk and compute chunk commands
------------------------------------------------------------------------- */
int ComputeChunkAtom::setup_chunks()
{
if (invoked_setup == update->ntimestep) return nchunk;
// check if setup needs to be done
// no if lock is in place
// no if nchunkflag = ONCE, and already done once
// otherwise yes
// even if no, check if need to re-compute bin volumes
// so that fix ave/chunk can do proper density normalization
int flag = 0;
if (lockfix) flag = 1;
if (nchunkflag == ONCE && invoked_setup >= 0) flag = 1;
if (flag) {
if (binflag && scaleflag == REDUCED && domain->box_change_size)
bin_volumes();
return nchunk;
}
invoked_setup = update->ntimestep;
// assign chunk IDs to atoms
// will exclude atoms not in group or in optional region
// for binning styles, need to setup bins and their volume first
// else chunk_volume_scalar = entire box volume
// IDs are needed to scan for max ID and for compress()
if (binflag) {
if (which == BIN1D || which == BIN2D || which == BIN3D)
nchunk = setup_xyz_bins();
else if (which == BINSPHERE) nchunk = setup_sphere_bins();
else if (which == BINCYLINDER) nchunk = setup_cylinder_bins();
bin_volumes();
} else {
chunk_volume_scalar = domain->xprd * domain->yprd;
if (domain->dimension == 3) chunk_volume_scalar *= domain->zprd;
}
assign_chunk_ids();
// set nchunk for chunk styles other than binning
// for styles other than TYPE, scan for max ID
if (which == TYPE) nchunk = atom->ntypes;
else if (!binflag) {
int nlocal = atom->nlocal;
int hi = -1;
for (int i = 0; i < nlocal; i++) {
if (exclude[i]) continue;
if (ichunk[i] > hi) hi = ichunk[i];
}
MPI_Allreduce(&hi,&nchunk,1,MPI_INT,MPI_MAX,world);
if (nchunk <= 0) nchunk = 1;
}
// apply limit setting as well as compression of chunks with no atoms
// if limit is set, there are 3 cases:
// no compression, limit specified before compression, or vice versa
if (limit && !binflag) {
if (!compress) {
if (limitstyle == LIMITMAX) nchunk = MIN(nchunk,limit);
else if (limitstyle == LIMITEXACT) nchunk = limit;
} else if (limitfirst) {
nchunk = MIN(nchunk,limit);
}
}
if (compress) compress_chunk_ids();
if (limit && !binflag && compress) {
if (limitstyle == LIMITMAX) nchunk = MIN(nchunk,limit);
else if (limitstyle == LIMITEXACT) nchunk = limit;
}
return nchunk;
}
/* ----------------------------------------------------------------------
assign chunk IDs for all atoms, via ichunk vector
except excluded atoms, their chunk IDs are set to 0 later
also set exclude vector to 0/1 for all atoms
excluded atoms are those not in group or in optional region
called from compute_ichunk() and setup_chunks()
------------------------------------------------------------------------- */
void ComputeChunkAtom::assign_chunk_ids()
{
int i;
// grow integer chunk index vector if necessary
if (atom->nmax > nmaxint) {
memory->destroy(ichunk);
memory->destroy(exclude);
nmaxint = atom->nmax;
memory->create(ichunk,nmaxint,"chunk/atom:ichunk");
memory->create(exclude,nmaxint,"chunk/atom:exclude");
}
// update region if necessary
if (regionflag) region->prematch();
// exclude = 1 if atom is not assigned to a chunk
// exclude atoms not in group or not in optional region
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (regionflag) {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit &&
region->match(x[i][0],x[i][1],x[i][2])) exclude[i] = 0;
else exclude[i] = 1;
}
} else {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) exclude[i] = 0;
else exclude[i] = 1;
}
}
// set ichunk to style value for included atoms
// binning styles apply discard rule, others do not yet
if (binflag) {
if (which == BIN1D) atom2bin1d();
else if (which == BIN2D) atom2bin2d();
else if (which == BIN3D) atom2bin3d();
else if (which == BINSPHERE) atom2binsphere();
else if (which == BINCYLINDER) atom2bincylinder();
} else if (which == TYPE) {
int *type = atom->type;
for (i = 0; i < nlocal; i++) {
if (exclude[i]) continue;
ichunk[i] = type[i];
}
} else if (which == MOLECULE) {
tagint *molecule = atom->molecule;
for (i = 0; i < nlocal; i++) {
if (exclude[i]) continue;
ichunk[i] = static_cast<int> (molecule[i]);
}
} else if (which == COMPUTE) {
if (!(cchunk->invoked_flag & INVOKED_PERATOM)) {
cchunk->compute_peratom();
cchunk->invoked_flag |= INVOKED_PERATOM;
}
if (argindex == 0) {
double *vec = cchunk->vector_atom;
for (i = 0; i < nlocal; i++) {
if (exclude[i]) continue;
ichunk[i] = static_cast<int> (vec[i]);
}
} else {
double **array = cchunk->array_atom;
int argm1 = argindex - 1;
for (i = 0; i < nlocal; i++) {
if (exclude[i]) continue;
ichunk[i] = static_cast<int> (array[i][argm1]);
}
}
} else if (which == FIX) {
if (update->ntimestep % fchunk->peratom_freq)
error->all(FLERR,"Fix used in compute chunk/atom not "
"computed at compatible time");
if (argindex == 0) {
double *vec = fchunk->vector_atom;
for (i = 0; i < nlocal; i++) {
if (exclude[i]) continue;
ichunk[i] = static_cast<int> (vec[i]);
}
} else {
double **array = fchunk->array_atom;
int argm1 = argindex - 1;
for (i = 0; i < nlocal; i++) {
if (exclude[i]) continue;
ichunk[i] = static_cast<int> (array[i][argm1]);
}
}
} else if (which == VARIABLE) {
if (atom->nmax > maxvar) {
maxvar = atom->nmax;
memory->destroy(varatom);
memory->create(varatom,maxvar,"chunk/atom:varatom");
}
input->variable->compute_atom(vchunk,igroup,varatom,1,0);
for (i = 0; i < nlocal; i++) {
if (exclude[i]) continue;
ichunk[i] = static_cast<int> (varatom[i]);
}
}
}
/* ----------------------------------------------------------------------
compress chunk IDs currently assigned to atoms across all processors
by removing those with no atoms assigned
current assignment excludes atoms not in group or in optional region
current Nchunk = max ID
operation:
use hash to store list of populated IDs that I own
add new IDs to populated lists communicated from all other procs
final hash has global list of populated ideas
reset Nchunk = length of global list
called by setup_chunks() when setting Nchunk
remapping of chunk IDs to smaller Nchunk occurs later in compute_ichunk()
------------------------------------------------------------------------- */
void ComputeChunkAtom::compress_chunk_ids()
{
hash->clear();
// put my IDs into hash
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (exclude[i]) continue;
if (hash->find(ichunk[i]) == hash->end()) (*hash)[ichunk[i]] = 0;
}
// n = # of my populated IDs
// nall = n summed across all procs
int n = hash->size();
bigint nbone = n;
bigint nball;
MPI_Allreduce(&nbone,&nball,1,MPI_LMP_BIGINT,MPI_SUM,world);
// create my list of populated IDs
int *list = NULL;
memory->create(list,n,"chunk/atom:list");
n = 0;
std::map<tagint,int>::iterator pos;
for (pos = hash->begin(); pos != hash->end(); ++pos)
list[n++] = pos->first;
// if nall < 1M, just allgather all ID lists on every proc
// else perform ring comm
// add IDs from all procs to my hash
if (nball <= IDMAX) {
// setup for allgatherv
int nprocs = comm->nprocs;
int nall = nball;
int *recvcounts,*displs,*listall;
memory->create(recvcounts,nprocs,"chunk/atom:recvcounts");
memory->create(displs,nprocs,"chunk/atom:displs");
memory->create(listall,nall,"chunk/atom:listall");
MPI_Allgather(&n,1,MPI_INT,recvcounts,1,MPI_INT,world);
displs[0] = 0;
for (int iproc = 1; iproc < nprocs; iproc++)
displs[iproc] = displs[iproc-1] + recvcounts[iproc-1];
// allgatherv acquires list of populated IDs from all procs
MPI_Allgatherv(list,n,MPI_INT,listall,recvcounts,displs,MPI_INT,world);
// add all unique IDs in listall to my hash
for (int i = 0; i < nall; i++)
if (hash->find(listall[i]) == hash->end()) (*hash)[listall[i]] = 0;
// clean up
memory->destroy(recvcounts);
memory->destroy(displs);
memory->destroy(listall);
} else {
comm->ring(n,sizeof(int),list,1,idring,NULL,(void *)this,0);
}
memory->destroy(list);
// nchunk = length of hash containing populated IDs from all procs
nchunk = hash->size();
// reset hash value of each original chunk ID to ordered index
// ordered index = new compressed chunk ID (1 to Nchunk)
// leverages fact that map stores keys in ascending order
// also allocate and set chunkID = list of original chunk IDs
// used by fix ave/chunk and compute property/chunk
memory->destroy(chunkID);
memory->create(chunkID,nchunk,"chunk/atom:chunkID");
n = 0;
for (pos = hash->begin(); pos != hash->end(); ++pos) {
chunkID[n] = pos->first;
(*hash)[pos->first] = ++n;
}
}
/* ----------------------------------------------------------------------
callback from comm->ring()
cbuf = list of N chunk IDs from another proc
loop over the list, add each to my hash
hash ends up storing all unique IDs across all procs
------------------------------------------------------------------------- */
void ComputeChunkAtom::idring(int n, char *cbuf, void *ptr)
{
ComputeChunkAtom *cptr = (ComputeChunkAtom *)ptr;
tagint *list = (tagint *) cbuf;
std::map<tagint,int> *hash = cptr->hash;
for (int i = 0; i < n; i++) (*hash)[list[i]] = 0;
}
/* ----------------------------------------------------------------------
one-time check for which = MOLECULE to check
if each chunk contains all atoms in the molecule
issue warning if not
note that this check is without regard to discard rule
if discard == NODISCARD, there is no easy way to check that all
atoms in an out-of-bounds molecule were added to a chunk,
some could have been excluded by group or region, others not
------------------------------------------------------------------------- */
void ComputeChunkAtom::check_molecules()
{
tagint *molecule = atom->molecule;
int nlocal = atom->nlocal;
int flag = 0;
if (!compress) {
for (int i = 0; i < nlocal; i++) {
if (molecule[i] > 0 && molecule[i] <= nchunk &&
ichunk[i] == 0) flag = 1;
}
} else {
int molid;
for (int i = 0; i < nlocal; i++) {
molid = static_cast<int> (molecule[i]);
if (hash->find(molid) != hash->end() && ichunk[i] == 0) flag = 1;
}
}
int flagall;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
if (flagall && comm->me == 0)
error->warning(FLERR,
"One or more chunks do not contain all atoms in molecule");
}
/* ----------------------------------------------------------------------
setup xyz spatial bins and their extent and coordinates
return nbins = # of bins, will become # of chunks
called from setup_chunks()
------------------------------------------------------------------------- */
int ComputeChunkAtom::setup_xyz_bins()
{
int i,j,k,m,n,idim;
double lo,hi,coord1,coord2;
// lo = bin boundary immediately below boxlo or minvalue
// hi = bin boundary immediately above boxhi or maxvalue
// allocate and initialize arrays based on new bin count
double binlo[3],binhi[3];
if (scaleflag == REDUCED) {
binlo[0] = domain->boxlo_lamda[0];
binlo[1] = domain->boxlo_lamda[1];
binlo[2] = domain->boxlo_lamda[2];
binhi[0] = domain->boxhi_lamda[0];
binhi[1] = domain->boxhi_lamda[1];
binhi[2] = domain->boxhi_lamda[2];
} else {
binlo[0] = domain->boxlo[0];
binlo[1] = domain->boxlo[1];
binlo[2] = domain->boxlo[2];
binhi[0] = domain->boxhi[0];
binhi[1] = domain->boxhi[1];
binhi[2] = domain->boxhi[2];
}
if (minflag[0] == COORD) binlo[0] = minvalue[0];
if (minflag[1] == COORD) binlo[1] = minvalue[1];
if (minflag[2] == COORD) binlo[2] = minvalue[2];
if (maxflag[0] == COORD) binhi[0] = maxvalue[0];
if (maxflag[1] == COORD) binhi[1] = maxvalue[1];
if (maxflag[2] == COORD) binhi[2] = maxvalue[2];
int nbins = 1;
for (m = 0; m < ndim; m++) {
idim = dim[m];
if (originflag[m] == LOWER) origin[m] = binlo[idim];
else if (originflag[m] == UPPER) origin[m] = binhi[idim];
else if (originflag[m] == CENTER)
origin[m] = 0.5 * (binlo[idim] + binhi[idim]);
if (origin[m] < binlo[idim]) {
n = static_cast<int> ((binlo[idim] - origin[m]) * invdelta[m]);
lo = origin[m] + n*delta[m];
} else {
n = static_cast<int> ((origin[m] - binlo[idim]) * invdelta[m]);
lo = origin[m] - n*delta[m];
if (lo > binlo[idim]) lo -= delta[m];
}
if (origin[m] < binhi[idim]) {
n = static_cast<int> ((binhi[idim] - origin[m]) * invdelta[m]);
hi = origin[m] + n*delta[m];
if (hi < binhi[idim]) hi += delta[m];
} else {
n = static_cast<int> ((origin[m] - binhi[idim]) * invdelta[m]);
hi = origin[m] - n*delta[m];
}
if (lo > hi) error->all(FLERR,"Invalid bin bounds in compute chunk/atom");
offset[m] = lo;
nlayers[m] = static_cast<int> ((hi-lo) * invdelta[m] + 0.5);
nbins *= nlayers[m];
}
// allocate and set bin coordinates
memory->destroy(coord);
memory->create(coord,nbins,ndim,"chunk/atom:coord");
if (ndim == 1) {
for (i = 0; i < nlayers[0]; i++)
coord[i][0] = offset[0] + (i+0.5)*delta[0];
} else if (ndim == 2) {
m = 0;
for (i = 0; i < nlayers[0]; i++) {
coord1 = offset[0] + (i+0.5)*delta[0];
for (j = 0; j < nlayers[1]; j++) {
coord[m][0] = coord1;
coord[m][1] = offset[1] + (j+0.5)*delta[1];
m++;
}
}
} else if (ndim == 3) {
m = 0;
for (i = 0; i < nlayers[0]; i++) {
coord1 = offset[0] + (i+0.5)*delta[0];
for (j = 0; j < nlayers[1]; j++) {
coord2 = offset[1] + (j+0.5)*delta[1];
for (k = 0; k < nlayers[2]; k++) {
coord[m][0] = coord1;
coord[m][1] = coord2;
coord[m][2] = offset[2] + (k+0.5)*delta[2];
m++;
}
}
}
}
return nbins;
}
/* ----------------------------------------------------------------------
setup spherical spatial bins and their single coordinate
return nsphere = # of bins, will become # of chunks
called from setup_chunks()
------------------------------------------------------------------------- */
int ComputeChunkAtom::setup_sphere_bins()
{
// convert sorigin_user to sorigin
// sorigin,srad are always in box units, for orthogonal or triclinic domains
// lamda2x works for either orthogonal or triclinic
if (scaleflag == REDUCED) {
domain->lamda2x(sorigin_user,sorigin);
sradmin = sradmin_user * (domain->boxhi[0]-domain->boxlo[0]);
sradmax = sradmax_user * (domain->boxhi[0]-domain->boxlo[0]);
} else {
sorigin[0] = sorigin_user[0];
sorigin[1] = sorigin_user[1];
sorigin[2] = sorigin_user[2];
sradmin = sradmin_user;
sradmax = sradmax_user;
}
// if pbcflag set, sradmax must be < 1/2 box in any periodic dim
// treat orthogonal and triclinic the same
// check every time bins are created
if (pbcflag) {
double *prd_half = domain->prd_half;
int *periodicity = domain->periodicity;
int flag = 0;
if (periodicity[0] && sradmax > prd_half[0]) flag = 1;
if (periodicity[1] && sradmax > prd_half[1]) flag = 1;
if (domain->dimension == 3 &&
periodicity[2] && sradmax > prd_half[2]) flag = 1;
if (flag)
error->all(FLERR,"Compute chunk/atom bin/sphere radius "
"is too large for periodic box");
}
sinvrad = nsbin / (sradmax-sradmin);
// allocate and set bin coordinates
// coord = midpt of radii for a spherical shell
memory->destroy(coord);
memory->create(coord,nsbin,1,"chunk/atom:coord");
double rlo,rhi;
for (int i = 0; i < nsbin; i++) {
rlo = sradmin + i * (sradmax-sradmin) / nsbin;
rhi = sradmin + (i+1) * (sradmax-sradmin) / nsbin;
if (i == nsbin-1) rhi = sradmax;
coord[i][0] = 0.5 * (rlo+rhi);
}
return nsbin;
}
/* ----------------------------------------------------------------------
setup cylindrical spatial bins and their single coordinate
return nsphere = # of bins, will become # of chunks
called from setup_chunks()
------------------------------------------------------------------------- */
int ComputeChunkAtom::setup_cylinder_bins()
{
// setup bins along cylinder axis
// ncplane = # of axis bins
ncplane = setup_xyz_bins();
// convert corigin_user to corigin
// corigin is always in box units, for orthogonal or triclinic domains
// lamda2x works for either orthogonal or triclinic
if (scaleflag == REDUCED) {
domain->lamda2x(corigin_user,corigin);
cradmin = cradmin_user * (domain->boxhi[cdim1]-domain->boxlo[cdim1]);
cradmax = cradmax_user * (domain->boxhi[cdim1]-domain->boxlo[cdim1]);
} else {
corigin[cdim1] = corigin_user[cdim1];
corigin[cdim2] = corigin_user[cdim2];
cradmin = cradmin_user;
cradmax = cradmax_user;
}
// if pbcflag set, sradmax must be < 1/2 box in any periodic non-axis dim
// treat orthogonal and triclinic the same
// check every time bins are created
if (pbcflag) {
double *prd_half = domain->prd_half;
int *periodicity = domain->periodicity;
int flag = 0;
if (periodicity[cdim1] && sradmax > prd_half[cdim1]) flag = 1;
if (periodicity[cdim2] && sradmax > prd_half[cdim2]) flag = 1;
if (flag)
error->all(FLERR,"Compute chunk/atom bin/cylinder radius "
"is too large for periodic box");
}
cinvrad = ncbin / (cradmax-cradmin);
// allocate and set radial bin coordinates
// radial coord = midpt of radii for a cylindrical shell
// axiscoord = saved bin coords along cylndrical axis
// radcoord = saved bin coords in radial direction
double **axiscoord = coord;
memory->create(coord,ncbin,1,"chunk/atom:coord");
double **radcoord = coord;
double rlo,rhi;
for (int i = 0; i < ncbin; i++) {
rlo = cradmin + i * (cradmax-cradmin) / ncbin;
rhi = cradmin + (i+1) * (cradmax-cradmin) / ncbin;
if (i == ncbin-1) rhi = cradmax;
coord[i][0] = 0.5 * (rlo+rhi);
}
// create array of combined coords for all bins
memory->create(coord,ncbin*ncplane,2,"chunk/atom:coord");
int m = 0;
for (int i = 0; i < ncbin; i++)
for (int j = 0; j < ncplane; j++) {
coord[m][0] = radcoord[i][0];
coord[m][1] = axiscoord[j][0];
m++;
}
memory->destroy(axiscoord);
memory->destroy(radcoord);
return ncbin*ncplane;
}
/* ----------------------------------------------------------------------
calculate chunk volumes = bin volumes
scalar if all bins have same volume
vector if per-bin volumes are different
------------------------------------------------------------------------- */
void ComputeChunkAtom::bin_volumes()
{
if (which == BIN1D || which == BIN2D || which == BIN3D) {
if (domain->dimension == 3)
chunk_volume_scalar = domain->xprd * domain->yprd * domain->zprd;
else chunk_volume_scalar = domain->xprd * domain->yprd;
double *prd;
if (scaleflag == REDUCED) prd = domain->prd_lamda;
else prd = domain->prd;
for (int m = 0; m < ndim; m++)
chunk_volume_scalar *= delta[m]/prd[dim[m]];
} else if (which == BINSPHERE) {
memory->destroy(chunk_volume_vec);
memory->create(chunk_volume_vec,nchunk,"chunk/atom:chunk_volume_vec");
double rlo,rhi,vollo,volhi;
for (int i = 0; i < nchunk; i++) {
rlo = sradmin + i * (sradmax-sradmin) / nsbin;
rhi = sradmin + (i+1) * (sradmax-sradmin) / nsbin;
if (i == nchunk-1) rhi = sradmax;
vollo = 4.0/3.0 * MY_PI * rlo*rlo*rlo;
volhi = 4.0/3.0 * MY_PI * rhi*rhi*rhi;
chunk_volume_vec[i] = volhi - vollo;
}
} else if (which == BINCYLINDER) {
memory->destroy(chunk_volume_vec);
memory->create(chunk_volume_vec,nchunk,"chunk/atom:chunk_volume_vec");
// slabthick = delta of bins along cylinder axis
double *prd;
if (scaleflag == REDUCED) prd = domain->prd_lamda;
else prd = domain->prd;
double slabthick = domain->prd[dim[0]] * delta[0]/prd[dim[0]];
// area lo/hi of concentric circles in radial direction
int iradbin;
double rlo,rhi,arealo,areahi;
for (int i = 0; i < nchunk; i++) {
iradbin = i / ncplane;
rlo = cradmin + iradbin * (cradmax-cradmin) / ncbin;
rhi = cradmin + (iradbin+1) * (cradmax-cradmin) / ncbin;
if (iradbin == ncbin-1) rhi = cradmax;
arealo = MY_PI * rlo*rlo;
areahi = MY_PI * rhi*rhi;
chunk_volume_vec[i] = (areahi-arealo) * slabthick;
}
}
}
/* ----------------------------------------------------------------------
assign each atom to a 1d spatial bin (layer)
------------------------------------------------------------------------- */
void ComputeChunkAtom::atom2bin1d()
{
int i,ibin;
double *boxlo,*boxhi,*prd;
double xremap;
double **x = atom->x;
int nlocal = atom->nlocal;
int idim = dim[0];
int nlayer1m1 = nlayers[0] - 1;
int periodicity = domain->periodicity[idim];
if (periodicity) {
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;
}
}
// remap each atom's relevant coord back into box via PBC if necessary
// if scaleflag = REDUCED, box coords -> lamda coords
// apply discard rule
if (scaleflag == REDUCED) domain->x2lamda(nlocal);
for (i = 0; i < nlocal; i++) {
if (exclude[i]) continue;
xremap = x[i][idim];
if (periodicity) {
if (xremap < boxlo[idim]) xremap += prd[idim];
if (xremap >= boxhi[idim]) xremap -= prd[idim];
}
ibin = static_cast<int> ((xremap - offset[0]) * invdelta[0]);
if (xremap < offset[0]) ibin--;
if (discard == MIXED) {
if (!minflag[idim]) ibin = MAX(ibin,0);
else if (ibin < 0) {
exclude[i] = 1;
continue;
}
if (!maxflag[idim]) ibin = MIN(ibin,nlayer1m1);
else if (ibin > nlayer1m1) {
exclude[i] = 1;
continue;
}
} else if (discard == NODISCARD) {
ibin = MAX(ibin,0);
ibin = MIN(ibin,nlayer1m1);
} else if (ibin < 0 || ibin > nlayer1m1) {
exclude[i] = 1;
continue;
}
ichunk[i] = ibin+1;
}
if (scaleflag == REDUCED) domain->lamda2x(nlocal);
}
/* ----------------------------------------------------------------------
assign each atom to a 2d spatial bin (pencil)
------------------------------------------------------------------------- */
void ComputeChunkAtom::atom2bin2d()
{
int i,ibin,i1bin,i2bin;
double *boxlo,*boxhi,*prd;
double xremap,yremap;
double **x = atom->x;
int nlocal = atom->nlocal;
int idim = dim[0];
int jdim = dim[1];
int nlayer1m1 = nlayers[0] - 1;
int nlayer2m1 = nlayers[1] - 1;
int *periodicity = domain->periodicity;
if (periodicity[idim] || periodicity[jdim]) {
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;
}
}
// remap each atom's relevant coord back into box via PBC if necessary
// if scaleflag = REDUCED, box coords -> lamda coords
// apply discard rule
if (scaleflag == REDUCED) domain->x2lamda(nlocal);
for (i = 0; i < nlocal; i++) {
if (exclude[i]) continue;
xremap = x[i][idim];
if (periodicity[idim]) {
if (xremap < boxlo[idim]) xremap += prd[idim];
if (xremap >= boxhi[idim]) xremap -= prd[idim];
}
i1bin = static_cast<int> ((xremap - offset[0]) * invdelta[0]);
if (xremap < offset[0]) i1bin--;
if (discard == MIXED) {
if (!minflag[idim]) i1bin = MAX(i1bin,0);
else if (i1bin < 0) {
exclude[i] = 1;
continue;
}
if (!maxflag[idim]) i1bin = MIN(i1bin,nlayer1m1);
else if (i1bin > nlayer1m1) {
exclude[i] = 1;
continue;
}
} else if (discard == NODISCARD) {
i1bin = MAX(i1bin,0);
i1bin = MIN(i1bin,nlayer1m1);
} else if (i1bin < 0 || i1bin > nlayer1m1) {
exclude[i] = 1;
continue;
}
yremap = x[i][jdim];
if (periodicity[jdim]) {
if (yremap < boxlo[jdim]) yremap += prd[jdim];
if (yremap >= boxhi[jdim]) yremap -= prd[jdim];
}
i2bin = static_cast<int> ((yremap - offset[1]) * invdelta[1]);
if (yremap < offset[1]) i2bin--;
if (discard == MIXED) {
if (!minflag[jdim]) i2bin = MAX(i2bin,0);
else if (i2bin < 0) {
exclude[i] = 1;
continue;
}
if (!maxflag[jdim]) i2bin = MIN(i2bin,nlayer2m1);
else if (i2bin > nlayer2m1) {
exclude[i] = 1;
continue;
}
} else if (discard == NODISCARD) {
i2bin = MAX(i2bin,0);
i2bin = MIN(i2bin,nlayer2m1);
} else if (i2bin < 0 || i2bin > nlayer2m1) {
exclude[i] = 1;
continue;
}
ibin = i1bin*nlayers[1] + i2bin;
ichunk[i] = ibin+1;
}
if (scaleflag == REDUCED) domain->lamda2x(nlocal);
}
/* ----------------------------------------------------------------------
assign each atom to a 3d spatial bin (brick)
------------------------------------------------------------------------- */
void ComputeChunkAtom::atom2bin3d()
{
int i,ibin,i1bin,i2bin,i3bin;
double *boxlo,*boxhi,*prd;
double xremap,yremap,zremap;
double **x = atom->x;
int nlocal = atom->nlocal;
int idim = dim[0];
int jdim = dim[1];
int kdim = dim[2];
int nlayer1m1 = nlayers[0] - 1;
int nlayer2m1 = nlayers[1] - 1;
int nlayer3m1 = nlayers[2] - 1;
int *periodicity = domain->periodicity;
if (periodicity[idim] || periodicity[jdim] || periodicity[kdim]) {
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;
}
}
// remap each atom's relevant coord back into box via PBC if necessary
// if scaleflag = REDUCED, box coords -> lamda coords
// apply discard rule
if (scaleflag == REDUCED) domain->x2lamda(nlocal);
for (i = 0; i < nlocal; i++) {
if (exclude[i]) continue;
xremap = x[i][idim];
if (periodicity[idim]) {
if (xremap < boxlo[idim]) xremap += prd[idim];
if (xremap >= boxhi[idim]) xremap -= prd[idim];
}
i1bin = static_cast<int> ((xremap - offset[0]) * invdelta[0]);
if (xremap < offset[0]) i1bin--;
if (discard == MIXED) {
if (!minflag[idim]) i1bin = MAX(i1bin,0);
else if (i1bin < 0) {
exclude[i] = 1;
continue;
}
if (!maxflag[idim]) i1bin = MIN(i1bin,nlayer1m1);
else if (i1bin > nlayer1m1) {
exclude[i] = 1;
continue;
}
} else if (discard == NODISCARD) {
i1bin = MAX(i1bin,0);
i1bin = MIN(i1bin,nlayer1m1);
} else if (i1bin < 0 || i1bin > nlayer1m1) {
exclude[i] = 1;
continue;
}
yremap = x[i][jdim];
if (periodicity[jdim]) {
if (yremap < boxlo[jdim]) yremap += prd[jdim];
if (yremap >= boxhi[jdim]) yremap -= prd[jdim];
}
i2bin = static_cast<int> ((yremap - offset[1]) * invdelta[1]);
if (yremap < offset[1]) i2bin--;
if (discard == MIXED) {
if (!minflag[jdim]) i2bin = MAX(i2bin,0);
else if (i2bin < 0) {
exclude[i] = 1;
continue;
}
if (!maxflag[jdim]) i2bin = MIN(i2bin,nlayer2m1);
else if (i2bin > nlayer2m1) {
exclude[i] = 1;
continue;
}
} else if (discard == NODISCARD) {
i2bin = MAX(i2bin,0);
i2bin = MIN(i2bin,nlayer2m1);
} else if (i2bin < 0 || i2bin > nlayer2m1) {
exclude[i] = 1;
continue;
}
zremap = x[i][kdim];
if (periodicity[kdim]) {
if (zremap < boxlo[kdim]) zremap += prd[kdim];
if (zremap >= boxhi[kdim]) zremap -= prd[kdim];
}
i3bin = static_cast<int> ((zremap - offset[2]) * invdelta[2]);
if (zremap < offset[2]) i3bin--;
if (discard == MIXED) {
if (!minflag[kdim]) i3bin = MAX(i3bin,0);
else if (i3bin < 0) {
exclude[i] = 1;
continue;
}
if (!maxflag[kdim]) i3bin = MIN(i3bin,nlayer3m1);
else if (i3bin > nlayer3m1) {
exclude[i] = 1;
continue;
}
} else if (discard == NODISCARD) {
i3bin = MAX(i3bin,0);
i3bin = MIN(i3bin,nlayer3m1);
} else if (i3bin < 0 || i3bin > nlayer3m1) {
exclude[i] = 1;
continue;
}
ibin = i1bin*nlayers[1]*nlayers[2] + i2bin*nlayers[2] + i3bin;
ichunk[i] = ibin+1;
}
if (scaleflag == REDUCED) domain->lamda2x(nlocal);
}
/* ----------------------------------------------------------------------
assign each atom to a spherical bin
------------------------------------------------------------------------- */
void ComputeChunkAtom::atom2binsphere()
{
int i,ibin;
double dx,dy,dz,r;
double xremap,yremap,zremap;
double *boxlo = domain->boxlo;
double *boxhi = domain->boxhi;
double *prd = domain->prd;
double *prd_half = domain->prd_half;
int *periodicity = domain->periodicity;
// remap each atom's relevant coords back into box via PBC if necessary
// apply discard rule based on rmin and rmax
double **x = atom->x;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) {
if (exclude[i]) continue;
xremap = x[i][0];
if (periodicity[0]) {
if (xremap < boxlo[0]) xremap += prd[0];
if (xremap >= boxhi[0]) xremap -= prd[0];
}
yremap = x[i][1];
if (periodicity[1]) {
if (xremap < boxlo[1]) yremap += prd[1];
if (xremap >= boxhi[1]) yremap -= prd[1];
}
zremap = x[i][2];
if (periodicity[2]) {
if (xremap < boxlo[2]) zremap += prd[2];
if (xremap >= boxhi[2]) zremap -= prd[2];
}
dx = xremap - sorigin[0];
dy = yremap - sorigin[1];
dz = zremap - sorigin[2];
// if requested, apply PBC to distance from sphere center
// treat orthogonal and triclinic the same
// with dx,dy,dz = lengths independent of each other
// so do not use domain->minimum_image() which couples for triclinic
if (pbcflag) {
if (periodicity[0]) {
if (fabs(dx) > prd_half[0]) {
if (dx < 0.0) dx += prd[0];
else dx -= prd[0];
}
}
if (periodicity[1]) {
if (fabs(dy) > prd_half[1]) {
if (dy < 0.0) dy += prd[1];
else dy -= prd[1];
}
}
if (periodicity[2]) {
if (fabs(dz) > prd_half[2]) {
if (dz < 0.0) dz += prd[2];
else dz -= prd[2];
}
}
}
r = sqrt(dx*dx + dy*dy + dz*dz);
ibin = static_cast<int> ((r - sradmin) * sinvrad);
if (r < sradmin) ibin--;
if (discard == MIXED || discard == NODISCARD) {
ibin = MAX(ibin,0);
ibin = MIN(ibin,nchunk-1);
} else if (ibin < 0 || ibin >= nchunk) {
exclude[i] = 1;
continue;
}
ichunk[i] = ibin+1;
}
}
/* ----------------------------------------------------------------------
assign each atom to a cylindrical bin
------------------------------------------------------------------------- */
void ComputeChunkAtom::atom2bincylinder()
{
int i,rbin,kbin;
double d1,d2,r;
double remap1,remap2;
// first use atom2bin1d() to bin all atoms along cylinder axis
atom2bin1d();
// now bin in radial direction
// kbin = bin along cylinder axis
// rbin = bin in radial direction
double *boxlo = domain->boxlo;
double *boxhi = domain->boxhi;
double *prd = domain->prd;
double *prd_half = domain->prd_half;
int *periodicity = domain->periodicity;
// remap each atom's relevant coords back into box via PBC if necessary
// apply discard rule based on rmin and rmax
double **x = atom->x;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) {
if (exclude[i]) continue;
kbin = ichunk[i] - 1;
remap1 = x[i][cdim1];
if (periodicity[cdim1]) {
if (remap1 < boxlo[cdim1]) remap1 += prd[cdim1];
if (remap1 >= boxhi[cdim1]) remap1 -= prd[cdim1];
}
remap2 = x[i][cdim2];
if (periodicity[cdim2]) {
if (remap2 < boxlo[cdim2]) remap2 += prd[cdim2];
if (remap2 >= boxhi[cdim2]) remap2 -= prd[cdim2];
}
d1 = remap1 - corigin[cdim1];
d2 = remap2 - corigin[cdim2];
// if requested, apply PBC to distance from cylinder axis
// treat orthogonal and triclinic the same
// with d1,d2 = lengths independent of each other
if (pbcflag) {
if (periodicity[cdim1]) {
if (fabs(d1) > prd_half[cdim1]) {
if (d1 < 0.0) d1 += prd[cdim1];
else d1 -= prd[cdim1];
}
}
if (periodicity[cdim2]) {
if (fabs(d2) > prd_half[cdim2]) {
if (d2 < 0.0) d2 += prd[cdim2];
else d2 -= prd[cdim2];
}
}
}
r = sqrt(d1*d1 + d2*d2);
rbin = static_cast<int> ((r - cradmin) * cinvrad);
if (r < cradmin) rbin--;
if (discard == MIXED || discard == NODISCARD) {
rbin = MAX(rbin,0);
rbin = MIN(rbin,ncbin-1);
} else if (rbin < 0 || rbin >= ncbin) {
exclude[i] = 1;
continue;
}
// combine axis and radial bin indices to set ichunk
ichunk[i] = rbin*ncplane + kbin + 1;
}
}
/* ----------------------------------------------------------------------
process args for one dimension of binning info
set dim, originflag, origin, delta
------------------------------------------------------------------------- */
void ComputeChunkAtom::readdim(int narg, char **arg, int iarg, int idim)
{
if (narg < iarg+3) error->all(FLERR,"Illegal compute chunk/atom command");
if (strcmp(arg[iarg],"x") == 0) dim[idim] = 0;
else if (strcmp(arg[iarg],"y") == 0) dim[idim] = 1;
else if (strcmp(arg[iarg],"z") == 0) dim[idim] = 2;
else error->all(FLERR,"Illegal compute chunk/atom command");
if (dim[idim] == 2 && domain->dimension == 2)
error->all(FLERR,"Cannot use compute chunk/atom bin z for 2d model");
if (strcmp(arg[iarg+1],"lower") == 0) originflag[idim] = LOWER;
else if (strcmp(arg[iarg+1],"center") == 0) originflag[idim] = CENTER;
else if (strcmp(arg[iarg+1],"upper") == 0) originflag[idim] = UPPER;
else originflag[idim] = COORD;
if (originflag[idim] == COORD)
origin[idim] = force->numeric(FLERR,arg[iarg+1]);
delta[idim] = force->numeric(FLERR,arg[iarg+2]);
}
/* ----------------------------------------------------------------------
initialize one atom's storage values, called when atom is created
just set chunkID to 0 for new atom
------------------------------------------------------------------------- */
void ComputeChunkAtom::set_arrays(int i)
{
if (!fixstore) return;
double *vstore = fixstore->vstore;
vstore[i] = 0.0;
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays and per-chunk arrays
note: nchunk is actually 0 until first call
------------------------------------------------------------------------- */
double ComputeChunkAtom::memory_usage()
{
double bytes = 2*MAX(nmaxint,0) * sizeof(int); // ichunk,exclude
bytes += nmax * sizeof(double); // chunk
bytes += ncoord*nchunk * sizeof(double); // coord
if (compress) bytes += nchunk * sizeof(int); // chunkID
return bytes;
}
c4science · Help