diff --git a/src/compute.cpp b/src/compute.cpp index 1ea1767d92..4beefac23c 100644 --- a/src/compute.cpp +++ b/src/compute.cpp @@ -95,9 +95,7 @@ Compute::Compute(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) ntime = maxtime = 0; tlist = NULL; - // setup map for molecule IDs - - molmap = NULL; + // data masks datamask = ALL_MASK; datamask_ext = ALL_MASK; @@ -114,9 +112,7 @@ Compute::~Compute() { delete [] id; delete [] style; - memory->destroy(tlist); - memory->destroy(molmap); } /* ---------------------------------------------------------------------- */ @@ -230,100 +226,3 @@ void Compute::clearstep() { ntime = 0; } - -/* ---------------------------------------------------------------------- - identify molecule IDs with atoms in group - warn if any atom in group has molecule ID = 0 - warn if any molecule has only some atoms in group - return Ncount = # of molecules with atoms in group - set molmap to NULL if molecule IDs include all in range from 1 to Ncount - else: molecule IDs range from idlo to idhi - set molmap to vector of length idhi-idlo+1 - molmap[id-idlo] = index from 0 to Ncount-1 - return idlo and idhi -------------------------------------------------------------------------- */ - -int Compute::molecules_in_group(tagint &idlo, tagint &idhi) -{ - int i; - - memory->destroy(molmap); - molmap = NULL; - - // find lo/hi molecule ID for any atom in group - // warn if atom in group has ID = 0 - - tagint *molecule = atom->molecule; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - tagint lo = BIG; - tagint hi = -BIG; - int flag = 0; - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (molecule[i] == 0) flag = 1; - lo = MIN(lo,molecule[i]); - hi = MAX(hi,molecule[i]); - } - - int flagall; - MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world); - if (flagall && comm->me == 0) - error->warning(FLERR,"Atom with molecule ID = 0 included in " - "compute molecule group"); - - MPI_Allreduce(&lo,&idlo,1,MPI_LMP_TAGINT,MPI_MIN,world); - MPI_Allreduce(&hi,&idhi,1,MPI_LMP_TAGINT,MPI_MAX,world); - if (idlo == BIG) return 0; - - // molmap = vector of length nlen - // set to 1 for IDs that appear in group across all procs, else 0 - - tagint nlen_tag = idhi-idlo+1; - if (nlen_tag > MAXSMALLINT) - error->all(FLERR,"Too many molecules for compute"); - int nlen = (int) nlen_tag; - - memory->create(molmap,nlen,"compute:molmap"); - for (i = 0; i < nlen; i++) molmap[i] = 0; - - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - molmap[molecule[i]-idlo] = 1; - - int *molmapall; - memory->create(molmapall,nlen,"compute:molmapall"); - MPI_Allreduce(molmap,molmapall,nlen,MPI_INT,MPI_MAX,world); - - // nmolecules = # of non-zero IDs in molmap - // molmap[i] = index of molecule, skipping molecules not in group with -1 - - int nmolecules = 0; - for (i = 0; i < nlen; i++) - if (molmapall[i]) molmap[i] = nmolecules++; - else molmap[i] = -1; - memory->destroy(molmapall); - - // warn if any molecule has some atoms in group and some not in group - - flag = 0; - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) continue; - if (molecule[i] < idlo || molecule[i] > idhi) continue; - if (molmap[molecule[i]-idlo] >= 0) flag = 1; - } - - MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world); - if (flagall && comm->me == 0) - error->warning(FLERR, - "One or more compute molecules has atoms not in group"); - - // if molmap simply stores 1 to Nmolecules, then free it - - if (idlo == 1 && idhi == nmolecules && nlen == nmolecules) { - memory->destroy(molmap); - molmap = NULL; - } - return nmolecules; -} diff --git a/src/compute.h b/src/compute.h index d776f224ee..57eee05aed 100644 --- a/src/compute.h +++ b/src/compute.h @@ -148,10 +148,6 @@ class Compute : protected Pointers { double **vbiasall; // stored velocity bias for all atoms int maxbias; // size of vbiasall array - int *molmap; // convert molecule ID to local index - - int molecules_in_group(tagint &, tagint &); - inline int sbmask(int j) { return j >> SBBITS & 3; } diff --git a/src/compute_atom_molecule.cpp b/src/compute_atom_molecule.cpp deleted file mode 100644 index 8b0d413c23..0000000000 --- a/src/compute_atom_molecule.cpp +++ /dev/null @@ -1,347 +0,0 @@ -/* ---------------------------------------------------------------------- - 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. -------------------------------------------------------------------------- */ - -#include "string.h" -#include "stdlib.h" -#include "compute_atom_molecule.h" -#include "atom.h" -#include "update.h" -#include "modify.h" -#include "compute.h" -#include "fix.h" -#include "input.h" -#include "variable.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; - -enum{COMPUTE,FIX,VARIABLE}; - -#define INVOKED_PERATOM 8 - -/* ---------------------------------------------------------------------- */ - -ComputeAtomMolecule:: -ComputeAtomMolecule(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg) -{ - if (narg < 4) error->all(FLERR,"Illegal compute atom/molecule command"); - - if (atom->molecular == 0) - error->all(FLERR,"Compute atom/molecule requires molecular atom style"); - - // parse args - - which = new int[narg-3]; - argindex = new int[narg-3]; - ids = new char*[narg-3]; - value2index = new int[narg-3]; - nvalues = 0; - - int iarg = 3; - while (iarg < narg) { - if (strncmp(arg[iarg],"c_",2) == 0 || - strncmp(arg[iarg],"f_",2) == 0 || - strncmp(arg[iarg],"v_",2) == 0) { - if (arg[iarg][0] == 'c') which[nvalues] = COMPUTE; - else if (arg[iarg][0] == 'f') which[nvalues] = FIX; - else if (arg[iarg][0] == 'v') which[nvalues] = VARIABLE; - - int n = strlen(arg[iarg]); - char *suffix = new char[n]; - strcpy(suffix,&arg[iarg][2]); - - char *ptr = strchr(suffix,'['); - if (ptr) { - if (suffix[strlen(suffix)-1] != ']') - error->all(FLERR,"Illegal compute reduce command"); - argindex[nvalues] = atoi(ptr+1); - *ptr = '\0'; - } else argindex[nvalues] = 0; - - n = strlen(suffix) + 1; - ids[nvalues] = new char[n]; - strcpy(ids[nvalues],suffix); - nvalues++; - delete [] suffix; - } else error->all(FLERR,"Illegal compute atom/molecule command"); - - iarg++; - } - - // setup and error check - - for (int i = 0; i < nvalues; i++) { - if (which[i] == COMPUTE) { - int icompute = modify->find_compute(ids[i]); - if (icompute < 0) - error->all(FLERR,"Compute ID for compute atom/molecule does not exist"); - if (modify->compute[icompute]->peratom_flag == 0) - error->all(FLERR,"Compute atom/molecule compute does not " - "calculate per-atom values"); - if (argindex[i] == 0 && - modify->compute[icompute]->size_peratom_cols != 0) - error->all(FLERR,"Compute atom/molecule compute does not " - "calculate a per-atom vector"); - if (argindex[i] && modify->compute[icompute]->size_peratom_cols == 0) - error->all(FLERR,"Compute atom/molecule compute does not " - "calculate a per-atom array"); - if (argindex[i] && - argindex[i] > modify->compute[icompute]->size_peratom_cols) - error->all(FLERR,"Compute atom/molecule compute array is " - "accessed out-of-range"); - - } else if (which[i] == FIX) { - int ifix = modify->find_fix(ids[i]); - if (ifix < 0) - error->all(FLERR,"Fix ID for compute atom/molecule does not exist"); - if (modify->fix[ifix]->peratom_flag == 0) - error->all(FLERR,"Compute atom/molecule fix does not " - "calculate per-atom values"); - if (argindex[i] == 0 && - modify->fix[ifix]->size_peratom_cols != 0) - error->all(FLERR,"Compute atom/molecule fix does not " - "calculate a per-atom vector"); - if (argindex[i] && modify->fix[ifix]->size_peratom_cols == 0) - error->all(FLERR,"Compute atom/molecule fix does not " - "calculate a per-atom array"); - if (argindex[i] && - argindex[i] > modify->fix[ifix]->size_peratom_cols) - error->all(FLERR, - "Compute atom/molecule fix array is accessed out-of-range"); - - } else if (which[i] == VARIABLE) { - int ivariable = input->variable->find(ids[i]); - if (ivariable < 0) - error->all(FLERR, - "Variable name for compute atom/molecule does not exist"); - if (input->variable->atomstyle(ivariable) == 0) - error->all(FLERR,"Compute atom/molecule variable is not " - "atom-style variable"); - } - } - - // setup molecule-based data - - nmolecules = molecules_in_group(idlo,idhi); - - vone = vector = NULL; - aone = array = NULL; - - if (nvalues == 1) { - memory->create(vone,nmolecules,"atom/molecule:vone"); - memory->create(vector,nmolecules,"atom/molecule:vector"); - vector_flag = 1; - size_vector = nmolecules; - extvector = 0; - } else { - memory->create(aone,nmolecules,nvalues,"atom/molecule:aone"); - memory->create(array,nmolecules,nvalues,"atom/molecule:array"); - array_flag = 1; - size_array_rows = nmolecules; - size_array_cols = nvalues; - extarray = 0; - } - - maxatom = 0; - scratch = NULL; -} - -/* ---------------------------------------------------------------------- */ - -ComputeAtomMolecule::~ComputeAtomMolecule() -{ - delete [] which; - delete [] argindex; - for (int m = 0; m < nvalues; m++) delete [] ids[m]; - delete [] ids; - delete [] value2index; - - memory->destroy(vone); - memory->destroy(vector); - memory->destroy(aone); - memory->destroy(array); - memory->destroy(scratch); -} - -/* ---------------------------------------------------------------------- */ - -void ComputeAtomMolecule::init() -{ - int ntmp = molecules_in_group(idlo,idhi); - if (ntmp != nmolecules) - error->all(FLERR,"Molecule count changed in compute atom/molecule"); - - // set indices and check validity of all computes,fixes,variables - - for (int m = 0; m < nvalues; m++) { - if (which[m] == COMPUTE) { - int icompute = modify->find_compute(ids[m]); - if (icompute < 0) - error->all(FLERR,"Compute ID for compute atom/molecule does not exist"); - value2index[m] = icompute; - - } else if (which[m] == FIX) { - int ifix = modify->find_fix(ids[m]); - if (ifix < 0) - error->all(FLERR,"Fix ID for compute atom/molecule does not exist"); - value2index[m] = ifix; - - } else if (which[m] == VARIABLE) { - int ivariable = input->variable->find(ids[m]); - if (ivariable < 0) - error->all(FLERR, - "Variable name for compute atom/molecule does not exist"); - value2index[m] = ivariable; - - } else value2index[m] = -1; - } -} - -/* ---------------------------------------------------------------------- */ - -void ComputeAtomMolecule::compute_vector() -{ - int i,j,n; - tagint imol; - - invoked_vector = update->ntimestep; - - for (n = 0; n < nmolecules; n++) vone[n] = 0.0; - compute_one(0); - - tagint *molecule = atom->molecule; - int *mask = atom->mask; - int nlocal = atom->nlocal; - j = 0; - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - imol = molecule[i]; - if (molmap) imol = molmap[imol-idlo]; - else imol--; - vone[imol] += peratom[j]; - } - j += nstride; - } - - int me; - MPI_Comm_rank(world,&me); - MPI_Allreduce(vone,vector,nmolecules,MPI_DOUBLE,MPI_SUM,world); -} - -/* ---------------------------------------------------------------------- */ - -void ComputeAtomMolecule::compute_array() -{ - int i,j,m,n; - tagint imol; - - invoked_array = update->ntimestep; - - tagint *molecule = atom->molecule; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - for (m = 0; m < nvalues; m++) { - for (n = 0; n < nmolecules; n++) aone[n][m] = 0.0; - compute_one(m); - - j = 0; - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - imol = molecule[i]; - if (molmap) imol = molmap[imol-idlo]; - else imol--; - aone[imol][m] += peratom[j]; - } - j += nstride; - } - } - - if (array) - MPI_Allreduce(&aone[0][0],&array[0][0],nvalues*nmolecules, - MPI_DOUBLE,MPI_SUM,world); -} - -/* ---------------------------------------------------------------------- - calculate per-atom values for one input M - invoke the appropriate compute,fix,variable - reallocate scratch if necessary for per-atom variable scratch space -------------------------------------------------------------------------- */ - -void ComputeAtomMolecule::compute_one(int m) -{ - int vidx = value2index[m]; - int aidx = argindex[m]; - - // invoke compute if not previously invoked - - if (which[m] == COMPUTE) { - Compute *compute = modify->compute[vidx]; - - if (!(compute->invoked_flag & INVOKED_PERATOM)) { - compute->compute_peratom(); - compute->invoked_flag |= INVOKED_PERATOM; - } - - if (aidx == 0) { - peratom = compute->vector_atom; - nstride = 1; - } else { - peratom = &compute->array_atom[0][aidx-1]; - nstride = compute->size_peratom_cols; - } - - // access fix fields, check if fix frequency is a match - - } else if (which[m] == FIX) { - if (update->ntimestep % modify->fix[vidx]->peratom_freq) - error->all(FLERR,"Fix used in compute atom/molecule not computed " - "at compatible time"); - Fix *fix = modify->fix[vidx]; - - if (aidx == 0) { - peratom = fix->vector_atom; - nstride = 1; - } else { - peratom = &fix->array_atom[0][aidx-1]; - nstride = fix->size_peratom_cols; - } - - // evaluate atom-style variable - - } else if (which[m] == VARIABLE) { - if (atom->nlocal > maxatom) { - maxatom = atom->nmax; - memory->destroy(scratch); - memory->create(scratch,maxatom,"atom/molecule:scratch"); - } - - peratom = scratch; - input->variable->compute_atom(vidx,igroup,peratom,1,0); - nstride = 1; - } -} - -/* ---------------------------------------------------------------------- - memory usage of local data -------------------------------------------------------------------------- */ - -double ComputeAtomMolecule::memory_usage() -{ - double bytes = (bigint) nmolecules * 2*nvalues * sizeof(double); - if (molmap) bytes += (idhi-idlo+1) * sizeof(int); - bytes += maxatom * sizeof(double); - return bytes; -} diff --git a/src/compute_atom_molecule.h b/src/compute_atom_molecule.h deleted file mode 100644 index 683433e17b..0000000000 --- a/src/compute_atom_molecule.h +++ /dev/null @@ -1,126 +0,0 @@ -/* -*- c++ -*- ---------------------------------------------------------- - 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. -------------------------------------------------------------------------- */ - -#ifdef COMPUTE_CLASS - -ComputeStyle(atom/molecule,ComputeAtomMolecule) - -#else - -#ifndef LMP_COMPUTE_ATOM_MOLECULE_H -#define LMP_COMPUTE_ATOM_MOLECULE_H - -#include "compute.h" - -namespace LAMMPS_NS { - -class ComputeAtomMolecule : public Compute { - public: - ComputeAtomMolecule(class LAMMPS *, int, char **); - ~ComputeAtomMolecule(); - void init(); - void compute_vector(); - void compute_array(); - double memory_usage(); - - private: - int nvalues,nmolecules; - tagint idlo,idhi; - - int *which,*argindex,*value2index; - char **ids; - - int nstride,maxatom; - double *vone; - double **aone; - double *scratch; - double *peratom; - - void compute_one(int); -}; - -} - -#endif -#endif - -/* ERROR/WARNING messages: - -E: Illegal ... command - -Self-explanatory. Check the input script syntax and compare to the -documentation for the command. You can use -echo screen as a -command-line option when running LAMMPS to see the offending line. - -E: Compute atom/molecule requires molecular atom style - -Self-explanatory. - -E: Compute ID for compute atom/molecule does not exist - -Self-explanatory. - -E: Compute atom/molecule compute does not calculate per-atom values - -Self-explanatory. - -E: Compute atom/molecule compute does not calculate a per-atom vector - -Self-explanatory. - -E: Compute atom/molecule compute does not calculate a per-atom array - -Self-explanatory. - -E: Compute atom/molecule compute array is accessed out-of-range - -Self-explanatory. - -E: Fix ID for compute atom/molecule does not exist - -Self-explanatory. - -E: Compute atom/molecule fix does not calculate per-atom values - -Self-explanatory. - -E: Compute atom/molecule fix does not calculate a per-atom vector - -Self-explanatory. - -E: Compute atom/molecule fix does not calculate a per-atom array - -Self-explanatory. - -E: Compute atom/molecule fix array is accessed out-of-range - -Self-explanatory. - -E: Variable name for compute atom/molecule does not exist - -Self-explanatory. - -E: Compute atom/molecule variable is not atom-style variable - -Self-explanatory. - -E: Molecule count changed in compute atom/molecule - -Number of molecules must remain constant over time. - -E: Fix used in compute atom/molecule not computed at compatible time - -The fix must produce per-atom quantities on timesteps that the compute -needs them. - -*/ diff --git a/src/compute_chunk_atom.cpp b/src/compute_chunk_atom.cpp new file mode 100644 index 0000000000..a8982ad18a --- /dev/null +++ b/src/compute_chunk_atom.cpp @@ -0,0 +1,1557 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#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 "memory.h" +#include "error.h" + +#include + +using namespace LAMMPS_NS; + +enum{BIN1D,BIN2D,BIN3D,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 + +// allocate space for static class variable + +ComputeChunkAtom *ComputeChunkAtom::cptr; + +/* ---------------------------------------------------------------------- */ + +ComputeChunkAtom::ComputeChunkAtom(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + 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],"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 fix ave/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; + + 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 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 == BIN3D && + (delta[0] <= 0.0 || delta[1] <= 0.0 || delta[2] <= 0.0)) + 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 + + if (binflag) { + double scale; + if (which == BIN1D) 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; + } + } + + // initialize chunk vector and per-chunk info + + nmax = 0; + chunk = NULL; + nmaxint = 0; + 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(); + 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 (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 wait until init() so that fix ave/chunk 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 (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*[5]; + newarg[0] = id_fix; + newarg[1] = group->names[igroup]; + newarg[2] = (char *) "STORE"; + newarg[3] = (char *) "1"; + newarg[4] = (char *) "1"; + modify->add_fix(5,newarg); + fixstore = (FixStore *) modify->fix[modify->nfix-1]; + delete [] newarg; + } + + if (!lockcount && fixstore) { + delete fixstore; + 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(); +} + +/* ---------------------------------------------------------------------- + only called by classes that use per-atom computes in standard way + dump, variable, thermo output, other computes, etc + not called by fix ave/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->nlocal > 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 ave/chunk for duration of its Nfreq epoch + OK if called by multiple fix ave/chunk commands + error if all callers do not have same duration + last caller holds the lock, so it can also unlock + lockstop 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 ave/chunk 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 ave/chunk 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 (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; + int old = ichunk[i]; + 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 + + int save = 0; + 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 ave/chunk and compute chunk commands +------------------------------------------------------------------------- */ + +int ComputeChunkAtom::setup_chunks() +{ + // check if setup needs to be done + // no if lock is in place + // no if already done on this timestep + // 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 (invoked_setup == update->ntimestep) 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 + // IDs are needed to scan for max ID and for compress() + + if (binflag) { + nchunk = setup_bins(); + bin_volumes(); + } + 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->nlocal > 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 == 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 (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 (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 (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; + int n = nlocal; + for (i = 0; i < nlocal; i++) { + if (exclude[i]) continue; + ichunk[i] = static_cast (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 (array[i][argm1]); + } + } + + } else if (which == VARIABLE) { + if (nlocal > 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 (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::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 { + cptr = this; + comm->ring(n,sizeof(int),list,1,idring,NULL,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) +{ + tagint *list = (tagint *) cbuf; + std::map *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() +{ + int *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 (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 spatial bins and their extent and coordinates + return nbins = # of bins, will become # of chunks + called from setup_chunks() +------------------------------------------------------------------------- */ + +int ComputeChunkAtom::setup_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]; + double *prd; + 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]; + prd = domain->prd_lamda; + } 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]; + prd = domain->prd; + } + + 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 ((binlo[idim] - origin[m]) * invdelta[m]); + lo = origin[m] + n*delta[m]; + } else { + n = static_cast ((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 ((binhi[idim] - origin[m]) * invdelta[m]); + hi = origin[m] + n*delta[m]; + if (hi < binhi[idim]) hi += delta[m]; + } else { + n = static_cast ((origin[m] - binhi[idim]) * invdelta[m]); + hi = origin[m] - n*delta[m]; + } + + if (lo > hi) error->all(FLERR,"Invalid bin bounds in fix ave/spatial"); + + offset[m] = lo; + nlayers[m] = static_cast ((hi-lo) * invdelta[m] + 0.5); + nbins *= nlayers[m]; + chunk_volume_scalar *= delta[m]/prd[idim]; + } + + // 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; +} + +/* ---------------------------------------------------------------------- + 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 { + memory->destroy(chunk_volume_vec); + memory->create(chunk_volume_vec,nchunk,"chunk/atom:chunk_volume_vec"); + // fill in the vector values + } +} + +/* ---------------------------------------------------------------------- + assign each atom to a 1d spatial bin (layer) +------------------------------------------------------------------------- */ + +void ComputeChunkAtom::atom2bin1d() +{ + int i,ibin; + double *boxlo,*boxhi,*prd; + double xremap; + double lamda[3]; + + double **x = atom->x; + int *mask = atom->mask; + 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 ((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 lamda[3]; + + double **x = atom->x; + int *mask = atom->mask; + 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 ((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 ((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 > nlayer1m1) { + 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 lamda[3]; + + double **x = atom->x; + int *mask = atom->mask; + 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 ((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 ((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 > nlayer1m1) { + 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 ((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); +} + +/* ---------------------------------------------------------------------- + process args for one dimension of binning info +------------------------------------------------------------------------- */ + +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; + + 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*nmaxint * sizeof(int); // ichunk,exclude + bytes += nmax * sizeof(double); // chunk + bytes += ncoord*nchunk * sizeof(double); // coord + if (compress) bytes += nchunk * sizeof(int); // chunkID + return bytes; +} diff --git a/src/compute_chunk_atom.h b/src/compute_chunk_atom.h new file mode 100644 index 0000000000..8f69735d92 --- /dev/null +++ b/src/compute_chunk_atom.h @@ -0,0 +1,125 @@ +/* -*- c++ -*- ---------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS + +ComputeStyle(chunk/atom,ComputeChunkAtom) + +#else + +#ifndef LMP_COMPUTE_CHUNK_ATOM_H +#define LMP_COMPUTE_CHUNK_ATOM_H + +#include "compute.h" +#include + +namespace LAMMPS_NS { + +class ComputeChunkAtom : public Compute { + public: + int nchunk,ncoord,compress,idsflag,lockcount; + int computeflag; // 1 if this compute invokes other computes + double chunk_volume_scalar; + double *chunk_volume_vec; + double **coord; + int *ichunk,*chunkID; + + ComputeChunkAtom(class LAMMPS *, int, char **); + ~ComputeChunkAtom(); + void init(); + void setup(); + void compute_peratom(); + void set_arrays(int); + double memory_usage(); + + void lock(class Fix *, bigint, bigint); + void unlock(class Fix *); + int setup_chunks(); + void compute_ichunk(); + + private: + int which,binflag; + int regionflag,nchunksetflag,nchunkflag,discard; + int limit,limitstyle,limitfirst; + int scaleflag; + double xscale,yscale,zscale; + int argindex; + char *cfvid; + + int ndim; + int dim[3],originflag[3],nlayers[3]; + int minflag[3],maxflag[3]; + double origin[3],delta[3]; + double offset[3],invdelta[3]; + double minvalue[3],maxvalue[3]; + + char *idregion; + class Region *region; + + class Compute *cchunk; + class Fix *fchunk; + int vchunk; + int maxvar; + double *varatom; + + char *id_fix; + class FixStore *fixstore; + + class Fix *lockfix; // ptr to FixAveChunk that is locking out setups + // NULL if no lock currently in place + bigint lockstart,lockstop; // timesteps for start and stop of locking + + bigint invoked_setup; // last timestep setup_chunks and nchunk calculated + bigint invoked_ichunk; // last timestep ichunk values calculated + int nmax,nmaxint; + double *chunk; + + int molcheck; // one-time check if all molecule atoms in chunk + int *exclude; // 1 if atom is not assigned to any chunk + std::map *hash; // store original chunks IDs before compression + + // static variable for ring communication callback to access class data + // callback functions for ring communication + + static ComputeChunkAtom *cptr; + static void idring(int, char *); + + void assign_chunk_ids(); + void compress_chunk_ids(); + void check_molecules(); + int setup_bins(); + void bin_volumes(); + void atom2bin1d(); + void atom2bin2d(); + void atom2bin3d(); + void readdim(int, char **, int, int); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +W: More than one compute ke/atom + +It is not efficient to use compute ke/atom more than once. + +*/ diff --git a/src/compute_com_chunk.cpp b/src/compute_com_chunk.cpp new file mode 100644 index 0000000000..b34310bc9c --- /dev/null +++ b/src/compute_com_chunk.cpp @@ -0,0 +1,242 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "string.h" +#include "compute_com_chunk.h" +#include "atom.h" +#include "update.h" +#include "modify.h" +#include "compute_chunk_atom.h" +#include "domain.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +enum{ONCE,NFREQ,EVERY}; + +/* ---------------------------------------------------------------------- */ + +ComputeCOMChunk::ComputeCOMChunk(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg != 4) error->all(FLERR,"Illegal compute com/chunk command"); + + array_flag = 1; + size_array_cols = 3; + size_array_rows = 0; + size_array_rows_variable = 1; + extarray = 0; + + // ID of compute chunk/atom + + int n = strlen(arg[3]) + 1; + idchunk = new char[n]; + strcpy(idchunk,arg[3]); + + init(); + + // chunk-based data + + nchunk = 1; + maxchunk = 0; + massproc = masstotal = NULL; + com = comall = NULL; + allocate(); + + firstflag = massneed = 1; +} + +/* ---------------------------------------------------------------------- */ + +ComputeCOMChunk::~ComputeCOMChunk() +{ + delete [] idchunk; + memory->destroy(massproc); + memory->destroy(masstotal); + memory->destroy(com); + memory->destroy(comall); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeCOMChunk::init() +{ + int icompute = modify->find_compute(idchunk); + if (icompute < 0) + error->all(FLERR,"Chunk/atom compute does not exist for compute com/chunk"); + cchunk = (ComputeChunkAtom *) modify->compute[icompute]; + if (strcmp(cchunk->style,"chunk/atom") != 0) + error->all(FLERR,"Compute com/chunk does not use chunk/atom compute"); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeCOMChunk::setup() +{ + // one-time calculation of per-chunk mass + // done in setup, so that ComputeChunkAtom::setup() is already called + + if (firstflag && cchunk->idsflag == ONCE) { + firstflag = massneed = 0; + compute_array(); + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeCOMChunk::compute_array() +{ + int index; + double massone; + double unwrap[3]; + + invoked_array = update->ntimestep; + + // compute chunk/atom assigns atoms to chunk IDs + // extract ichunk index vector from compute + // ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms + + nchunk = cchunk->setup_chunks(); + cchunk->compute_ichunk(); + int *ichunk = cchunk->ichunk; + + if (nchunk > maxchunk) allocate(); + + // zero local per-chunk values + + for (int i = 0; i < nchunk; i++) + com[i][0] = com[i][1] = com[i][2] = 0.0; + if (massneed) + for (int i = 0; i < nchunk; i++) massproc[i] = 0.0; + + // compute COM for each chunk + + double **x = atom->x; + int *mask = atom->mask; + int *type = atom->type; + imageint *image = atom->image; + double *mass = atom->mass; + double *rmass = atom->rmass; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + index = ichunk[i]-1; + if (index < 0) continue; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + domain->unmap(x[i],image[i],unwrap); + com[index][0] += unwrap[0] * massone; + com[index][1] += unwrap[1] * massone; + com[index][2] += unwrap[2] * massone; + if (massneed) massproc[index] += massone; + } + + MPI_Allreduce(&com[0][0],&comall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world); + if (massneed) + MPI_Allreduce(massproc,masstotal,nchunk,MPI_DOUBLE,MPI_SUM,world); + + for (int i = 0; i < nchunk; i++) { + if (masstotal[i] > 0.0) { + comall[i][0] /= masstotal[i]; + comall[i][1] /= masstotal[i]; + comall[i][2] /= masstotal[i]; + } else comall[i][0] = comall[i][1] = comall[i][2] = 0.0; + } +} + +/* ---------------------------------------------------------------------- + lock methods: called by fix ave/time + these methods insure vector/array size is locked for Nfreq epoch + by passing lock info along to compute chunk/atom +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + increment lock counter +------------------------------------------------------------------------- */ + +void ComputeCOMChunk::lock_enable() +{ + cchunk->lockcount++; +} + +/* ---------------------------------------------------------------------- + decrement lock counter in compute chunk/atom, it if still exists +------------------------------------------------------------------------- */ + +void ComputeCOMChunk::lock_disable() +{ + int icompute = modify->find_compute(idchunk); + if (icompute >= 0) { + cchunk = (ComputeChunkAtom *) modify->compute[icompute]; + cchunk->lockcount--; + } +} + +/* ---------------------------------------------------------------------- + calculate and return # of chunks = length of vector/array +------------------------------------------------------------------------- */ + +int ComputeCOMChunk::lock_length() +{ + nchunk = cchunk->setup_chunks(); + return nchunk; +} + +/* ---------------------------------------------------------------------- + set the lock from startstep to stopstep +------------------------------------------------------------------------- */ + +void ComputeCOMChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep) +{ + cchunk->lock(fixptr,startstep,stopstep); +} + +/* ---------------------------------------------------------------------- + unset the lock +------------------------------------------------------------------------- */ + +void ComputeCOMChunk::unlock(Fix *fixptr) +{ + cchunk->unlock(fixptr); +} + +/* ---------------------------------------------------------------------- + free and reallocate per-chunk arrays +------------------------------------------------------------------------- */ + +void ComputeCOMChunk::allocate() +{ + memory->destroy(massproc); + memory->destroy(masstotal); + memory->destroy(com); + memory->destroy(comall); + size_array_rows = maxchunk = nchunk; + memory->create(massproc,maxchunk,"com/chunk:massproc"); + memory->create(masstotal,maxchunk,"com/chunk:masstotal"); + memory->create(com,maxchunk,3,"com/chunk:com"); + memory->create(comall,maxchunk,3,"com/chunk:comall"); + array = comall; +} + +/* ---------------------------------------------------------------------- + memory usage of local data +------------------------------------------------------------------------- */ + +double ComputeCOMChunk::memory_usage() +{ + double bytes = (bigint) maxchunk * 2 * sizeof(double); + bytes += (bigint) maxchunk * 2*3 * sizeof(double); + return bytes; +} diff --git a/src/compute_com_molecule.h b/src/compute_com_chunk.h similarity index 72% rename from src/compute_com_molecule.h rename to src/compute_com_chunk.h index 207e619670..543a565900 100644 --- a/src/compute_com_molecule.h +++ b/src/compute_com_chunk.h @@ -13,31 +13,43 @@ #ifdef COMPUTE_CLASS -ComputeStyle(com/molecule,ComputeCOMMolecule) +ComputeStyle(com/chunk,ComputeCOMChunk) #else -#ifndef LMP_COMPUTE_COM_MOLECULE_H -#define LMP_COMPUTE_COM_MOLECULE_H +#ifndef LMP_COMPUTE_COM_CHUNK_H +#define LMP_COMPUTE_COM_CHUNK_H #include "compute.h" namespace LAMMPS_NS { -class ComputeCOMMolecule : public Compute { +class ComputeCOMChunk : public Compute { public: - ComputeCOMMolecule(class LAMMPS *, int, char **); - ~ComputeCOMMolecule(); + ComputeCOMChunk(class LAMMPS *, int, char **); + ~ComputeCOMChunk(); void init(); + void setup(); void compute_array(); + + void lock_enable(); + void lock_disable(); + int lock_length(); + void lock(class Fix *, bigint, bigint); + void unlock(class Fix *); + double memory_usage(); private: - int nmolecules; - tagint idlo,idhi; + int nchunk,maxchunk; + int firstflag,massneed; + char *idchunk; + class ComputeChunkAtom *cchunk; double *massproc,*masstotal; double **com,**comall; + + void allocate(); }; } diff --git a/src/compute_com_molecule.cpp b/src/compute_com_molecule.cpp deleted file mode 100644 index 74b2cb5ac3..0000000000 --- a/src/compute_com_molecule.cpp +++ /dev/null @@ -1,148 +0,0 @@ -/* ---------------------------------------------------------------------- - 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. -------------------------------------------------------------------------- */ - -#include "compute_com_molecule.h" -#include "atom.h" -#include "update.h" -#include "domain.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; - -/* ---------------------------------------------------------------------- */ - -ComputeCOMMolecule::ComputeCOMMolecule(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg) -{ - if (narg != 3) error->all(FLERR,"Illegal compute com/molecule command"); - - if (atom->molecular == 0) - error->all(FLERR,"Compute com/molecule requires molecular atom style"); - - array_flag = 1; - size_array_cols = 3; - extarray = 0; - - // setup molecule-based data - - nmolecules = molecules_in_group(idlo,idhi); - size_array_rows = nmolecules; - - memory->create(massproc,nmolecules,"com/molecule:massproc"); - memory->create(masstotal,nmolecules,"com/molecule:masstotal"); - memory->create(com,nmolecules,3,"com/molecule:com"); - memory->create(comall,nmolecules,3,"com/molecule:comall"); - array = comall; - - // compute masstotal for each molecule - - int *mask = atom->mask; - tagint *molecule = atom->molecule; - int *type = atom->type; - double *mass = atom->mass; - double *rmass = atom->rmass; - int nlocal = atom->nlocal; - - tagint imol; - double massone; - - for (int i = 0; i < nmolecules; i++) massproc[i] = 0.0; - - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - imol = molecule[i]; - if (molmap) imol = molmap[imol-idlo]; - else imol--; - massproc[imol] += massone; - } - - MPI_Allreduce(massproc,masstotal,nmolecules,MPI_DOUBLE,MPI_SUM,world); -} - -/* ---------------------------------------------------------------------- */ - -ComputeCOMMolecule::~ComputeCOMMolecule() -{ - memory->destroy(massproc); - memory->destroy(masstotal); - memory->destroy(com); - memory->destroy(comall); -} - -/* ---------------------------------------------------------------------- */ - -void ComputeCOMMolecule::init() -{ - int ntmp = molecules_in_group(idlo,idhi); - if (ntmp != nmolecules) - error->all(FLERR,"Molecule count changed in compute com/molecule"); -} - -/* ---------------------------------------------------------------------- */ - -void ComputeCOMMolecule::compute_array() -{ - tagint imol; - double massone; - double unwrap[3]; - - invoked_array = update->ntimestep; - - for (int i = 0; i < nmolecules; i++) - com[i][0] = com[i][1] = com[i][2] = 0.0; - - double **x = atom->x; - int *mask = atom->mask; - tagint *molecule = atom->molecule; - int *type = atom->type; - imageint *image = atom->image; - double *mass = atom->mass; - double *rmass = atom->rmass; - int nlocal = atom->nlocal; - - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - imol = molecule[i]; - if (molmap) imol = molmap[imol-idlo]; - else imol--; - domain->unmap(x[i],image[i],unwrap); - com[imol][0] += unwrap[0] * massone; - com[imol][1] += unwrap[1] * massone; - com[imol][2] += unwrap[2] * massone; - } - - MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules, - MPI_DOUBLE,MPI_SUM,world); - for (int i = 0; i < nmolecules; i++) { - comall[i][0] /= masstotal[i]; - comall[i][1] /= masstotal[i]; - comall[i][2] /= masstotal[i]; - } -} - -/* ---------------------------------------------------------------------- - memory usage of local data -------------------------------------------------------------------------- */ - -double ComputeCOMMolecule::memory_usage() -{ - double bytes = (bigint) nmolecules * 2 * sizeof(double); - if (molmap) bytes += (idhi-idlo+1) * sizeof(int); - bytes += (bigint) nmolecules * 2*3 * sizeof(double); - return bytes; -} diff --git a/src/compute_gyration_chunk.cpp b/src/compute_gyration_chunk.cpp new file mode 100644 index 0000000000..91bcdf97c7 --- /dev/null +++ b/src/compute_gyration_chunk.cpp @@ -0,0 +1,357 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "math.h" +#include "string.h" +#include "compute_gyration_chunk.h" +#include "atom.h" +#include "update.h" +#include "modify.h" +#include "compute_chunk_atom.h" +#include "domain.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputeGyrationChunk::ComputeGyrationChunk(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg < 4) error->all(FLERR,"Illegal compute gyration/chunk command"); + + // ID of compute chunk/atom + + int n = strlen(arg[6]) + 1; + idchunk = new char[n]; + strcpy(idchunk,arg[6]); + + init(); + + // optional args + + tensor = 0; + int iarg = 3; + while (iarg < narg) { + if (strcmp(arg[iarg],"tensor") == 0) { + tensor = 1; + iarg++; + } else error->all(FLERR,"Illegal compute gyration/chunk command"); + } + + if (tensor) { + array_flag = 1; + size_array_cols = 6; + size_array_rows = 0; + size_array_rows_variable = 1; + extarray = 0; + } else { + vector_flag = 1; + size_vector = 0; + size_vector_variable = 1; + extvector = 0; + } + + // chunk-based data + + nchunk = 1; + maxchunk = 0; + massproc = masstotal = NULL; + com = comall = NULL; + rg = rgall = NULL; + rgt = rgtall = NULL; + allocate(); +} + +/* ---------------------------------------------------------------------- */ + +ComputeGyrationChunk::~ComputeGyrationChunk() +{ + delete [] idchunk; + memory->destroy(massproc); + memory->destroy(masstotal); + memory->destroy(com); + memory->destroy(comall); + memory->destroy(rg); + memory->destroy(rgall); + memory->destroy(rgt); + memory->destroy(rgtall); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeGyrationChunk::init() +{ + int icompute = modify->find_compute(idchunk); + if (icompute < 0) + error->all(FLERR,"Chunk/atom compute does not exist for " + "compute gyration/chunk"); + cchunk = (ComputeChunkAtom *) modify->compute[icompute]; + if (strcmp(cchunk->style,"chunk/atom") != 0) + error->all(FLERR,"Compute gyration/chunk does not use chunk/atom compute"); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeGyrationChunk::compute_vector() +{ + int i,index; + double dx,dy,dz,massone; + double unwrap[3]; + + invoked_array = update->ntimestep; + + com_chunk(); + int *ichunk = cchunk->ichunk; + + for (i = 0; i < nchunk; i++) rg[i] = 0.0; + + // compute Rg for each chunk + + double **x = atom->x; + int *mask = atom->mask; + int *type = atom->type; + imageint *image = atom->image; + double *mass = atom->mass; + double *rmass = atom->rmass; + int nlocal = atom->nlocal; + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + index = ichunk[i]-1; + if (index < 0) continue; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - comall[index][0]; + dy = unwrap[1] - comall[index][1]; + dz = unwrap[2] - comall[index][2]; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + rg[index] += (dx*dx + dy*dy + dz*dz) * massone; + } + + MPI_Allreduce(rg,rgall,nchunk,MPI_DOUBLE,MPI_SUM,world); + + for (int i = 0; i < nchunk; i++) + rgall[i] = sqrt(rgall[i]/masstotal[i]); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeGyrationChunk::compute_array() +{ + int i,j,index; + double dx,dy,dz,massone; + double unwrap[3]; + + invoked_array = update->ntimestep; + + com_chunk(); + int *ichunk = cchunk->ichunk; + + for (i = 0; i < nchunk; i++) + for (j = 0; j < 6; j++) rgt[i][j] = 0.0; + + double **x = atom->x; + int *mask = atom->mask; + int *type = atom->type; + imageint *image = atom->image; + double *mass = atom->mass; + double *rmass = atom->rmass; + int nlocal = atom->nlocal; + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + index = ichunk[i]-1; + if (index < 0) continue; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - comall[index][0]; + dy = unwrap[1] - comall[index][1]; + dz = unwrap[2] - comall[index][2]; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + rgt[index][0] += dx*dx * massone; + rgt[index][1] += dy*dy * massone; + rgt[index][2] += dz*dz * massone; + rgt[index][3] += dx*dy * massone; + rgt[index][4] += dx*dz * massone; + rgt[index][5] += dy*dz * massone; + } + + if (nchunk) + MPI_Allreduce(&rgt[0][0],&rgtall[0][0],nchunk*6,MPI_DOUBLE,MPI_SUM,world); + + for (i = 0; i < nchunk; i++) + for (j = 0; j < 6; j++) + rgtall[i][j] = rgtall[i][j]/masstotal[i]; +} + + +/* ---------------------------------------------------------------------- + calculate per-chunk COM, used by both scalar and tensor +------------------------------------------------------------------------- */ + +void ComputeGyrationChunk::com_chunk() +{ + int index; + double massone; + double unwrap[3]; + + // compute chunk/atom assigns atoms to chunk IDs + // extract ichunk index vector from compute + // ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms + + nchunk = cchunk->setup_chunks(); + cchunk->compute_ichunk(); + int *ichunk = cchunk->ichunk; + + if (nchunk > maxchunk) allocate(); + + // zero local per-chunk values + + for (int i = 0; i < nchunk; i++) { + massproc[i] = 0.0; + com[i][0] = com[i][1] = com[i][2] = 0.0; + } + + // compute COM for each chunk + + double **x = atom->x; + int *mask = atom->mask; + int *type = atom->type; + imageint *image = atom->image; + double *mass = atom->mass; + double *rmass = atom->rmass; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + index = ichunk[i]-1; + if (index < 0) continue; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + domain->unmap(x[i],image[i],unwrap); + massproc[index] += massone; + com[index][0] += unwrap[0] * massone; + com[index][1] += unwrap[1] * massone; + com[index][2] += unwrap[2] * massone; + } + + MPI_Allreduce(massproc,masstotal,nchunk,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&com[0][0],&comall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world); + + for (int i = 0; i < nchunk; i++) { + comall[i][0] /= masstotal[i]; + comall[i][1] /= masstotal[i]; + comall[i][2] /= masstotal[i]; + } +} + +/* ---------------------------------------------------------------------- + lock methods: called by fix ave/time + these methods insure vector/array size is locked for Nfreq epoch + by passing lock info along to compute chunk/atom +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + increment lock counter +------------------------------------------------------------------------- */ + +void ComputeGyrationChunk::lock_enable() +{ + cchunk->lockcount++; +} + +/* ---------------------------------------------------------------------- + decrement lock counter in compute chunk/atom, it if still exists +------------------------------------------------------------------------- */ + +void ComputeGyrationChunk::lock_disable() +{ + int icompute = modify->find_compute(idchunk); + if (icompute >= 0) { + cchunk = (ComputeChunkAtom *) modify->compute[icompute]; + cchunk->lockcount--; + } +} + +/* ---------------------------------------------------------------------- + calculate and return # of chunks = length of vector/array +------------------------------------------------------------------------- */ + +int ComputeGyrationChunk::lock_length() +{ + nchunk = cchunk->setup_chunks(); + return nchunk; +} + +/* ---------------------------------------------------------------------- + set the lock from startstep to stopstep +------------------------------------------------------------------------- */ + +void ComputeGyrationChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep) +{ + cchunk->lock(fixptr,startstep,stopstep); +} + +/* ---------------------------------------------------------------------- + unset the lock +------------------------------------------------------------------------- */ + +void ComputeGyrationChunk::unlock(Fix *fixptr) +{ + cchunk->unlock(fixptr); +} + +/* ---------------------------------------------------------------------- + free and reallocate per-chunk arrays +------------------------------------------------------------------------- */ + +void ComputeGyrationChunk::allocate() +{ + memory->destroy(massproc); + memory->destroy(masstotal); + memory->destroy(com); + memory->destroy(comall); + memory->destroy(rg); + memory->destroy(rgall); + memory->destroy(rgt); + memory->destroy(rgtall); + size_array_rows = maxchunk = nchunk; + memory->create(massproc,maxchunk,"gyration/chunk:massproc"); + memory->create(masstotal,maxchunk,"gyration/chunk:masstotal"); + memory->create(com,maxchunk,3,"gyration/chunk:com"); + memory->create(comall,maxchunk,3,"gyration/chunk:comall"); + if (tensor) { + memory->create(rgt,maxchunk,6,"gyration/chunk:rgt"); + memory->create(rgtall,maxchunk,6,"gyration/chunk:rgtall"); + array = rgtall; + } else { + memory->create(rg,maxchunk,"gyration/chunk:rg"); + memory->create(rgall,maxchunk,"gyration/chunk:rgall"); + vector = rgall; + } +} + +/* ---------------------------------------------------------------------- + memory usage of local data +------------------------------------------------------------------------- */ + +double ComputeGyrationChunk::memory_usage() +{ + double bytes = (bigint) maxchunk * 2 * sizeof(double); + bytes += (bigint) maxchunk * 2*3 * sizeof(double); + if (tensor) bytes += (bigint) maxchunk * 2*6 * sizeof(double); + else bytes += (bigint) maxchunk * 2 * sizeof(double); + return bytes; +} diff --git a/src/compute_gyration_molecule.h b/src/compute_gyration_chunk.h similarity index 71% rename from src/compute_gyration_molecule.h rename to src/compute_gyration_chunk.h index 5fec77db56..85ed59e14f 100644 --- a/src/compute_gyration_molecule.h +++ b/src/compute_gyration_chunk.h @@ -13,37 +13,47 @@ #ifdef COMPUTE_CLASS -ComputeStyle(gyration/molecule,ComputeGyrationMolecule) +ComputeStyle(gyration/chunk,ComputeGyrationChunk) #else -#ifndef LMP_COMPUTE_GYRATION_MOLECULE_H -#define LMP_COMPUTE_GYRATION_MOLECULE_H +#ifndef LMP_COMPUTE_GYRATION_CHUNK_H +#define LMP_COMPUTE_GYRATION_CHUNK_H #include "compute.h" namespace LAMMPS_NS { -class ComputeGyrationMolecule : public Compute { +class ComputeGyrationChunk : public Compute { public: - ComputeGyrationMolecule(class LAMMPS *, int, char **); - ~ComputeGyrationMolecule(); + ComputeGyrationChunk(class LAMMPS *, int, char **); + ~ComputeGyrationChunk(); void init(); void compute_vector(); void compute_array(); + + void lock_enable(); + void lock_disable(); + int lock_length(); + void lock(class Fix *, bigint, bigint); + void unlock(class Fix *); + double memory_usage(); private: + int nchunk,maxchunk; + char *idchunk; + class ComputeChunkAtom *cchunk; + int tensor; - int nmolecules; - tagint idlo,idhi; double *massproc,*masstotal; double **com,**comall; - double *rg; - double **rgt; + double *rg,*rgall; + double **rgt,**rgtall; - void molcom(); + void com_chunk(); + void allocate(); }; } diff --git a/src/compute_gyration_molecule.cpp b/src/compute_gyration_molecule.cpp deleted file mode 100644 index cd36fe2021..0000000000 --- a/src/compute_gyration_molecule.cpp +++ /dev/null @@ -1,275 +0,0 @@ -/* ---------------------------------------------------------------------- - 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. -------------------------------------------------------------------------- */ - -#include "math.h" -#include "string.h" -#include "compute_gyration_molecule.h" -#include "atom.h" -#include "update.h" -#include "domain.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; - -/* ---------------------------------------------------------------------- */ - -ComputeGyrationMolecule::ComputeGyrationMolecule(LAMMPS *lmp, - int narg, char **arg) : - Compute(lmp, narg, arg) -{ - if (narg < 3) error->all(FLERR,"Illegal compute gyration/molecule command"); - - if (atom->molecular == 0) - error->all(FLERR,"Compute gyration/molecule requires molecular atom style"); - - tensor = 0; - - int iarg = 3; - while (iarg < narg) { - if (strcmp(arg[iarg],"tensor") == 0) { - tensor = 1; - iarg++; - } else error->all(FLERR,"Illegal compute gyration/molecule command"); - } - - // setup molecule-based data - - nmolecules = molecules_in_group(idlo,idhi); - - memory->create(massproc,nmolecules,"gyration/molecule:massproc"); - memory->create(masstotal,nmolecules,"gyration/molecule:masstotal"); - memory->create(com,nmolecules,3,"gyration/molecule:com"); - memory->create(comall,nmolecules,3,"gyration/molecule:comall"); - - rg = vector = NULL; - rgt = array = NULL; - - if (tensor) { - memory->create(rgt,nmolecules,6,"gyration/molecule:rgt"); - memory->create(array,nmolecules,6,"gyration/molecule:array"); - array_flag = 1; - size_array_rows = nmolecules; - size_array_cols = 6; - extarray = 0; - } else { - memory->create(rg,nmolecules,"gyration/molecule:rg"); - memory->create(vector,nmolecules,"gyration/molecule:vector"); - vector_flag = 1; - size_vector = nmolecules; - extvector = 0; - } - - // compute masstotal for each molecule - - int *mask = atom->mask; - tagint *molecule = atom->molecule; - int *type = atom->type; - double *mass = atom->mass; - double *rmass = atom->rmass; - int nlocal = atom->nlocal; - - tagint imol; - double massone; - - for (int i = 0; i < nmolecules; i++) massproc[i] = 0.0; - - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - imol = molecule[i]; - if (molmap) imol = molmap[imol-idlo]; - else imol--; - massproc[imol] += massone; - } - - MPI_Allreduce(massproc,masstotal,nmolecules,MPI_DOUBLE,MPI_SUM,world); -} - -/* ---------------------------------------------------------------------- */ - -ComputeGyrationMolecule::~ComputeGyrationMolecule() -{ - memory->destroy(massproc); - memory->destroy(masstotal); - memory->destroy(com); - memory->destroy(comall); - memory->destroy(rg); - memory->destroy(rgt); -} - -/* ---------------------------------------------------------------------- */ - -void ComputeGyrationMolecule::init() -{ - int ntmp = molecules_in_group(idlo,idhi); - if (ntmp != nmolecules) - error->all(FLERR,"Molecule count changed in compute gyration/molecule"); -} - -/* ---------------------------------------------------------------------- */ - -void ComputeGyrationMolecule::compute_vector() -{ - tagint imol; - double dx,dy,dz,massone; - double unwrap[3]; - - invoked_array = update->ntimestep; - - molcom(); - - for (int i = 0; i < nmolecules; i++) rg[i] = 0.0; - - double **x = atom->x; - int *mask = atom->mask; - tagint *molecule = atom->molecule; - int *type = atom->type; - imageint *image = atom->image; - double *mass = atom->mass; - double *rmass = atom->rmass; - int nlocal = atom->nlocal; - - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - imol = molecule[i]; - if (molmap) imol = molmap[imol-idlo]; - else imol--; - domain->unmap(x[i],image[i],unwrap); - dx = unwrap[0] - comall[imol][0]; - dy = unwrap[1] - comall[imol][1]; - dz = unwrap[2] - comall[imol][2]; - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - rg[imol] += (dx*dx + dy*dy + dz*dz) * massone; - } - - MPI_Allreduce(rg,vector,nmolecules,MPI_DOUBLE,MPI_SUM,world); - - for (int i = 0; i < nmolecules; i++) - vector[i] = sqrt(vector[i]/masstotal[i]); -} - -/* ---------------------------------------------------------------------- */ - -void ComputeGyrationMolecule::compute_array() -{ - int i,j; - tagint imol; - double dx,dy,dz,massone; - double unwrap[3]; - - invoked_array = update->ntimestep; - - molcom(); - - for (i = 0; i < nmolecules; i++) - for (j = 0; j < 6; j++) - rgt[i][j] = 0.0; - - double **x = atom->x; - int *mask = atom->mask; - tagint *molecule = atom->molecule; - int *type = atom->type; - imageint *image = atom->image; - double *mass = atom->mass; - double *rmass = atom->rmass; - int nlocal = atom->nlocal; - - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - imol = molecule[i]; - if (molmap) imol = molmap[imol-idlo]; - else imol--; - domain->unmap(x[i],image[i],unwrap); - dx = unwrap[0] - comall[imol][0]; - dy = unwrap[1] - comall[imol][1]; - dz = unwrap[2] - comall[imol][2]; - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - rgt[imol][0] += dx*dx * massone; - rgt[imol][1] += dy*dy * massone; - rgt[imol][2] += dz*dz * massone; - rgt[imol][3] += dx*dy * massone; - rgt[imol][4] += dx*dz * massone; - rgt[imol][5] += dy*dz * massone; - } - - if (nmolecules) - MPI_Allreduce(&rgt[0][0],&array[0][0],nmolecules*6, - MPI_DOUBLE,MPI_SUM,world); - - for (i = 0; i < nmolecules; i++) - for (j = 0; j < 6; j++) - array[i][j] = array[i][j]/masstotal[i]; -} - - -/* ---------------------------------------------------------------------- - calculate per-molecule COM -------------------------------------------------------------------------- */ - -void ComputeGyrationMolecule::molcom() -{ - tagint imol; - double massone; - double unwrap[3]; - - for (int i = 0; i < nmolecules; i++) - com[i][0] = com[i][1] = com[i][2] = 0.0; - - double **x = atom->x; - int *mask = atom->mask; - tagint *molecule = atom->molecule; - int *type = atom->type; - imageint *image = atom->image; - double *mass = atom->mass; - double *rmass = atom->rmass; - int nlocal = atom->nlocal; - - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - imol = molecule[i]; - if (molmap) imol = molmap[imol-idlo]; - else imol--; - domain->unmap(x[i],image[i],unwrap); - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - com[imol][0] += unwrap[0] * massone; - com[imol][1] += unwrap[1] * massone; - com[imol][2] += unwrap[2] * massone; - } - - MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules, - MPI_DOUBLE,MPI_SUM,world); - for (int i = 0; i < nmolecules; i++) { - comall[i][0] /= masstotal[i]; - comall[i][1] /= masstotal[i]; - comall[i][2] /= masstotal[i]; - } -} - -/* ---------------------------------------------------------------------- - memory usage of local data -------------------------------------------------------------------------- */ - -double ComputeGyrationMolecule::memory_usage() -{ - double bytes = (bigint) nmolecules * 2 * sizeof(double); - if (molmap) bytes += (idhi-idlo+1) * sizeof(int); - bytes += (bigint) nmolecules * 2*3 * sizeof(double); - if (tensor) bytes += (bigint) nmolecules * 2*6 * sizeof(double); - else bytes += (bigint) nmolecules * 2 * sizeof(double); - return bytes; -} diff --git a/src/compute_inertia_chunk.cpp b/src/compute_inertia_chunk.cpp new file mode 100644 index 0000000000..f04b291471 --- /dev/null +++ b/src/compute_inertia_chunk.cpp @@ -0,0 +1,256 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "string.h" +#include "compute_inertia_chunk.h" +#include "atom.h" +#include "update.h" +#include "modify.h" +#include "compute_chunk_atom.h" +#include "domain.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputeInertiaChunk::ComputeInertiaChunk(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg != 4) error->all(FLERR,"Illegal compute inertia/chunk command"); + + array_flag = 1; + size_array_cols = 6; + size_array_rows = 0; + size_array_rows_variable = 1; + extarray = 0; + + // ID of compute chunk/atom + + int n = strlen(arg[6]) + 1; + idchunk = new char[n]; + strcpy(idchunk,arg[6]); + + init(); + + // chunk-based data + + nchunk = 1; + maxchunk = 0; + massproc = masstotal = NULL; + com = comall = NULL; + inertia = inertiaall = NULL; + allocate(); +} + +/* ---------------------------------------------------------------------- */ + +ComputeInertiaChunk::~ComputeInertiaChunk() +{ + delete [] idchunk; + memory->destroy(massproc); + memory->destroy(masstotal); + memory->destroy(com); + memory->destroy(comall); + memory->destroy(inertia); + memory->destroy(inertiaall); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeInertiaChunk::init() +{ + int icompute = modify->find_compute(idchunk); + if (icompute < 0) + error->all(FLERR,"Chunk/atom compute does not exist for " + "compute inertia/chunk"); + cchunk = (ComputeChunkAtom *) modify->compute[icompute]; + if (strcmp(cchunk->style,"chunk/atom") != 0) + error->all(FLERR,"Compute inertia/chunk does not use chunk/atom compute"); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeInertiaChunk::compute_array() +{ + int i,j,index; + double dx,dy,dz,massone; + double unwrap[3]; + + invoked_array = update->ntimestep; + + // compute chunk/atom assigns atoms to chunk IDs + // extract ichunk index vector from compute + // ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms + + nchunk = cchunk->setup_chunks(); + cchunk->compute_ichunk(); + int *ichunk = cchunk->ichunk; + + if (nchunk > maxchunk) allocate(); + + // zero local per-chunk values + + for (int i = 0; i < nchunk; i++) { + massproc[i] = 0.0; + com[i][0] = com[i][1] = com[i][2] = 0.0; + for (j = 0; j < 6; j++) inertia[i][j] = 0.0; + } + + // compute COM for each chunk + + double **x = atom->x; + int *mask = atom->mask; + int *type = atom->type; + imageint *image = atom->image; + double *mass = atom->mass; + double *rmass = atom->rmass; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + index = ichunk[i]-1; + if (index < 0) continue; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + domain->unmap(x[i],image[i],unwrap); + massproc[index] += massone; + com[index][0] += unwrap[0] * massone; + com[index][1] += unwrap[1] * massone; + com[index][2] += unwrap[2] * massone; + } + + MPI_Allreduce(massproc,masstotal,nchunk,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&com[0][0],&comall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world); + + for (int i = 0; i < nchunk; i++) { + comall[i][0] /= masstotal[i]; + comall[i][1] /= masstotal[i]; + comall[i][2] /= masstotal[i]; + } + + // compute inertia tensor for each chunk + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + index = ichunk[i]-1; + if (index < 0) continue; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - comall[index][0]; + dy = unwrap[1] - comall[index][1]; + dz = unwrap[2] - comall[index][2]; + inertia[index][0] += massone * (dy*dy + dz*dz); + inertia[index][1] += massone * (dx*dx + dz*dz); + inertia[index][2] += massone * (dx*dx + dy*dy); + inertia[index][3] -= massone * dx*dy; + inertia[index][4] -= massone * dy*dz; + inertia[index][5] -= massone * dx*dz; + } + + MPI_Allreduce(&inertia[0][0],&inertiaall[0][0],6*nchunk, + MPI_DOUBLE,MPI_SUM,world); +} + +/* ---------------------------------------------------------------------- + lock methods: called by fix ave/time + these methods insure vector/array size is locked for Nfreq epoch + by passing lock info along to compute chunk/atom +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + increment lock counter +------------------------------------------------------------------------- */ + +void ComputeInertiaChunk::lock_enable() +{ + cchunk->lockcount++; +} + +/* ---------------------------------------------------------------------- + decrement lock counter in compute chunk/atom, it if still exists +------------------------------------------------------------------------- */ + +void ComputeInertiaChunk::lock_disable() +{ + int icompute = modify->find_compute(idchunk); + if (icompute >= 0) { + cchunk = (ComputeChunkAtom *) modify->compute[icompute]; + cchunk->lockcount--; + } +} + +/* ---------------------------------------------------------------------- + calculate and return # of chunks = length of vector/array +------------------------------------------------------------------------- */ + +int ComputeInertiaChunk::lock_length() +{ + nchunk = cchunk->setup_chunks(); + return nchunk; +} + +/* ---------------------------------------------------------------------- + set the lock from startstep to stopstep +------------------------------------------------------------------------- */ + +void ComputeInertiaChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep) +{ + cchunk->lock(fixptr,startstep,stopstep); +} + +/* ---------------------------------------------------------------------- + unset the lock +------------------------------------------------------------------------- */ + +void ComputeInertiaChunk::unlock(Fix *fixptr) +{ + cchunk->unlock(fixptr); +} + + +/* ---------------------------------------------------------------------- + free and reallocate per-chunk arrays +------------------------------------------------------------------------- */ + +void ComputeInertiaChunk::allocate() +{ + memory->destroy(massproc); + memory->destroy(masstotal); + memory->destroy(com); + memory->destroy(comall); + memory->destroy(inertia); + memory->destroy(inertiaall); + size_array_rows = maxchunk = nchunk; + memory->create(massproc,maxchunk,"inertia/chunk:massproc"); + memory->create(masstotal,maxchunk,"inertia/chunk:masstotal"); + memory->create(com,maxchunk,3,"inertia/chunk:com"); + memory->create(comall,maxchunk,3,"inertia/chunk:comall"); + memory->create(inertia,maxchunk,6,"inertia/chunk:inertia"); + memory->create(inertiaall,maxchunk,6,"inertia/chunk:inertiaall"); + array = inertiaall; +} + +/* ---------------------------------------------------------------------- + memory usage of local data +------------------------------------------------------------------------- */ + +double ComputeInertiaChunk::memory_usage() +{ + double bytes = (bigint) maxchunk * 2 * sizeof(double); + bytes += (bigint) maxchunk * 2*3 * sizeof(double); + bytes += (bigint) maxchunk * 2*6 * sizeof(double); + return bytes; +} diff --git a/src/compute_inertia_molecule.h b/src/compute_inertia_chunk.h similarity index 73% rename from src/compute_inertia_molecule.h rename to src/compute_inertia_chunk.h index 2e3e0c31c0..a42b18149e 100644 --- a/src/compute_inertia_molecule.h +++ b/src/compute_inertia_chunk.h @@ -13,32 +13,42 @@ #ifdef COMPUTE_CLASS -ComputeStyle(inertia/molecule,ComputeInertiaMolecule) +ComputeStyle(inertia/chunk,ComputeInertiaChunk) #else -#ifndef LMP_COMPUTE_INERTIA_MOLECULE_H -#define LMP_COMPUTE_INERTIA_MOLECULE_H +#ifndef LMP_COMPUTE_INERTIA_CHUNK_H +#define LMP_COMPUTE_INERTIA_CHUNK_H #include "compute.h" namespace LAMMPS_NS { -class ComputeInertiaMolecule : public Compute { +class ComputeInertiaChunk : public Compute { public: - ComputeInertiaMolecule(class LAMMPS *, int, char **); - ~ComputeInertiaMolecule(); + ComputeInertiaChunk(class LAMMPS *, int, char **); + ~ComputeInertiaChunk(); void init(); void compute_array(); + + void lock_enable(); + void lock_disable(); + int lock_length(); + void lock(class Fix *, bigint, bigint); + void unlock(class Fix *); + double memory_usage(); private: - int nmolecules; - tagint idlo,idhi; + int nchunk,maxchunk; + char *idchunk; + class ComputeChunkAtom *cchunk; double *massproc,*masstotal; double **com,**comall; double **inertia,**inertiaall; + + void allocate(); }; } diff --git a/src/compute_inertia_molecule.cpp b/src/compute_inertia_molecule.cpp deleted file mode 100644 index eacdf455a1..0000000000 --- a/src/compute_inertia_molecule.cpp +++ /dev/null @@ -1,185 +0,0 @@ -/* ---------------------------------------------------------------------- - 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. -------------------------------------------------------------------------- */ - -#include "compute_inertia_molecule.h" -#include "atom.h" -#include "update.h" -#include "domain.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; - -/* ---------------------------------------------------------------------- */ - -ComputeInertiaMolecule:: -ComputeInertiaMolecule(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg) -{ - if (narg != 3) error->all(FLERR,"Illegal compute inertia/molecule command"); - - if (atom->molecular == 0) - error->all(FLERR,"Compute inertia/molecule requires molecular atom style"); - - array_flag = 1; - size_array_cols = 6; - extarray = 0; - - // setup molecule-based data - - nmolecules = molecules_in_group(idlo,idhi); - size_array_rows = nmolecules; - - memory->create(massproc,nmolecules,"inertia/molecule:massproc"); - memory->create(masstotal,nmolecules,"inertia/molecule:masstotal"); - memory->create(com,nmolecules,3,"inertia/molecule:com"); - memory->create(comall,nmolecules,3,"inertia/molecule:comall"); - memory->create(inertia,nmolecules,6,"inertia/molecule:inertia"); - memory->create(inertiaall,nmolecules,6,"inertia/molecule:inertiaall"); - array = inertiaall; - - // compute masstotal for each molecule - - int *mask = atom->mask; - tagint *molecule = atom->molecule; - int *type = atom->type; - double *mass = atom->mass; - double *rmass = atom->rmass; - int nlocal = atom->nlocal; - - tagint imol; - double massone; - - for (int i = 0; i < nmolecules; i++) massproc[i] = 0.0; - - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - imol = molecule[i]; - if (molmap) imol = molmap[imol-idlo]; - else imol--; - massproc[imol] += massone; - } - - MPI_Allreduce(massproc,masstotal,nmolecules,MPI_DOUBLE,MPI_SUM,world); -} - -/* ---------------------------------------------------------------------- */ - -ComputeInertiaMolecule::~ComputeInertiaMolecule() -{ - memory->destroy(massproc); - memory->destroy(masstotal); - memory->destroy(com); - memory->destroy(comall); - memory->destroy(inertia); - memory->destroy(inertiaall); -} - -/* ---------------------------------------------------------------------- */ - -void ComputeInertiaMolecule::init() -{ - int ntmp = molecules_in_group(idlo,idhi); - if (ntmp != nmolecules) - error->all(FLERR,"Molecule count changed in compute inertia/molecule"); -} - -/* ---------------------------------------------------------------------- */ - -void ComputeInertiaMolecule::compute_array() -{ - int i,j; - tagint imol; - double dx,dy,dz,massone; - double unwrap[3]; - - invoked_array = update->ntimestep; - - double **x = atom->x; - int *mask = atom->mask; - tagint *molecule = atom->molecule; - int *type = atom->type; - imageint *image = atom->image; - double *mass = atom->mass; - double *rmass = atom->rmass; - int nlocal = atom->nlocal; - - // center-of-mass for each molecule - - for (i = 0; i < nmolecules; i++) - com[i][0] = com[i][1] = com[i][2] = 0.0; - - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - imol = molecule[i]; - if (molmap) imol = molmap[imol-idlo]; - else imol--; - domain->unmap(x[i],image[i],unwrap); - com[imol][0] += unwrap[0] * massone; - com[imol][1] += unwrap[1] * massone; - com[imol][2] += unwrap[2] * massone; - } - - MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules, - MPI_DOUBLE,MPI_SUM,world); - for (i = 0; i < nmolecules; i++) { - comall[i][0] /= masstotal[i]; - comall[i][1] /= masstotal[i]; - comall[i][2] /= masstotal[i]; - } - - // inertia tensor for each molecule - - for (i = 0; i < nmolecules; i++) - for (j = 0; j < 6; j++) - inertia[i][j] = 0.0; - - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - imol = molecule[i]; - if (molmap) imol = molmap[imol-idlo]; - else imol--; - domain->unmap(x[i],image[i],unwrap); - dx = unwrap[0] - comall[imol][0]; - dy = unwrap[1] - comall[imol][1]; - dz = unwrap[2] - comall[imol][2]; - inertia[imol][0] += massone * (dy*dy + dz*dz); - inertia[imol][1] += massone * (dx*dx + dz*dz); - inertia[imol][2] += massone * (dx*dx + dy*dy); - inertia[imol][3] -= massone * dx*dy; - inertia[imol][4] -= massone * dy*dz; - inertia[imol][5] -= massone * dx*dz; - } - - MPI_Allreduce(&inertia[0][0],&inertiaall[0][0],6*nmolecules, - MPI_DOUBLE,MPI_SUM,world); -} - -/* ---------------------------------------------------------------------- - memory usage of local data -------------------------------------------------------------------------- */ - -double ComputeInertiaMolecule::memory_usage() -{ - double bytes = (bigint) nmolecules * 2 * sizeof(double); - if (molmap) bytes += (idhi-idlo+1) * sizeof(int); - bytes += (bigint) nmolecules * 2*3 * sizeof(double); - bytes += (bigint) nmolecules * 2*6 * sizeof(double); - return bytes; -} diff --git a/src/compute_msd_chunk.cpp b/src/compute_msd_chunk.cpp new file mode 100644 index 0000000000..180b420b33 --- /dev/null +++ b/src/compute_msd_chunk.cpp @@ -0,0 +1,261 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "string.h" +#include "compute_msd_chunk.h" +#include "atom.h" +#include "update.h" +#include "modify.h" +#include "compute_chunk_atom.h" +#include "domain.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputeMSDChunk::ComputeMSDChunk(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg != 4) error->all(FLERR,"Illegal compute msd/chunk command"); + + array_flag = 1; + size_array_cols = 4; + size_array_rows = 0; + size_array_rows_variable = 1; + extarray = 0; + + // ID of compute chunk/atom + + int n = strlen(arg[6]) + 1; + idchunk = new char[n]; + strcpy(idchunk,arg[6]); + + init(); + + massproc = masstotal = NULL; + com = comall = cominit = NULL; + msd = NULL; + + firstflag = 1; +} + +/* ---------------------------------------------------------------------- */ + +ComputeMSDChunk::~ComputeMSDChunk() +{ + delete [] idchunk; + memory->destroy(massproc); + memory->destroy(masstotal); + memory->destroy(com); + memory->destroy(comall); + memory->destroy(cominit); + memory->destroy(msd); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeMSDChunk::init() +{ + int icompute = modify->find_compute(idchunk); + if (icompute < 0) + error->all(FLERR,"Chunk/atom compute does not exist for compute msd/chunk"); + cchunk = (ComputeChunkAtom *) modify->compute[icompute]; + if (strcmp(cchunk->style,"chunk/atom") != 0) + error->all(FLERR,"Compute msd/chunk does not use chunk/atom compute"); +} + +/* ---------------------------------------------------------------------- + compute initial COM for each chunk + only once on timestep compute is defined, when firstflag = 1 +------------------------------------------------------------------------- */ + +void ComputeMSDChunk::setup() +{ + if (!firstflag) return; + firstflag = 0; + + compute_array(); + for (int i = 0; i < nchunk; i++) { + cominit[i][0] = comall[i][0]; + cominit[i][1] = comall[i][1]; + cominit[i][2] = comall[i][2]; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeMSDChunk::compute_array() +{ + int index; + double massone; + double unwrap[3]; + + invoked_array = update->ntimestep; + + // compute chunk/atom assigns atoms to chunk IDs + // extract ichunk index vector from compute + // ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms + + int n = cchunk->setup_chunks(); + cchunk->compute_ichunk(); + int *ichunk = cchunk->ichunk; + + // first time call, allocate per-chunk arrays + // thereafter, require nchunk remain the same + + if (firstflag) allocate(n); + else if (n != nchunk) + error->all(FLERR,"Compute msd/chunk nchunk is not static"); + + // zero local per-chunk values + + for (int i = 0; i < nchunk; i++) { + massproc[i] = 0.0; + com[i][0] = com[i][1] = com[i][2] = 0.0; + } + + // compute current COM for each chunk + + double **x = atom->x; + int *mask = atom->mask; + int *type = atom->type; + imageint *image = atom->image; + double *mass = atom->mass; + double *rmass = atom->rmass; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + index = ichunk[i]-1; + if (index < 0) continue; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + domain->unmap(x[i],image[i],unwrap); + massproc[index] += massone; + com[index][0] += unwrap[0] * massone; + com[index][1] += unwrap[1] * massone; + com[index][2] += unwrap[2] * massone; + } + + MPI_Allreduce(massproc,masstotal,nchunk,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&com[0][0],&comall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world); + + for (int i = 0; i < nchunk; i++) { + comall[i][0] /= masstotal[i]; + comall[i][1] /= masstotal[i]; + comall[i][2] /= masstotal[i]; + } + + // MSD is difference between current and initial COM + // cominit does not yet exist when called from constructor + + if (firstflag) return; + + double dx,dy,dz; + + for (int i = 0; i < nchunk; i++) { + dx = comall[i][0] - cominit[i][0]; + dy = comall[i][1] - cominit[i][1]; + dz = comall[i][2] - cominit[i][2]; + msd[i][0] = dx*dx; + msd[i][1] = dy*dy; + msd[i][2] = dz*dz; + msd[i][3] = dx*dx + dy*dy + dz*dz; + } +} + +/* ---------------------------------------------------------------------- + lock methods: called by fix ave/time + these methods insure vector/array size is locked for Nfreq epoch + by passing lock info along to compute chunk/atom +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + increment lock counter +------------------------------------------------------------------------- */ + +void ComputeMSDChunk::lock_enable() +{ + cchunk->lockcount++; +} + +/* ---------------------------------------------------------------------- + decrement lock counter in compute chunk/atom, it if still exists +------------------------------------------------------------------------- */ + +void ComputeMSDChunk::lock_disable() +{ + int icompute = modify->find_compute(idchunk); + if (icompute >= 0) { + cchunk = (ComputeChunkAtom *) modify->compute[icompute]; + cchunk->lockcount--; + } +} + +/* ---------------------------------------------------------------------- + calculate and return # of chunks = length of vector/array +------------------------------------------------------------------------- */ + +int ComputeMSDChunk::lock_length() +{ + nchunk = cchunk->setup_chunks(); + return nchunk; +} + +/* ---------------------------------------------------------------------- + set the lock from startstep to stopstep +------------------------------------------------------------------------- */ + +void ComputeMSDChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep) +{ + cchunk->lock(fixptr,startstep,stopstep); +} + +/* ---------------------------------------------------------------------- + unset the lock +------------------------------------------------------------------------- */ + +void ComputeMSDChunk::unlock(Fix *fixptr) +{ + cchunk->unlock(fixptr); +} + +/* ---------------------------------------------------------------------- + one-time allocate of per-chunk arrays +------------------------------------------------------------------------- */ + +void ComputeMSDChunk::allocate(int n) +{ + size_array_rows = nchunk = n; + memory->create(massproc,nchunk,"msd/chunk:massproc"); + memory->create(masstotal,nchunk,"msd/chunk:masstotal"); + memory->create(com,nchunk,3,"msd/chunk:com"); + memory->create(comall,nchunk,3,"msd/chunk:comall"); + memory->create(cominit,nchunk,3,"msd/chunk:cominit"); + memory->create(msd,nchunk,4,"msd/chunk:msd"); + array = msd; +} + +/* ---------------------------------------------------------------------- + memory usage of local data +------------------------------------------------------------------------- */ + +double ComputeMSDChunk::memory_usage() +{ + double bytes = (bigint) nchunk * 2 * sizeof(double); + bytes += (bigint) nchunk * 3*3 * sizeof(double); + bytes += (bigint) nchunk * 4 * sizeof(double); + return bytes; +} diff --git a/src/compute_msd_molecule.h b/src/compute_msd_chunk.h similarity index 68% rename from src/compute_msd_molecule.h rename to src/compute_msd_chunk.h index 454dfe42a8..8562aa976e 100644 --- a/src/compute_msd_molecule.h +++ b/src/compute_msd_chunk.h @@ -13,33 +13,45 @@ #ifdef COMPUTE_CLASS -ComputeStyle(msd/molecule,ComputeMSDMolecule) +ComputeStyle(msd/chunk,ComputeMSDChunk) #else -#ifndef LMP_COMPUTE_MSD_MOLECULE_H -#define LMP_COMPUTE_MSD_MOLECULE_H +#ifndef LMP_COMPUTE_MSD_CHUNK_H +#define LMP_COMPUTE_MSD_CHUNK_H #include "compute.h" namespace LAMMPS_NS { -class ComputeMSDMolecule : public Compute { +class ComputeMSDChunk : public Compute { public: - ComputeMSDMolecule(class LAMMPS *, int, char **); - ~ComputeMSDMolecule(); + ComputeMSDChunk(class LAMMPS *, int, char **); + ~ComputeMSDChunk(); void init(); + void setup(); void compute_array(); + + void lock_enable(); + void lock_disable(); + int lock_length(); + void lock(class Fix *, bigint, bigint); + void unlock(class Fix *); + double memory_usage(); private: - int nmolecules; - tagint idlo,idhi; - int firstflag; + int nchunk; + char *idchunk; + class ComputeChunkAtom *cchunk; double *massproc,*masstotal; double **com,**comall,**cominit; double **msd; + + int firstflag; + + void allocate(int); }; } @@ -55,11 +67,11 @@ Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. -E: Compute msd/molecule requires molecular atom style +E: Compute com/molecule requires molecular atom style Self-explanatory. -E: Molecule count changed in compute msd/molecule +E: Molecule count changed in compute com/molecule Number of molecules must remain constant over time. diff --git a/src/compute_msd_molecule.cpp b/src/compute_msd_molecule.cpp deleted file mode 100644 index 515436ba3a..0000000000 --- a/src/compute_msd_molecule.cpp +++ /dev/null @@ -1,181 +0,0 @@ -/* ---------------------------------------------------------------------- - 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. -------------------------------------------------------------------------- */ - -#include "compute_msd_molecule.h" -#include "atom.h" -#include "update.h" -#include "domain.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; - -/* ---------------------------------------------------------------------- */ - -ComputeMSDMolecule::ComputeMSDMolecule(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg) -{ - if (narg != 3) error->all(FLERR,"Illegal compute msd/molecule command"); - - if (atom->molecular == 0) - error->all(FLERR,"Compute msd/molecule requires molecular atom style"); - - array_flag = 1; - size_array_cols = 4; - extarray = 0; - - // setup molecule-based data and initial COM positions - - nmolecules = molecules_in_group(idlo,idhi); - size_array_rows = nmolecules; - - memory->create(massproc,nmolecules,"msd/molecule:massproc"); - memory->create(masstotal,nmolecules,"msd/molecule:masstotal"); - memory->create(com,nmolecules,3,"msd/molecule:com"); - memory->create(comall,nmolecules,3,"msd/molecule:comall"); - memory->create(cominit,nmolecules,3,"msd/molecule:cominit"); - memory->create(msd,nmolecules,4,"msd/molecule:msd"); - array = msd; - - // compute masstotal for each molecule - - int *mask = atom->mask; - tagint *molecule = atom->molecule; - int *type = atom->type; - double *mass = atom->mass; - double *rmass = atom->rmass; - int nlocal = atom->nlocal; - - tagint imol; - double massone; - - for (int i = 0; i < nmolecules; i++) massproc[i] = 0.0; - - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - imol = molecule[i]; - if (molmap) imol = molmap[imol-idlo]; - else imol--; - massproc[imol] += massone; - } - - MPI_Allreduce(massproc,masstotal,nmolecules,MPI_DOUBLE,MPI_SUM,world); - - // compute initial COM for each molecule - - firstflag = 1; - compute_array(); - for (int i = 0; i < nmolecules; i++) { - cominit[i][0] = comall[i][0]; - cominit[i][1] = comall[i][1]; - cominit[i][2] = comall[i][2]; - } - firstflag = 0; -} - -/* ---------------------------------------------------------------------- */ - -ComputeMSDMolecule::~ComputeMSDMolecule() -{ - memory->destroy(massproc); - memory->destroy(masstotal); - memory->destroy(com); - memory->destroy(comall); - memory->destroy(cominit); - memory->destroy(msd); -} - -/* ---------------------------------------------------------------------- */ - -void ComputeMSDMolecule::init() -{ - int ntmp = molecules_in_group(idlo,idhi); - if (ntmp != nmolecules) - error->all(FLERR,"Molecule count changed in compute msd/molecule"); -} - -/* ---------------------------------------------------------------------- */ - -void ComputeMSDMolecule::compute_array() -{ - tagint imol; - double dx,dy,dz,massone; - double unwrap[3]; - - invoked_array = update->ntimestep; - - // compute current COM positions - - for (int i = 0; i < nmolecules; i++) - com[i][0] = com[i][1] = com[i][2] = 0.0; - - double **x = atom->x; - int *mask = atom->mask; - tagint *molecule = atom->molecule; - int *type = atom->type; - imageint *image = atom->image; - double *mass = atom->mass; - double *rmass = atom->rmass; - int nlocal = atom->nlocal; - - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - imol = molecule[i]; - if (molmap) imol = molmap[imol-idlo]; - else imol--; - domain->unmap(x[i],image[i],unwrap); - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - com[imol][0] += unwrap[0] * massone; - com[imol][1] += unwrap[1] * massone; - com[imol][2] += unwrap[2] * massone; - } - - MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules, - MPI_DOUBLE,MPI_SUM,world); - for (int i = 0; i < nmolecules; i++) { - comall[i][0] /= masstotal[i]; - comall[i][1] /= masstotal[i]; - comall[i][2] /= masstotal[i]; - } - - // MSD is difference between current and initial COM - // cominit does not yet exist when called from constructor - - if (firstflag) return; - - for (int i = 0; i < nmolecules; i++) { - dx = comall[i][0] - cominit[i][0]; - dy = comall[i][1] - cominit[i][1]; - dz = comall[i][2] - cominit[i][2]; - msd[i][0] = dx*dx; - msd[i][1] = dy*dy; - msd[i][2] = dz*dz; - msd[i][3] = dx*dx + dy*dy + dz*dz; - } -} - -/* ---------------------------------------------------------------------- - memory usage of local data -------------------------------------------------------------------------- */ - -double ComputeMSDMolecule::memory_usage() -{ - double bytes = (bigint) nmolecules * 2 * sizeof(double); - if (molmap) bytes += (idhi-idlo+1) * sizeof(int); - bytes += (bigint) nmolecules * 2*3 * sizeof(double); - bytes += (bigint) nmolecules * 4 * sizeof(double); - return bytes; -} diff --git a/src/compute_property_chunk.cpp b/src/compute_property_chunk.cpp new file mode 100644 index 0000000000..9f61657716 --- /dev/null +++ b/src/compute_property_chunk.cpp @@ -0,0 +1,335 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "string.h" +#include "compute_property_chunk.h" +#include "atom.h" +#include "update.h" +#include "modify.h" +#include "compute_chunk_atom.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputePropertyChunk::ComputePropertyChunk(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg < 5) error->all(FLERR,"Illegal compute property/chunk command"); + + // ID of compute chunk/atom + + int n = strlen(arg[3]) + 1; + idchunk = new char[n]; + strcpy(idchunk,arg[3]); + + init(); + + // parse values + + nvalues = narg - 4; + pack_choice = new FnPtrPack[nvalues]; + countflag = 0; + + int i; + for (int iarg = 4; iarg < narg; iarg++) { + i = iarg-4; + + if (strcmp(arg[iarg],"count") == 0) { + pack_choice[i] = &ComputePropertyChunk::pack_count; + countflag = 1; + } else if (strcmp(arg[iarg],"id") == 0) { + if (!cchunk->compress) + error->all(FLERR,"Compute chunk/atom stores no IDs for " + "compute property/chunk"); + pack_choice[i] = &ComputePropertyChunk::pack_id; + } else if (strcmp(arg[iarg],"coord1") == 0) { + if (cchunk->ncoord < 1) + error->all(FLERR,"Compute chunk/atom stores no coord1 for " + "compute property/chunk"); + pack_choice[i] = &ComputePropertyChunk::pack_coord1; + } else if (strcmp(arg[iarg],"coord2") == 0) { + if (cchunk->ncoord < 2) + error->all(FLERR,"Compute chunk/atom stores no coord2 for " + "compute property/chunk"); + pack_choice[i] = &ComputePropertyChunk::pack_coord2; + } else if (strcmp(arg[iarg],"coord3") == 0) { + if (cchunk->ncoord < 3) + error->all(FLERR,"Compute chunk/atom stores no coord3 for " + "compute property/chunk"); + pack_choice[i] = &ComputePropertyChunk::pack_coord3; + } else error->all(FLERR, + "Invalid keyword in compute property/chunk command"); + } + + // initialization + + nchunk = 1; + maxchunk = 0; + vector = NULL; + array = NULL; + count_one = count_all = NULL; + allocate(); + + if (nvalues == 1) { + vector_flag = 1; + size_vector = 0; + size_vector_variable = 1; + extvector = 0; + } else { + array_flag = 1; + size_array_cols = nvalues; + size_array_rows = 0; + size_array_rows_variable = 1; + extarray = 0; + } +} + +/* ---------------------------------------------------------------------- */ + +ComputePropertyChunk::~ComputePropertyChunk() +{ + delete [] idchunk; + delete [] pack_choice; + memory->destroy(vector); + memory->destroy(array); + memory->destroy(count_one); + memory->destroy(count_all); +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyChunk::init() +{ + int icompute = modify->find_compute(idchunk); + if (icompute < 0) + error->all(FLERR,"Chunk/atom compute does not exist for " + "compute property/chunk"); + cchunk = (ComputeChunkAtom *) modify->compute[icompute]; + if (strcmp(cchunk->style,"chunk/atom") != 0) + error->all(FLERR,"Compute property/chunk does not use chunk/atom compute"); +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyChunk::compute_vector() +{ + invoked_vector = update->ntimestep; + + // compute chunk/atom assigns atoms to chunk IDs + // if need count, extract ichunk index vector from compute + // ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms + + nchunk = cchunk->setup_chunks(); + if (nchunk > maxchunk) allocate(); + + if (countflag) { + cchunk->compute_ichunk(); + ichunk = cchunk->ichunk; + } + + // fill vector + + buf = vector; + (this->*pack_choice[0])(0); +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyChunk::compute_array() +{ + invoked_array = update->ntimestep; + + // compute chunk/atom assigns atoms to chunk IDs + // if need count, extract ichunk index vector from compute + // ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms + + nchunk = cchunk->setup_chunks(); + if (nchunk > maxchunk) allocate(); + + if (countflag) { + cchunk->compute_ichunk(); + ichunk = cchunk->ichunk; + } + + // fill array + + if (array) buf = &array[0][0]; + for (int n = 0; n < nvalues; n++) + (this->*pack_choice[n])(n); +} + +/* ---------------------------------------------------------------------- + lock methods: called by fix ave/time + these methods insure vector/array size is locked for Nfreq epoch + by passing lock info along to compute chunk/atom +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + increment lock counter +------------------------------------------------------------------------- */ + +void ComputePropertyChunk::lock_enable() +{ + cchunk->lockcount++; +} + +/* ---------------------------------------------------------------------- + decrement lock counter in compute chunk/atom, it if still exists +------------------------------------------------------------------------- */ + +void ComputePropertyChunk::lock_disable() +{ + int icompute = modify->find_compute(idchunk); + if (icompute >= 0) { + cchunk = (ComputeChunkAtom *) modify->compute[icompute]; + cchunk->lockcount--; + } +} + +/* ---------------------------------------------------------------------- + calculate and return # of chunks = length of vector/array +------------------------------------------------------------------------- */ + +int ComputePropertyChunk::lock_length() +{ + nchunk = cchunk->setup_chunks(); + return nchunk; +} + +/* ---------------------------------------------------------------------- + set the lock from startstep to stopstep +------------------------------------------------------------------------- */ + +void ComputePropertyChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep) +{ + cchunk->lock(fixptr,startstep,stopstep); +} + +/* ---------------------------------------------------------------------- + unset the lock +------------------------------------------------------------------------- */ + +void ComputePropertyChunk::unlock(Fix *fixptr) +{ + cchunk->unlock(fixptr); +} + +/* ---------------------------------------------------------------------- + free and reallocate per-chunk arrays +------------------------------------------------------------------------- */ + +void ComputePropertyChunk::allocate() +{ + memory->destroy(vector); + memory->destroy(array); + memory->destroy(count_one); + memory->destroy(count_all); + size_array_rows = maxchunk = nchunk; + if (nvalues == 1) memory->create(vector,maxchunk,"property/chunk:vector"); + else memory->create(array,maxchunk,nvalues,"property/chunk:array"); + if (countflag) { + memory->create(count_one,maxchunk,"property/chunk:count_one"); + memory->create(count_one,maxchunk,"property/chunk:count_all"); + } +} + +/* ---------------------------------------------------------------------- + memory usage of local data +------------------------------------------------------------------------- */ + +double ComputePropertyChunk::memory_usage() +{ + double bytes = (bigint) nchunk * nvalues * sizeof(double); + if (countflag) bytes += (bigint) nchunk * 2 * sizeof(int); + return bytes; +} + + +/* ---------------------------------------------------------------------- + one method for every keyword compute property/chunk can output + the property is packed into buf starting at n with stride nvalues + customize a new keyword by adding a method +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyChunk::pack_count(int n) +{ + int index; + + for (int m = 0; m < nchunk; m++) count_one[m] = 0; + + int *mask = atom->mask; + int nlocal = atom->nlocal; + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + index = ichunk[i]-1; + if (index < 0) continue; + count_one[index]++; + } + } + + MPI_Allreduce(count_one,count_all,nchunk,MPI_INT,MPI_SUM,world); + + for (int m = 0; m < nchunk; m++) { + buf[n] = count_all[m]; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyChunk::pack_id(int n) +{ + int *origID = cchunk->chunkID; + for (int m = 0; m < nchunk; m++) { + buf[n] = origID[m]; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyChunk::pack_coord1(int n) +{ + double **coord = cchunk->coord; + for (int m = 0; m < nchunk; m++) { + buf[n] = coord[m][0]; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyChunk::pack_coord2(int n) +{ + double **coord = cchunk->coord; + for (int m = 0; m < nchunk; m++) { + buf[n] = coord[m][1]; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyChunk::pack_coord3(int n) +{ + double **coord = cchunk->coord; + for (int m = 0; m < nchunk; m++) { + buf[n] = coord[m][2]; + n += nvalues; + } +} diff --git a/src/compute_property_molecule.h b/src/compute_property_chunk.h similarity index 67% rename from src/compute_property_molecule.h rename to src/compute_property_chunk.h index c810c3231b..21e9078d5d 100644 --- a/src/compute_property_molecule.h +++ b/src/compute_property_chunk.h @@ -13,37 +13,53 @@ #ifdef COMPUTE_CLASS -ComputeStyle(property/molecule,ComputePropertyMolecule) +ComputeStyle(property/chunk,ComputePropertyChunk) #else -#ifndef LMP_COMPUTE_PROPERTY_MOLECULE_H -#define LMP_COMPUTE_PROPERTY_MOLECULE_H +#ifndef LMP_COMPUTE_CHUNK_MOLECULE_H +#define LMP_COMPUTE_CHUNK_MOLECULE_H #include "compute.h" namespace LAMMPS_NS { -class ComputePropertyMolecule : public Compute { +class ComputePropertyChunk : public Compute { public: - ComputePropertyMolecule(class LAMMPS *, int, char **); - ~ComputePropertyMolecule(); + ComputePropertyChunk(class LAMMPS *, int, char **); + ~ComputePropertyChunk(); void init(); void compute_vector(); void compute_array(); + + void lock_enable(); + void lock_disable(); + int lock_length(); + void lock(class Fix *, bigint, bigint); + void unlock(class Fix *); + double memory_usage(); private: - int nvalues,nmolecules; - tagint idlo,idhi; + int nchunk,maxchunk; + char *idchunk; + class ComputeChunkAtom *cchunk; + int *ichunk; + int nvalues,countflag; double *buf; + int *count_one,*count_all; - typedef void (ComputePropertyMolecule::*FnPtrPack)(int); + void allocate(); + + typedef void (ComputePropertyChunk::*FnPtrPack)(int); FnPtrPack *pack_choice; // ptrs to pack functions - void pack_mol(int); void pack_count(int); + void pack_id(int); + void pack_coord1(int); + void pack_coord2(int); + void pack_coord3(int); }; } diff --git a/src/compute_property_molecule.cpp b/src/compute_property_molecule.cpp deleted file mode 100644 index dafe2db680..0000000000 --- a/src/compute_property_molecule.cpp +++ /dev/null @@ -1,173 +0,0 @@ -/* ---------------------------------------------------------------------- - 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. -------------------------------------------------------------------------- */ - -#include "string.h" -#include "compute_property_molecule.h" -#include "atom.h" -#include "update.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; - -/* ---------------------------------------------------------------------- */ - -ComputePropertyMolecule:: -ComputePropertyMolecule(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg) -{ - if (narg < 4) error->all(FLERR,"Illegal compute property/molecule command"); - - if (atom->molecular == 0) - error->all(FLERR,"Compute property/molecule requires molecular atom style"); - - nvalues = narg - 3; - - pack_choice = new FnPtrPack[nvalues]; - - int i; - for (int iarg = 3; iarg < narg; iarg++) { - i = iarg-3; - - if (strcmp(arg[iarg],"mol") == 0) - pack_choice[i] = &ComputePropertyMolecule::pack_mol; - else if (strcmp(arg[iarg],"count") == 0) - pack_choice[i] = &ComputePropertyMolecule::pack_count; - else error->all(FLERR, - "Invalid keyword in compute property/molecule command"); - } - - // setup molecule-based data - - nmolecules = molecules_in_group(idlo,idhi); - - vector = NULL; - array = NULL; - - if (nvalues == 1) { - memory->create(vector,nmolecules,"property/molecule:vector"); - vector_flag = 1; - size_vector = nmolecules; - extvector = 0; - } else { - memory->create(array,nmolecules,nvalues,"property/molecule:array"); - array_flag = 1; - size_array_rows = nmolecules; - size_array_cols = nvalues; - extarray = 0; - } - - // fill vector or array with molecule values - - if (nvalues == 1) { - buf = vector; - (this->*pack_choice[0])(0); - } else { - if (array) buf = &array[0][0]; - for (int n = 0; n < nvalues; n++) - (this->*pack_choice[n])(n); - } -} - -/* ---------------------------------------------------------------------- */ - -ComputePropertyMolecule::~ComputePropertyMolecule() -{ - delete [] pack_choice; - memory->destroy(vector); - memory->destroy(array); -} - -/* ---------------------------------------------------------------------- */ - -void ComputePropertyMolecule::init() -{ - int ntmp = molecules_in_group(idlo,idhi); - if (ntmp != nmolecules) - error->all(FLERR,"Molecule count changed in compute property/molecule"); -} - -/* ---------------------------------------------------------------------- */ - -void ComputePropertyMolecule::compute_vector() -{ - invoked_vector = update->ntimestep; -} - -/* ---------------------------------------------------------------------- */ - -void ComputePropertyMolecule::compute_array() -{ - invoked_array = update->ntimestep; -} - -/* ---------------------------------------------------------------------- - memory usage of local data -------------------------------------------------------------------------- */ - -double ComputePropertyMolecule::memory_usage() -{ - double bytes = (bigint) nmolecules * nvalues * sizeof(double); - if (molmap) bytes += (idhi-idlo+1) * sizeof(int); - return bytes; -} - -/* ---------------------------------------------------------------------- - one method for every keyword compute property/molecule can output - the atom property is packed into buf starting at n with stride nvalues - customize a new keyword by adding a method -------------------------------------------------------------------------- */ - -void ComputePropertyMolecule::pack_mol(int n) -{ - for (tagint m = idlo; m <= idhi; m++) - if (molmap == NULL || molmap[m-idlo] >= 0) { - buf[n] = m; - n += nvalues; - } -} - -/* ---------------------------------------------------------------------- */ - -void ComputePropertyMolecule::pack_count(int n) -{ - int i,m; - tagint imol; - - int *count_one = new int[nmolecules]; - for (m = 0; m < nmolecules; m++) count_one[m] = 0; - - tagint *molecule = atom->molecule; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - imol = molecule[i]; - if (molmap) imol = molmap[imol-idlo]; - else imol--; - count_one[imol]++; - } - - int *count_all = new int[nmolecules]; - MPI_Allreduce(count_one,count_all,nmolecules,MPI_INT,MPI_SUM,world); - - for (m = 0; m < nmolecules; m++) - if (molmap == NULL || molmap[m] >= 0) { - buf[n] = count_all[m]; - n += nvalues; - } - - delete [] count_one; - delete [] count_all; -} diff --git a/src/compute_temp_chunk.cpp b/src/compute_temp_chunk.cpp new file mode 100644 index 0000000000..66b16ad2e4 --- /dev/null +++ b/src/compute_temp_chunk.cpp @@ -0,0 +1,407 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "string.h" +#include "compute_temp_chunk.h" +#include "atom.h" +#include "update.h" +#include "force.h" +#include "modify.h" +#include "compute_chunk_atom.h" +#include "domain.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputeTempChunk::ComputeTempChunk(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg < 4) error->all(FLERR,"Illegal compute temp/chunk command"); + + vector_flag = 1; + size_vector = 1; + size_vector_variable = 1; + extvector = 0; + + // ID of compute chunk/atom + + int n = strlen(arg[3]) + 1; + idchunk = new char[n]; + strcpy(idchunk,arg[3]); + + biasflag = 0; + init(); + + // optional args + + comflag = 0; + biasflag = 0; + id_bias = NULL; + adof = domain->dimension; + cdof = 0.0; + + int iarg = 4; + while (iarg < narg) { + if (strcmp(arg[iarg],"com") == 0) { + if (iarg+2 > narg) + error->all(FLERR,"Illegal compute temp/chunk command"); + if (strcmp(arg[iarg+1],"yes") == 0) comflag = 1; + else if (strcmp(arg[iarg+1],"no") == 0) comflag = 0; + else error->all(FLERR,"Illegal compute temp/chunk command"); + iarg += 2; + } else if (strcmp(arg[iarg],"bias") == 0) { + if (iarg+2 > narg) + error->all(FLERR,"Illegal compute temp/chunk command"); + biasflag = 1; + int n = strlen(arg[iarg+1]) + 1; + id_bias = new char[n]; + strcpy(id_bias,arg[iarg+1]); + iarg += 2; + } else if (strcmp(arg[iarg],"adof") == 0) { + if (iarg+2 > narg) + error->all(FLERR,"Illegal compute temp/chunk command"); + adof = force->numeric(FLERR,arg[iarg+1]); + iarg += 2; + } else if (strcmp(arg[iarg],"cdof") == 0) { + if (iarg+2 > narg) + error->all(FLERR,"Illegal compute temp/chunk command"); + cdof = force->numeric(FLERR,arg[iarg+1]); + iarg += 2; + } else error->all(FLERR,"Illegal compute temp/chunk command"); + } + + if (comflag && biasflag) + error->all(FLERR,"Cannot use both com and bias with compute temp/chunk"); + + // error check on bias compute + + if (biasflag) { + int i = modify->find_compute(id_bias); + if (i < 0) + error->all(FLERR,"Could not find compute ID for temperature bias"); + tbias = modify->compute[i]; + if (tbias->tempflag == 0) + error->all(FLERR,"Bias compute does not calculate temperature"); + if (tbias->tempbias == 0) + error->all(FLERR,"Bias compute does not calculate a velocity bias"); + } + + // chunk-based data + + nchunk = 1; + maxchunk = 0; + ke = keall = NULL; + count = countall = NULL; + massproc = masstotal = NULL; + vcm = vcmall = NULL; + allocate(); +} + +/* ---------------------------------------------------------------------- */ + +ComputeTempChunk::~ComputeTempChunk() +{ + delete [] idchunk; + delete [] id_bias; + memory->destroy(ke); + memory->destroy(keall); + memory->destroy(count); + memory->destroy(countall); + memory->destroy(massproc); + memory->destroy(masstotal); + memory->destroy(vcm); + memory->destroy(vcmall); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeTempChunk::init() +{ + int icompute = modify->find_compute(idchunk); + if (icompute < 0) + error->all(FLERR,"Chunk/atom compute does not exist for " + "compute temp/chunk"); + cchunk = (ComputeChunkAtom *) modify->compute[icompute]; + if (strcmp(cchunk->style,"chunk/atom") != 0) + error->all(FLERR,"Compute temp/chunk does not use chunk/atom compute"); + + if (biasflag) { + int i = modify->find_compute(id_bias); + if (i < 0) + error->all(FLERR,"Could not find compute ID for temperature bias"); + tbias = modify->compute[i]; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeTempChunk::compute_vector() +{ + int index; + + invoked_vector = update->ntimestep; + + // compute chunk/atom assigns atoms to chunk IDs + // extract ichunk index vector from compute + // ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms + + nchunk = cchunk->setup_chunks(); + cchunk->compute_ichunk(); + int *ichunk = cchunk->ichunk; + + if (nchunk > maxchunk) allocate(); + + // calculate COM velocity for each chunk + + if (comflag) vcm_compute(); + + // remove velocity bias + + if (biasflag) { + if (tbias->invoked_scalar != update->ntimestep) tbias->compute_scalar(); + tbias->remove_bias_all(); + } + + // zero local per-chunk values + + for (int i = 0; i < nchunk; i++) { + count[i] = 0; + ke[i] = 0.0; + } + + // compute temperature for each chunk + // option for removing COM velocity + + double **v = atom->v; + double *mass = atom->mass; + double *rmass = atom->rmass; + int *mask = atom->mask; + int *type = atom->type; + int nlocal = atom->nlocal; + + if (!comflag) { + if (rmass) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + index = ichunk[i]-1; + if (index < 0) continue; + ke[index] += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * + rmass[i]; + count[index]++; + } + } else { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + index = ichunk[i]-1; + if (index < 0) continue; + ke[index] += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * + mass[type[i]]; + count[index]++; + } + } + + } else { + double vx,vy,vz; + if (rmass) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + index = ichunk[i]-1; + if (index < 0) continue; + vx = v[i][0] - vcmall[index][0]; + vy = v[i][1] - vcmall[index][1]; + vz = v[i][2] - vcmall[index][2]; + ke[index] += (vx*vx + vy*vy + vz*vz) * rmass[i]; + count[index]++; + } + } else { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + index = ichunk[i]-1; + if (index < 0) continue; + vx = v[i][0] - vcmall[index][0]; + vy = v[i][1] - vcmall[index][1]; + vz = v[i][2] - vcmall[index][2]; + ke[index] += (vx*vx + vy*vy + vz*vz) * mass[type[i]]; + count[index]++; + } + } + } + + MPI_Allreduce(ke,keall,nchunk,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(count,countall,nchunk,MPI_INT,MPI_SUM,world); + + // restore velocity bias + + if (biasflag) tbias->restore_bias_all(); + + // normalize temperatures by per-chunk DOF + + double dof,tfactor; + double mvv2e = force->mvv2e; + double boltz = force->boltz; + + for (int i = 0; i < nchunk; i++) { + dof = cdof + adof*countall[i]; + if (dof > 0.0) tfactor = mvv2e / (dof * boltz); + else tfactor = 0.0; + keall[i] *= tfactor; + } +} + +/* ---------------------------------------------------------------------- + calculate velocity of COM for each chunk +------------------------------------------------------------------------- */ + +void ComputeTempChunk::vcm_compute() +{ + int index; + double massone; + + int *ichunk = cchunk->ichunk; + + for (int i = 0; i < nchunk; i++) { + vcm[i][0] = vcm[i][1] = vcm[i][2] = 0.0; + massproc[i] = 0.0; + } + + double **v = atom->v; + int *mask = atom->mask; + int *type = atom->type; + double *mass = atom->mass; + double *rmass = atom->rmass; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + index = ichunk[i]-1; + if (index < 0) continue; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + vcm[index][0] += v[i][0] * massone; + vcm[index][1] += v[i][1] * massone; + vcm[index][2] += v[i][2] * massone; + massproc[index] += massone; + } + + MPI_Allreduce(&vcm[0][0],&vcmall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(massproc,masstotal,nchunk,MPI_DOUBLE,MPI_SUM,world); + + for (int i = 0; i < nchunk; i++) { + vcmall[i][0] /= masstotal[i]; + vcmall[i][1] /= masstotal[i]; + vcmall[i][2] /= masstotal[i]; + } +} + +/* ---------------------------------------------------------------------- + lock methods: called by fix ave/time + these methods insure vector/array size is locked for Nfreq epoch + by passing lock info along to compute chunk/atom +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + increment lock counter +------------------------------------------------------------------------- */ + +void ComputeTempChunk::lock_enable() +{ + cchunk->lockcount++; +} + +/* ---------------------------------------------------------------------- + decrement lock counter in compute chunk/atom, it if still exists +------------------------------------------------------------------------- */ + +void ComputeTempChunk::lock_disable() +{ + int icompute = modify->find_compute(idchunk); + if (icompute >= 0) { + cchunk = (ComputeChunkAtom *) modify->compute[icompute]; + cchunk->lockcount--; + } +} + +/* ---------------------------------------------------------------------- + calculate and return # of chunks = length of vector/array +------------------------------------------------------------------------- */ + +int ComputeTempChunk::lock_length() +{ + nchunk = cchunk->setup_chunks(); + return nchunk; +} + +/* ---------------------------------------------------------------------- + set the lock from startstep to stopstep +------------------------------------------------------------------------- */ + +void ComputeTempChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep) +{ + cchunk->lock(fixptr,startstep,stopstep); +} + +/* ---------------------------------------------------------------------- + unset the lock +------------------------------------------------------------------------- */ + +void ComputeTempChunk::unlock(Fix *fixptr) +{ + cchunk->unlock(fixptr); +} + +/* ---------------------------------------------------------------------- + free and reallocate per-chunk arrays +------------------------------------------------------------------------- */ + +void ComputeTempChunk::allocate() +{ + memory->destroy(ke); + memory->destroy(keall); + memory->destroy(count); + memory->destroy(countall); + size_vector = maxchunk = nchunk; + memory->create(ke,maxchunk,"temp/chunk:ke"); + memory->create(keall,maxchunk,"temp/chunk:keall"); + memory->create(count,maxchunk,"temp/chunk:count"); + memory->create(countall,maxchunk,"temp/chunk:countall"); + vector = keall; + + if (comflag) { + memory->destroy(massproc); + memory->destroy(masstotal); + memory->destroy(vcm); + memory->destroy(vcmall); + memory->create(massproc,maxchunk,"vcm/chunk:massproc"); + memory->create(masstotal,maxchunk,"vcm/chunk:masstotal"); + memory->create(vcm,maxchunk,3,"vcm/chunk:vcm"); + memory->create(vcmall,maxchunk,3,"vcm/chunk:vcmall"); + } +} + +/* ---------------------------------------------------------------------- + memory usage of local data +------------------------------------------------------------------------- */ + +double ComputeTempChunk::memory_usage() +{ + double bytes = (bigint) maxchunk * 2 * sizeof(double); + bytes = (bigint) maxchunk * 2 * sizeof(int); + if (comflag) { + bytes += (bigint) maxchunk * 2 * sizeof(double); + bytes += (bigint) maxchunk * 2*3 * sizeof(double); + } + return bytes; +} diff --git a/src/compute_temp_chunk.h b/src/compute_temp_chunk.h new file mode 100644 index 0000000000..ddd5e52214 --- /dev/null +++ b/src/compute_temp_chunk.h @@ -0,0 +1,80 @@ +/* -*- c++ -*- ---------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS + +ComputeStyle(temp/chunk,ComputeTempChunk) + +#else + +#ifndef LMP_COMPUTE_TEMP_CHUNK_H +#define LMP_COMPUTE_TEMP_CHUNK_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeTempChunk : public Compute { + public: + ComputeTempChunk(class LAMMPS *, int, char **); + ~ComputeTempChunk(); + void init(); + void compute_vector(); + + void lock_enable(); + void lock_disable(); + int lock_length(); + void lock(class Fix *, bigint, bigint); + void unlock(class Fix *); + + double memory_usage(); + + private: + int nchunk,maxchunk,comflag,biasflag; + char *idchunk; + class ComputeChunkAtom *cchunk; + double adof,cdof; + char *id_bias; + class Compute *tbias; // ptr to additional bias compute + + double *ke,*keall; + int *count,*countall; + double *massproc,*masstotal; + double **vcm,**vcmall; + + void vcm_compute(); + void allocate(); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Compute com/molecule requires molecular atom style + +Self-explanatory. + +E: Molecule count changed in compute com/molecule + +Number of molecules must remain constant over time. + +*/ diff --git a/src/compute_torque_chunk.cpp b/src/compute_torque_chunk.cpp new file mode 100644 index 0000000000..8b29cc690a --- /dev/null +++ b/src/compute_torque_chunk.cpp @@ -0,0 +1,252 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "string.h" +#include "compute_torque_chunk.h" +#include "atom.h" +#include "update.h" +#include "modify.h" +#include "compute_chunk_atom.h" +#include "domain.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputeTorqueChunk::ComputeTorqueChunk(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg != 4) error->all(FLERR,"Illegal compute inertia/chunk command"); + + array_flag = 1; + size_array_cols = 3; + size_array_rows = 0; + size_array_rows_variable = 1; + extarray = 0; + + // ID of compute chunk/atom + + int n = strlen(arg[6]) + 1; + idchunk = new char[n]; + strcpy(idchunk,arg[6]); + + init(); + + // chunk-based data + + nchunk = 1; + maxchunk = 0; + massproc = masstotal = NULL; + com = comall = NULL; + torque = torqueall = NULL; + allocate(); +} + +/* ---------------------------------------------------------------------- */ + +ComputeTorqueChunk::~ComputeTorqueChunk() +{ + delete [] idchunk; + memory->destroy(massproc); + memory->destroy(masstotal); + memory->destroy(com); + memory->destroy(comall); + memory->destroy(torque); + memory->destroy(torqueall); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeTorqueChunk::init() +{ + int icompute = modify->find_compute(idchunk); + if (icompute < 0) + error->all(FLERR,"Chunk/atom compute does not exist for " + "compute torque/chunk"); + cchunk = (ComputeChunkAtom *) modify->compute[icompute]; + if (strcmp(cchunk->style,"chunk/atom") != 0) + error->all(FLERR,"Compute torque/chunk does not use chunk/atom compute"); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeTorqueChunk::compute_array() +{ + int i,j,index; + double dx,dy,dz,massone; + double unwrap[3]; + + invoked_array = update->ntimestep; + + // compute chunk/atom assigns atoms to chunk IDs + // extract ichunk index vector from compute + // ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms + + nchunk = cchunk->setup_chunks(); + cchunk->compute_ichunk(); + int *ichunk = cchunk->ichunk; + + if (nchunk > maxchunk) allocate(); + + // zero local per-chunk values + + for (int i = 0; i < nchunk; i++) { + massproc[i] = 0.0; + com[i][0] = com[i][1] = com[i][2] = 0.0; + torque[i][0] = torque[i][1] = torque[i][2] = 0.0; + } + + // compute COM for each chunk + + double **x = atom->x; + int *mask = atom->mask; + int *type = atom->type; + imageint *image = atom->image; + double *mass = atom->mass; + double *rmass = atom->rmass; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + index = ichunk[i]-1; + if (index < 0) continue; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + domain->unmap(x[i],image[i],unwrap); + massproc[index] += massone; + com[index][0] += unwrap[0] * massone; + com[index][1] += unwrap[1] * massone; + com[index][2] += unwrap[2] * massone; + } + + MPI_Allreduce(massproc,masstotal,nchunk,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&com[0][0],&comall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world); + + for (int i = 0; i < nchunk; i++) { + comall[i][0] /= masstotal[i]; + comall[i][1] /= masstotal[i]; + comall[i][2] /= masstotal[i]; + } + + // compute torque on each chunk + + double **f = atom->f; + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + index = ichunk[i]-1; + if (index < 0) continue; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - comall[index][0]; + dy = unwrap[1] - comall[index][1]; + dz = unwrap[2] - comall[index][2]; + torque[i][0] += dy*f[i][2] - dz*f[i][1]; + torque[i][1] += dz*f[i][0] - dx*f[i][2]; + torque[i][2] += dx*f[i][1] - dy*f[i][0]; + } + + MPI_Allreduce(&torque[0][0],&torqueall[0][0],3*nchunk, + MPI_DOUBLE,MPI_SUM,world); +} + +/* ---------------------------------------------------------------------- + lock methods: called by fix ave/time + these methods insure vector/array size is locked for Nfreq epoch + by passing lock info along to compute chunk/atom +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + increment lock counter +------------------------------------------------------------------------- */ + +void ComputeTorqueChunk::lock_enable() +{ + cchunk->lockcount++; +} + +/* ---------------------------------------------------------------------- + decrement lock counter in compute chunk/atom, it if still exists +------------------------------------------------------------------------- */ + +void ComputeTorqueChunk::lock_disable() +{ + int icompute = modify->find_compute(idchunk); + if (icompute >= 0) { + cchunk = (ComputeChunkAtom *) modify->compute[icompute]; + cchunk->lockcount--; + } +} + +/* ---------------------------------------------------------------------- + calculate and return # of chunks = length of vector/array +------------------------------------------------------------------------- */ + +int ComputeTorqueChunk::lock_length() +{ + nchunk = cchunk->setup_chunks(); + return nchunk; +} + +/* ---------------------------------------------------------------------- + set the lock from startstep to stopstep +------------------------------------------------------------------------- */ + +void ComputeTorqueChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep) +{ + cchunk->lock(fixptr,startstep,stopstep); +} + +/* ---------------------------------------------------------------------- + unset the lock +------------------------------------------------------------------------- */ + +void ComputeTorqueChunk::unlock(Fix *fixptr) +{ + cchunk->unlock(fixptr); +} + +/* ---------------------------------------------------------------------- + free and reallocate per-chunk arrays +------------------------------------------------------------------------- */ + +void ComputeTorqueChunk::allocate() +{ + memory->destroy(massproc); + memory->destroy(masstotal); + memory->destroy(com); + memory->destroy(comall); + memory->destroy(torque); + memory->destroy(torqueall); + size_array_rows = maxchunk = nchunk; + memory->create(massproc,maxchunk,"torque/chunk:massproc"); + memory->create(masstotal,maxchunk,"torque/chunk:masstotal"); + memory->create(com,maxchunk,3,"torque/chunk:com"); + memory->create(comall,maxchunk,3,"torque/chunk:comall"); + memory->create(torque,maxchunk,6,"torque/chunk:torque"); + memory->create(torqueall,maxchunk,6,"torque/chunk:torqueall"); + array = torqueall; +} + +/* ---------------------------------------------------------------------- + memory usage of local data +------------------------------------------------------------------------- */ + +double ComputeTorqueChunk::memory_usage() +{ + double bytes = (bigint) maxchunk * 2 * sizeof(double); + bytes += (bigint) maxchunk * 2*3 * sizeof(double); + bytes += (bigint) maxchunk * 2*3 * sizeof(double); + return bytes; +} diff --git a/src/compute_torque_chunk.h b/src/compute_torque_chunk.h new file mode 100644 index 0000000000..2a684939f1 --- /dev/null +++ b/src/compute_torque_chunk.h @@ -0,0 +1,75 @@ +/* -*- c++ -*- ---------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS + +ComputeStyle(torque/chunk,ComputeTorqueChunk) + +#else + +#ifndef LMP_COMPUTE_TORQUE_CHUNK_H +#define LMP_COMPUTE_TORQUE_CHUNK_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeTorqueChunk : public Compute { + public: + ComputeTorqueChunk(class LAMMPS *, int, char **); + ~ComputeTorqueChunk(); + void init(); + void compute_array(); + + void lock_enable(); + void lock_disable(); + int lock_length(); + void lock(class Fix *, bigint, bigint); + void unlock(class Fix *); + + double memory_usage(); + + private: + int nchunk,maxchunk; + char *idchunk; + class ComputeChunkAtom *cchunk; + + double *massproc,*masstotal; + double **com,**comall; + double **torque,**torqueall; + + void allocate(); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Compute inertia/molecule requires molecular atom style + +Self-explanatory. + +E: Molecule count changed in compute inertia/molecule + +Number of molecules must remain constant over time. + +*/ diff --git a/src/compute_vcm_chunk.cpp b/src/compute_vcm_chunk.cpp new file mode 100644 index 0000000000..68e08c7ab3 --- /dev/null +++ b/src/compute_vcm_chunk.cpp @@ -0,0 +1,237 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "string.h" +#include "compute_vcm_chunk.h" +#include "atom.h" +#include "update.h" +#include "modify.h" +#include "compute_chunk_atom.h" +#include "domain.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +enum{ONCE,NFREQ,EVERY}; + +/* ---------------------------------------------------------------------- */ + +ComputeVCMChunk::ComputeVCMChunk(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg != 4) error->all(FLERR,"Illegal compute vcm/chunk command"); + + array_flag = 1; + size_array_cols = 3; + size_array_rows = 0; + size_array_rows_variable = 1; + extarray = 0; + + // ID of compute chunk/atom + + int n = strlen(arg[6]) + 1; + idchunk = new char[n]; + strcpy(idchunk,arg[6]); + + init(); + + // chunk-based data + + nchunk = 1; + maxchunk = 0; + massproc = masstotal = NULL; + vcm = vcmall = NULL; + allocate(); + + firstflag = massneed = 1; +} + +/* ---------------------------------------------------------------------- */ + +ComputeVCMChunk::~ComputeVCMChunk() +{ + delete [] idchunk; + memory->destroy(massproc); + memory->destroy(masstotal); + memory->destroy(vcm); + memory->destroy(vcmall); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeVCMChunk::init() +{ + int icompute = modify->find_compute(idchunk); + if (icompute < 0) + error->all(FLERR,"Chunk/atom compute does not exist for compute vcm/chunk"); + cchunk = (ComputeChunkAtom *) modify->compute[icompute]; + if (strcmp(cchunk->style,"chunk/atom") != 0) + error->all(FLERR,"Compute vcm/chunk does not use chunk/atom compute"); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeVCMChunk::setup() +{ + // one-time calculation of per-chunk mass + // done in setup, so that ComputeChunkAtom::setup() is already called + + if (firstflag && cchunk->idsflag == ONCE) { + firstflag = massneed = 0; + compute_array(); + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeVCMChunk::compute_array() +{ + int index; + double massone; + + invoked_array = update->ntimestep; + + // compute chunk/atom assigns atoms to chunk IDs + // extract ichunk index vector from compute + // ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms + + nchunk = cchunk->setup_chunks(); + cchunk->compute_ichunk(); + int *ichunk = cchunk->ichunk; + + if (nchunk > maxchunk) allocate(); + + // zero local per-chunk values + + for (int i = 0; i < nchunk; i++) + vcm[i][0] = vcm[i][1] = vcm[i][2] = 0.0; + if (massneed) + for (int i = 0; i < nchunk; i++) massproc[i] = 0.0; + + // compute VCM for each chunk + + double **v = atom->v; + int *mask = atom->mask; + int *type = atom->type; + double *mass = atom->mass; + double *rmass = atom->rmass; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + index = ichunk[i]-1; + if (index < 0) continue; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + vcm[index][0] += v[i][0] * massone; + vcm[index][1] += v[i][1] * massone; + vcm[index][2] += v[i][2] * massone; + if (massneed) massproc[index] += massone; + } + + MPI_Allreduce(&vcm[0][0],&vcmall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world); + if (massneed) + MPI_Allreduce(massproc,masstotal,nchunk,MPI_DOUBLE,MPI_SUM,world); + + for (int i = 0; i < nchunk; i++) { + vcmall[i][0] /= masstotal[i]; + vcmall[i][1] /= masstotal[i]; + vcmall[i][2] /= masstotal[i]; + } +} + +/* ---------------------------------------------------------------------- + lock methods: called by fix ave/time + these methods insure vector/array size is locked for Nfreq epoch + by passing lock info along to compute chunk/atom +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + increment lock counter +------------------------------------------------------------------------- */ + +void ComputeVCMChunk::lock_enable() +{ + cchunk->lockcount++; +} + +/* ---------------------------------------------------------------------- + decrement lock counter in compute chunk/atom, it if still exists +------------------------------------------------------------------------- */ + +void ComputeVCMChunk::lock_disable() +{ + int icompute = modify->find_compute(idchunk); + if (icompute >= 0) { + cchunk = (ComputeChunkAtom *) modify->compute[icompute]; + cchunk->lockcount--; + } +} + +/* ---------------------------------------------------------------------- + calculate and return # of chunks = length of vector/array +------------------------------------------------------------------------- */ + +int ComputeVCMChunk::lock_length() +{ + nchunk = cchunk->setup_chunks(); + return nchunk; +} + +/* ---------------------------------------------------------------------- + set the lock from startstep to stopstep +------------------------------------------------------------------------- */ + +void ComputeVCMChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep) +{ + cchunk->lock(fixptr,startstep,stopstep); +} + +/* ---------------------------------------------------------------------- + unset the lock +------------------------------------------------------------------------- */ + +void ComputeVCMChunk::unlock(Fix *fixptr) +{ + cchunk->unlock(fixptr); +} + +/* ---------------------------------------------------------------------- + free and reallocate per-chunk arrays +------------------------------------------------------------------------- */ + +void ComputeVCMChunk::allocate() +{ + memory->destroy(massproc); + memory->destroy(masstotal); + memory->destroy(vcm); + memory->destroy(vcmall); + size_array_rows = maxchunk = nchunk; + memory->create(massproc,maxchunk,"vcm/chunk:massproc"); + memory->create(masstotal,maxchunk,"vcm/chunk:masstotal"); + memory->create(vcm,maxchunk,3,"vcm/chunk:vcm"); + memory->create(vcmall,maxchunk,3,"vcm/chunk:vcmall"); + array = vcmall; +} + +/* ---------------------------------------------------------------------- + memory usage of local data +------------------------------------------------------------------------- */ + +double ComputeVCMChunk::memory_usage() +{ + double bytes = (bigint) maxchunk * 2 * sizeof(double); + bytes += (bigint) maxchunk * 2*3 * sizeof(double); + return bytes; +} diff --git a/src/compute_vcm_molecule.h b/src/compute_vcm_chunk.h similarity index 72% rename from src/compute_vcm_molecule.h rename to src/compute_vcm_chunk.h index b5bce320ca..dd49ab8dcf 100644 --- a/src/compute_vcm_molecule.h +++ b/src/compute_vcm_chunk.h @@ -13,31 +13,43 @@ #ifdef COMPUTE_CLASS -ComputeStyle(vcm/molecule,ComputeVCMMolecule) +ComputeStyle(vcm/chunk,ComputeVCMChunk) #else -#ifndef LMP_COMPUTE_VCM_MOLECULE_H -#define LMP_COMPUTE_VCM_MOLECULE_H +#ifndef LMP_COMPUTE_VCM_CHUNK_H +#define LMP_COMPUTE_VCM_CHUNK_H #include "compute.h" namespace LAMMPS_NS { -class ComputeVCMMolecule : public Compute { +class ComputeVCMChunk : public Compute { public: - ComputeVCMMolecule(class LAMMPS *, int, char **); - ~ComputeVCMMolecule(); + ComputeVCMChunk(class LAMMPS *, int, char **); + ~ComputeVCMChunk(); void init(); + void setup(); void compute_array(); + + void lock_enable(); + void lock_disable(); + int lock_length(); + void lock(class Fix *, bigint, bigint); + void unlock(class Fix *); + double memory_usage(); private: - int nmolecules; - tagint idlo,idhi; + int nchunk,maxchunk; + int firstflag,massneed; + char *idchunk; + class ComputeChunkAtom *cchunk; double *massproc,*masstotal; double **vcm,**vcmall; + + void allocate(); }; } diff --git a/src/compute_vcm_molecule.cpp b/src/compute_vcm_molecule.cpp deleted file mode 100644 index b61ed75dfc..0000000000 --- a/src/compute_vcm_molecule.cpp +++ /dev/null @@ -1,145 +0,0 @@ -/* ---------------------------------------------------------------------- - 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. -------------------------------------------------------------------------- */ - -#include "compute_vcm_molecule.h" -#include "atom.h" -#include "update.h" -#include "domain.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; - -/* ---------------------------------------------------------------------- */ - -ComputeVCMMolecule::ComputeVCMMolecule(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg) -{ - if (narg != 3) error->all(FLERR,"Illegal compute vcm/molecule command"); - - if (atom->molecular == 0) - error->all(FLERR,"Compute vcm/molecule requires molecular atom style"); - - array_flag = 1; - size_array_cols = 3; - extarray = 0; - - // setup molecule-based data - - nmolecules = molecules_in_group(idlo,idhi); - size_array_rows = nmolecules; - - memory->create(massproc,nmolecules,"vcm/molecule:massproc"); - memory->create(masstotal,nmolecules,"vcm/molecule:masstotal"); - memory->create(vcm,nmolecules,3,"vcm/molecule:vcm"); - memory->create(vcmall,nmolecules,3,"vcm/molecule:vcmall"); - array = vcmall; - - // compute masstotal for each molecule - - int *mask = atom->mask; - tagint *molecule = atom->molecule; - int *type = atom->type; - double *mass = atom->mass; - double *rmass = atom->rmass; - int nlocal = atom->nlocal; - - tagint imol; - double massone; - - for (int i = 0; i < nmolecules; i++) massproc[i] = 0.0; - - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - imol = molecule[i]; - if (molmap) imol = molmap[imol-idlo]; - else imol--; - massproc[imol] += massone; - } - - MPI_Allreduce(massproc,masstotal,nmolecules,MPI_DOUBLE,MPI_SUM,world); -} - -/* ---------------------------------------------------------------------- */ - -ComputeVCMMolecule::~ComputeVCMMolecule() -{ - memory->destroy(massproc); - memory->destroy(masstotal); - memory->destroy(vcm); - memory->destroy(vcmall); -} - -/* ---------------------------------------------------------------------- */ - -void ComputeVCMMolecule::init() -{ - int ntmp = molecules_in_group(idlo,idhi); - if (ntmp != nmolecules) - error->all(FLERR,"Molecule count changed in compute vcm/molecule"); -} - -/* ---------------------------------------------------------------------- */ - -void ComputeVCMMolecule::compute_array() -{ - tagint imol; - double massone; - - invoked_array = update->ntimestep; - - for (int i = 0; i < nmolecules; i++) - vcm[i][0] = vcm[i][1] = vcm[i][2] = 0.0; - - double **v = atom->v; - int *mask = atom->mask; - tagint *molecule = atom->molecule; - int *type = atom->type; - double *mass = atom->mass; - double *rmass = atom->rmass; - int nlocal = atom->nlocal; - - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - imol = molecule[i]; - if (molmap) imol = molmap[imol-idlo]; - else imol--; - vcm[imol][0] += v[i][0] * massone; - vcm[imol][1] += v[i][1] * massone; - vcm[imol][2] += v[i][2] * massone; - } - - MPI_Allreduce(&vcm[0][0],&vcmall[0][0],3*nmolecules, - MPI_DOUBLE,MPI_SUM,world); - for (int i = 0; i < nmolecules; i++) { - vcmall[i][0] /= masstotal[i]; - vcmall[i][1] /= masstotal[i]; - vcmall[i][2] /= masstotal[i]; - } -} - -/* ---------------------------------------------------------------------- - memory usage of local data -------------------------------------------------------------------------- */ - -double ComputeVCMMolecule::memory_usage() -{ - double bytes = (bigint) nmolecules * 2 * sizeof(double); - if (molmap) bytes += (idhi-idlo+1) * sizeof(int); - bytes += (bigint) nmolecules * 2*3 * sizeof(double); - return bytes; -} diff --git a/src/fix_ave_chunk.cpp b/src/fix_ave_chunk.cpp new file mode 100644 index 0000000000..74380d96e3 --- /dev/null +++ b/src/fix_ave_chunk.cpp @@ -0,0 +1,1038 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "string.h" +#include "fix_ave_chunk.h" +#include "atom.h" +#include "update.h" +#include "force.h" +#include "domain.h" +#include "modify.h" +#include "compute.h" +#include "compute_chunk_atom.h" +#include "input.h" +#include "variable.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +enum{V,F,DENSITY_NUMBER,DENSITY_MASS,MASS,TEMPERATURE,COMPUTE,FIX,VARIABLE}; +enum{SAMPLE,ALL}; +enum{ONE,RUNNING,WINDOW}; +enum{NOSCALE,ATOM}; + +#define INVOKED_PERATOM 8 + +/* ---------------------------------------------------------------------- */ + +FixAveChunk::FixAveChunk(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + if (narg < 7) error->all(FLERR,"Illegal fix ave/chunk command"); + + MPI_Comm_rank(world,&me); + + nevery = force->inumeric(FLERR,arg[3]); + nrepeat = force->inumeric(FLERR,arg[4]); + nfreq = force->inumeric(FLERR,arg[5]); + + int n = strlen(arg[6]) + 1; + idchunk = new char[n]; + strcpy(idchunk,arg[6]); + + global_freq = nfreq; + no_change_box = 1; + time_depend = 1; + + // parse values until one isn't recognized + + int iarg = 7; + which = new int[narg-iarg]; + argindex = new int[narg-iarg]; + ids = new char*[narg-iarg]; + value2index = new int[narg-iarg]; + nvalues = 0; + + while (iarg < narg) { + ids[nvalues] = NULL; + + if (strcmp(arg[iarg],"vx") == 0) { + which[nvalues] = V; + argindex[nvalues++] = 0; + } else if (strcmp(arg[iarg],"vy") == 0) { + which[nvalues] = V; + argindex[nvalues++] = 1; + } else if (strcmp(arg[iarg],"vz") == 0) { + which[nvalues] = V; + argindex[nvalues++] = 2; + + } else if (strcmp(arg[iarg],"fx") == 0) { + which[nvalues] = F; + argindex[nvalues++] = 0; + } else if (strcmp(arg[iarg],"fy") == 0) { + which[nvalues] = F; + argindex[nvalues++] = 1; + } else if (strcmp(arg[iarg],"fz") == 0) { + which[nvalues] = F; + argindex[nvalues++] = 2; + + } else if (strcmp(arg[iarg],"density/number") == 0) { + which[nvalues] = DENSITY_NUMBER; + argindex[nvalues++] = 0; + } else if (strcmp(arg[iarg],"density/mass") == 0) { + which[nvalues] = DENSITY_MASS; + argindex[nvalues++] = 0; + } else if (strcmp(arg[iarg],"mass") == 0) { + which[nvalues] = MASS; + argindex[nvalues++] = 0; + } else if (strcmp(arg[iarg],"temp") == 0) { + which[nvalues] = TEMPERATURE; + argindex[nvalues++] = 0; + + } else if (strncmp(arg[iarg],"c_",2) == 0 || + strncmp(arg[iarg],"f_",2) == 0 || + strncmp(arg[iarg],"v_",2) == 0) { + if (arg[iarg][0] == 'c') which[nvalues] = COMPUTE; + else if (arg[iarg][0] == 'f') which[nvalues] = FIX; + else if (arg[iarg][0] == 'v') which[nvalues] = VARIABLE; + + int n = strlen(arg[iarg]); + char *suffix = new char[n]; + strcpy(suffix,&arg[iarg][2]); + + char *ptr = strchr(suffix,'['); + if (ptr) { + if (suffix[strlen(suffix)-1] != ']') + error->all(FLERR,"Illegal fix ave/chunk command"); + argindex[nvalues] = atoi(ptr+1); + *ptr = '\0'; + } else argindex[nvalues] = 0; + + n = strlen(suffix) + 1; + ids[nvalues] = new char[n]; + strcpy(ids[nvalues],suffix); + nvalues++; + delete [] suffix; + + } else break; + + iarg++; + } + + // optional args + + normflag = ALL; + ave = ONE; + scaleflag = ATOM; + fp = NULL; + nwindow = 0; + biasflag = 0; + id_bias = NULL; + adof = domain->dimension; + cdof = 0.0; + overwrite = 0; + char *title1 = NULL; + char *title2 = NULL; + char *title3 = NULL; + + while (iarg < narg) { + if (strcmp(arg[iarg],"norm") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/chunk command"); + if (strcmp(arg[iarg+1],"all") == 0) normflag = ALL; + else if (strcmp(arg[iarg+1],"sample") == 0) normflag = SAMPLE; + else error->all(FLERR,"Illegal fix ave/chunk command"); + iarg += 2; + } else if (strcmp(arg[iarg],"ave") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/chunk command"); + if (strcmp(arg[iarg+1],"one") == 0) ave = ONE; + else if (strcmp(arg[iarg+1],"running") == 0) ave = RUNNING; + else if (strcmp(arg[iarg+1],"window") == 0) ave = WINDOW; + else error->all(FLERR,"Illegal fix ave/chunk command"); + if (ave == WINDOW) { + if (iarg+3 > narg) error->all(FLERR,"Illegal fix ave/chunk command"); + nwindow = force->inumeric(FLERR,arg[iarg+2]); + if (nwindow <= 0) error->all(FLERR,"Illegal fix ave/chunk command"); + } + iarg += 2; + if (ave == WINDOW) iarg++; + } else if (strcmp(arg[iarg],"scale") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/chunk command"); + if (strcmp(arg[iarg+1],"none") == 0) scaleflag = NOSCALE; + else if (strcmp(arg[iarg+1],"atom") == 0) scaleflag = ATOM; + else error->all(FLERR,"Illegal fix ave/chunk command"); + iarg += 2; + + } else if (strcmp(arg[iarg],"bias") == 0) { + if (iarg+2 > narg) + error->all(FLERR,"Illegal fix ave/chunk command"); + biasflag = 1; + int n = strlen(arg[iarg+1]) + 1; + id_bias = new char[n]; + strcpy(id_bias,arg[iarg+1]); + iarg += 2; + } else if (strcmp(arg[iarg],"adof") == 0) { + if (iarg+2 > narg) + error->all(FLERR,"Illegal fix ave/chunk command"); + adof = force->numeric(FLERR,arg[iarg+1]); + iarg += 2; + } else if (strcmp(arg[iarg],"cdof") == 0) { + if (iarg+2 > narg) + error->all(FLERR,"Illegal fix ave/chunk command"); + cdof = force->numeric(FLERR,arg[iarg+1]); + iarg += 2; + + } else if (strcmp(arg[iarg],"file") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/chunk command"); + if (me == 0) { + fp = fopen(arg[iarg+1],"w"); + if (fp == NULL) { + char str[128]; + sprintf(str,"Cannot open fix ave/chunk file %s",arg[iarg+1]); + error->one(FLERR,str); + } + } + iarg += 2; + } else if (strcmp(arg[iarg],"overwrite") == 0) { + overwrite = 1; + iarg += 1; + } else if (strcmp(arg[iarg],"title1") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/chunk command"); + delete [] title1; + int n = strlen(arg[iarg+1]) + 1; + title1 = new char[n]; + strcpy(title1,arg[iarg+1]); + iarg += 2; + } else if (strcmp(arg[iarg],"title2") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/chunk command"); + delete [] title2; + int n = strlen(arg[iarg+1]) + 1; + title2 = new char[n]; + strcpy(title2,arg[iarg+1]); + iarg += 2; + } else if (strcmp(arg[iarg],"title3") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/chunk command"); + delete [] title3; + int n = strlen(arg[iarg+1]) + 1; + title3 = new char[n]; + strcpy(title3,arg[iarg+1]); + iarg += 2; + } else error->all(FLERR,"Illegal fix ave/chunk command"); + } + + // setup and error check + + if (nevery <= 0 || nrepeat <= 0 || nfreq <= 0) + error->all(FLERR,"Illegal fix ave/chunk command"); + if (nfreq % nevery || (nrepeat-1)*nevery >= nfreq) + error->all(FLERR,"Illegal fix ave/chunk command"); + if (ave != RUNNING && overwrite) + error->all(FLERR,"Illegal fix ave/chunk command"); + + if (biasflag) { + int i = modify->find_compute(id_bias); + if (i < 0) + error->all(FLERR,"Could not find compute ID for temperature bias"); + tbias = modify->compute[i]; + if (tbias->tempflag == 0) + error->all(FLERR,"Bias compute does not calculate temperature"); + if (tbias->tempbias == 0) + error->all(FLERR,"Bias compute does not calculate a velocity bias"); + } + + for (int i = 0; i < nvalues; i++) { + if (which[i] == COMPUTE) { + int icompute = modify->find_compute(ids[i]); + if (icompute < 0) + error->all(FLERR,"Compute ID for fix ave/chunk does not exist"); + if (modify->compute[icompute]->peratom_flag == 0) + error->all(FLERR,"Fix ave/chunk compute does not " + "calculate per-atom values"); + if (argindex[i] == 0 && + modify->compute[icompute]->size_peratom_cols != 0) + error->all(FLERR,"Fix ave/chunk compute does not " + "calculate a per-atom vector"); + if (argindex[i] && modify->compute[icompute]->size_peratom_cols == 0) + error->all(FLERR,"Fix ave/chunk compute does not " + "calculate a per-atom array"); + if (argindex[i] && + argindex[i] > modify->compute[icompute]->size_peratom_cols) + error->all(FLERR, + "Fix ave/chunk compute vector is accessed out-of-range"); + + } else if (which[i] == FIX) { + int ifix = modify->find_fix(ids[i]); + if (ifix < 0) + error->all(FLERR,"Fix ID for fix ave/chunk does not exist"); + if (modify->fix[ifix]->peratom_flag == 0) + error->all(FLERR, + "Fix ave/chunk fix does not calculate per-atom values"); + if (argindex[i] && modify->fix[ifix]->size_peratom_cols != 0) + error->all(FLERR, + "Fix ave/chunk fix does not calculate a per-atom vector"); + if (argindex[i] && modify->fix[ifix]->size_peratom_cols == 0) + error->all(FLERR, + "Fix ave/chunk fix does not calculate a per-atom array"); + if (argindex[i] && argindex[i] > modify->fix[ifix]->size_peratom_cols) + error->all(FLERR,"Fix ave/chunk fix vector is accessed out-of-range"); + } else if (which[i] == VARIABLE) { + int ivariable = input->variable->find(ids[i]); + if (ivariable < 0) + error->all(FLERR,"Variable name for fix ave/chunk does not exist"); + if (input->variable->atomstyle(ivariable) == 0) + error->all(FLERR,"Fix ave/chunk variable is not atom-style variable"); + } + } + + // increment lock counter in compute chunk/atom + // only if nrepeat > 1, so that locking spans multiple timesteps + + int icompute = modify->find_compute(idchunk); + if (icompute < 0) + error->all(FLERR,"Chunk/atom compute does not exist for fix ave/chunk"); + cchunk = (ComputeChunkAtom *) modify->compute[icompute]; + if (strcmp(cchunk->style,"chunk/atom") != 0) + error->all(FLERR,"Fix ave/chunk does not use chunk/atom compute"); + if (nrepeat > 1) cchunk->lockcount++; + + // print file comment lines + + if (fp && me == 0) { + if (title1) fprintf(fp,"%s\n",title1); + else fprintf(fp,"# Chunk-averaged data for fix %s and group %s\n", + id,arg[1]); + if (title2) fprintf(fp,"%s\n",title2); + else fprintf(fp,"# Timestep Number-of-chunks Total-count\n"); + if (title3) fprintf(fp,"%s\n",title3); + else { + int compress = cchunk->compress; + int ncoord = cchunk->ncoord; + if (!compress) { + if (ncoord == 0) fprintf(fp,"# Chunk Ncount"); + else if (ncoord == 1) fprintf(fp,"# Chunk Coord1 Ncount"); + else if (ncoord == 2) fprintf(fp,"# Chunk Coord1 Coord2 Ncount"); + else if (ncoord == 3) + fprintf(fp,"# Chunk Coord1 Coord2 Coord3 Ncount"); + } else { + if (ncoord == 0) fprintf(fp,"# Chunk OrigID Ncount"); + else if (ncoord == 1) fprintf(fp,"# Chunk OrigID Coord1 Ncount"); + else if (ncoord == 2) fprintf(fp,"# Chunk OrigID Coord1 Coord2 Ncount"); + else if (ncoord == 3) + fprintf(fp,"# Chunk OrigID Coord1 Coord2 Coord3 Ncount"); + } + for (int i = 0; i < nvalues; i++) fprintf(fp," %s",arg[7+i]); + fprintf(fp,"\n"); + } + filepos = ftell(fp); + } + + delete [] title1; + delete [] title2; + delete [] title3; + + // this fix produces a global array + // size_array_rows is variable and set by allocate() + + int compress = cchunk->compress; + int ncoord = cchunk->ncoord; + colextra = compress + ncoord; + + array_flag = 1; + size_array_cols = colextra + 1 + nvalues; + size_array_rows_variable = 1; + extarray = 0; + + // initializations + + irepeat = 0; + iwindow = window_limit = 0; + normcount = 0; + + maxvar = 0; + varatom = NULL; + + count_one = count_many = count_sum = count_total = NULL; + count_list = NULL; + values_one = values_many = values_sum = values_total = NULL; + values_list = NULL; + + maxchunk = 0; + nchunk = 1; + allocate(); + + // nvalid = next step on which end_of_step does something + // add nvalid to all computes that store invocation times + // since don't know a priori which are invoked by this fix + // once in end_of_step() can set timestep for ones actually invoked + + nvalid = nextvalid(); + modify->addstep_compute_all(nvalid); +} + +/* ---------------------------------------------------------------------- */ + +FixAveChunk::~FixAveChunk() +{ + delete [] which; + delete [] argindex; + for (int i = 0; i < nvalues; i++) delete [] ids[i]; + delete [] ids; + delete [] value2index; + + if (fp && me == 0) fclose(fp); + + memory->destroy(varatom); + + memory->destroy(count_one); + memory->destroy(count_many); + memory->destroy(count_sum); + memory->destroy(count_total); + memory->destroy(count_list); + memory->destroy(values_one); + memory->destroy(values_many); + memory->destroy(values_sum); + memory->destroy(values_total); + memory->destroy(values_list); + + // decrement lock counter in compute chunk/atom, it if still exists + + if (nrepeat > 1) { + int icompute = modify->find_compute(idchunk); + if (icompute >= 0) { + cchunk = (ComputeChunkAtom *) modify->compute[icompute]; + cchunk->lockcount--; + } + } + + delete [] idchunk; +} + +/* ---------------------------------------------------------------------- */ + +int FixAveChunk::setmask() +{ + int mask = 0; + mask |= END_OF_STEP; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixAveChunk::init() +{ + // set indices and check validity of all computes,fixes,variables + // check that fix frequency is acceptable + + int icompute = modify->find_compute(idchunk); + if (icompute < 0) + error->all(FLERR,"Chunk/atom compute does not exist for fix ave/chunk"); + cchunk = (ComputeChunkAtom *) modify->compute[icompute]; + + if (biasflag) { + int i = modify->find_compute(id_bias); + if (i < 0) + error->all(FLERR,"Could not find compute ID for temperature bias"); + tbias = modify->compute[i]; + } + + for (int m = 0; m < nvalues; m++) { + if (which[m] == COMPUTE) { + int icompute = modify->find_compute(ids[m]); + if (icompute < 0) + error->all(FLERR,"Compute ID for fix ave/chunk does not exist"); + value2index[m] = icompute; + + } else if (which[m] == FIX) { + int ifix = modify->find_fix(ids[m]); + if (ifix < 0) + error->all(FLERR,"Fix ID for fix ave/chunk does not exist"); + value2index[m] = ifix; + + if (nevery % modify->fix[ifix]->peratom_freq) + error->all(FLERR, + "Fix for fix ave/chunk not computed at compatible time"); + + } else if (which[m] == VARIABLE) { + int ivariable = input->variable->find(ids[m]); + if (ivariable < 0) + error->all(FLERR,"Variable name for fix ave/chunk does not exist"); + value2index[m] = ivariable; + + } else value2index[m] = -1; + } + + // need to reset nvalid if nvalid < ntimestep b/c minimize was performed + + if (nvalid < update->ntimestep) { + irepeat = 0; + nvalid = nextvalid(); + modify->addstep_compute_all(nvalid); + } +} + +/* ---------------------------------------------------------------------- + only does averaging if nvalid = current timestep + do not call setup_chunks(), even though fix ave/spatial called setup_bins() + b/c could cause nchunk to change if Nfreq epoch crosses 2 runs + does mean that if change_box is used between runs to change box size, + that nchunk may not track it +------------------------------------------------------------------------- */ + +void FixAveChunk::setup(int vflag) +{ + end_of_step(); +} + +/* ---------------------------------------------------------------------- */ + +void FixAveChunk::end_of_step() +{ + int i,j,m,n,index; + + // skip if not step which requires doing something + + bigint ntimestep = update->ntimestep; + if (ntimestep != nvalid) return; + + // first sample within single Nfreq epoch + // zero out arrays that accumulate over many samples, but not across epochs + // invoke setup_chunks() to determine current nchunk + // re-allocate per-chunk arrays if needed + // then invoke lock() so nchunk cannot change until Nfreq epoch is over + // use final arg = -1 for infinite-time window + // wrap setup_chunks in clearstep/addstep b/c it may invoke computes + // both nevery and nfreq are future steps, + // since call below to cchunk->ichunk() + // does not re-invoke internal cchunk compute on this same step + + if (irepeat == 0) { + if (cchunk->computeflag) modify->clearstep_compute(); + nchunk = cchunk->setup_chunks(); + if (cchunk->computeflag) { + modify->addstep_compute(ntimestep+nevery); + modify->addstep_compute(ntimestep+nfreq); + } + allocate(); + if (ave == RUNNING || ave == WINDOW) cchunk->lock(this,ntimestep,-1); + else cchunk->lock(this,ntimestep,ntimestep+(nrepeat-1)*nevery); + for (m = 0; m < nchunk; m++) { + count_many[m] = count_sum[m] = 0.0; + for (i = 0; i < nvalues; i++) values_many[m][i] = 0.0; + } + } + + // zero out arrays for one sample + + for (m = 0; m < nchunk; m++) { + count_one[m] = 0.0; + for (i = 0; i < nvalues; i++) values_one[m][i] = 0.0; + } + + // compute chunk/atom assigns atoms to chunk IDs + // extract ichunk index vector from compute + // ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms + // wrap compute_ichunk in clearstep/addstep b/c it may invoke computes + + if (cchunk->computeflag) modify->clearstep_compute(); + + cchunk->compute_ichunk(); + int *ichunk = cchunk->ichunk; + + if (cchunk->computeflag) modify->addstep_compute(ntimestep+nevery); + + // perform the computation for one sample + // count # of atoms in each bin + // accumulate results of attributes,computes,fixes,variables to local copy + // sum within each chunk, only include atoms in fix group + // compute/fix/variable may invoke computes so wrap with clear/add + + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit && ichunk[i] > 0) + count_one[ichunk[i]-1]++; + + modify->clearstep_compute(); + + for (m = 0; m < nvalues; m++) { + n = value2index[m]; + j = argindex[m]; + + // V,F adds velocities,forces to values + + if (which[m] == V || which[m] == F) { + double **attribute; + if (which[m] == V) attribute = atom->v; + else attribute = atom->f; + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit && ichunk[i] > 0) { + index = ichunk[i]-1; + values_one[index][m] += attribute[i][j]; + } + + // DENSITY_NUMBER adds 1 to values + + } else if (which[m] == DENSITY_NUMBER) { + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit && ichunk[i] > 0) { + index = ichunk[i]-1; + values_one[index][m] += 1.0; + } + + // DENSITY_MASS or MASS adds mass to values + + } else if (which[m] == DENSITY_MASS || which[m] == MASS) { + int *type = atom->type; + double *mass = atom->mass; + double *rmass = atom->rmass; + + if (rmass) { + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit && ichunk[i] > 0) { + index = ichunk[i]-1; + values_one[index][m] += rmass[i]; + } + } else { + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit && ichunk[i] > 0) { + index = ichunk[i]-1; + values_one[index][m] += mass[type[i]]; + } + } + + // TEMPERATURE adds KE to values + // subtract and restore velocity bias if requested + + } else if (which[m] == TEMPERATURE) { + + if (biasflag) { + if (tbias->invoked_scalar != ntimestep) tbias->compute_scalar(); + tbias->remove_bias_all(); + } + + double **v = atom->v; + int *type = atom->type; + double *mass = atom->mass; + double *rmass = atom->rmass; + + if (rmass) { + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit && ichunk[i] > 0) { + index = ichunk[i]-1; + values_one[index][m] += + (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i]; + } + } else { + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit && ichunk[i] > 0) { + index = ichunk[i]-1; + values_one[index][m] += + (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * + mass[type[i]]; + } + } + + if (biasflag) tbias->restore_bias_all(); + + // COMPUTE adds its scalar or vector component to values + // invoke compute if not previously invoked + + } else if (which[m] == COMPUTE) { + Compute *compute = modify->compute[n]; + if (!(compute->invoked_flag & INVOKED_PERATOM)) { + compute->compute_peratom(); + compute->invoked_flag |= INVOKED_PERATOM; + } + double *vector = compute->vector_atom; + double **array = compute->array_atom; + int jm1 = j - 1; + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit && ichunk[i] > 0) { + index = ichunk[i]-1; + if (j == 0) values_one[index][m] += vector[i]; + else values_one[index][m] += array[i][jm1]; + } + + // FIX adds its scalar or vector component to values + // access fix fields, guaranteed to be ready + + } else if (which[m] == FIX) { + double *vector = modify->fix[n]->vector_atom; + double **array = modify->fix[n]->array_atom; + int jm1 = j - 1; + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit && ichunk[i] > 0) { + index = ichunk[i]-1; + if (j == 0) values_one[index][m] += vector[i]; + else values_one[index][m] += array[i][jm1]; + } + + // VARIABLE adds its per-atom quantities to values + // evaluate atom-style variable + + } else if (which[m] == VARIABLE) { + if (nlocal > maxvar) { + maxvar = atom->nmax; + memory->destroy(varatom); + memory->create(varatom,maxvar,"ave/chunk:varatom"); + } + + input->variable->compute_atom(n,igroup,varatom,1,0); + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit && ichunk[i] > 0) { + index = ichunk[i]-1; + values_one[index][m] += varatom[i]; + } + } + } + + // process the current sample + // if normflag = ALL, accumulate values,count separately to many + // if normflag = SAMPLE, one = value/count, accumulate one to many + // count is MPI summed here, value is MPI summed below across samples + // exception is TEMPERATURE: normalize by DOF + // exception is DENSITYs: no normalize by atom count + // exception is scaleflag = NOSCALE : no normalize by atom count + // check last so other options can take precedence + + if (normflag == ALL) { + for (m = 0; m < nchunk; m++) { + count_many[m] += count_one[m]; + for (j = 0; j < nvalues; j++) + values_many[m][j] += values_one[m][j]; + } + } else if (normflag == SAMPLE) { + MPI_Allreduce(count_one,count_many,nchunk,MPI_DOUBLE,MPI_SUM,world); + for (m = 0; m < nchunk; m++) { + if (count_many[m] > 0.0) + for (j = 0; j < nvalues; j++) { + if (which[j] == TEMPERATURE) + values_many[m][j] += values_one[m][j] / (cdof + adof*count_many[m]); + else if (which[j] == DENSITY_NUMBER || which[j] == DENSITY_MASS || + scaleflag == NOSCALE) + values_many[m][j] += values_one[m][j]; + else + values_many[m][j] += values_one[m][j]/count_many[m]; + } + count_sum[m] += count_many[m]; + } + } + + // done if irepeat < nrepeat + // else reset irepeat and nvalid + + irepeat++; + if (irepeat < nrepeat) { + nvalid += nevery; + modify->addstep_compute(nvalid); + return; + } + + irepeat = 0; + nvalid = ntimestep+nfreq - (nrepeat-1)*nevery; + modify->addstep_compute(nvalid); + + // unlock compute chunk/atom at end of Nfreq epoch + // do not unlock if ave = RUNNING or WINDOW + + if (ave == ONE) cchunk->unlock(this); + + // time average across samples + // if normflag = ALL, final is total value / total count + // exception is TEMPERATURE: normalize by DOF for total count + // exception is DENSITYs: normalize by repeat, not total count + // exception is scaleflag == NOSCALE: normalize by repeat, not total count + // check last so other options can take precedence + // if normflag = SAMPLE, final is sum of ave / repeat + + double repeat = nrepeat; + double mv2d = force->mv2d; + + if (normflag == ALL) { + MPI_Allreduce(count_many,count_sum,nchunk,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&values_many[0][0],&values_sum[0][0],nchunk*nvalues, + MPI_DOUBLE,MPI_SUM,world); + for (m = 0; m < nchunk; m++) { + if (count_sum[m] > 0.0) + for (j = 0; j < nvalues; j++) { + if (which[j] == TEMPERATURE) + values_sum[m][j] /= (cdof + adof*count_sum[m]); + else if (which[j] == DENSITY_MASS) + values_sum[m][j] *= mv2d/repeat; + else if (which[j] == DENSITY_NUMBER || scaleflag == NOSCALE) + values_sum[m][j] /= repeat; + else values_sum[m][j] /= count_sum[m]; + } + count_sum[m] /= repeat; + } + } else if (normflag == SAMPLE) { + MPI_Allreduce(&values_many[0][0],&values_sum[0][0],nchunk*nvalues, + MPI_DOUBLE,MPI_SUM,world); + for (m = 0; m < nchunk; m++) { + for (j = 0; j < nvalues; j++) values_sum[m][j] /= repeat; + count_sum[m] /= repeat; + } + } + + // DENSITYs are additionally normalized by chunk volume + // only relevant if chunks are spatial bins + + for (j = 0; j < nvalues; j++) + if (which[j] == DENSITY_NUMBER || which[j] == DENSITY_MASS) { + double chunk_volume = cchunk->chunk_volume_scalar; + for (m = 0; m < nchunk; m++) + values_sum[m][j] /= chunk_volume; + } + + // if ave = ONE, only single Nfreq timestep value is needed + // if ave = RUNNING, combine with all previous Nfreq timestep values + // if ave = WINDOW, comine with nwindow most recent Nfreq timestep values + + if (ave == ONE) { + for (m = 0; m < nchunk; m++) { + for (i = 0; i < nvalues; i++) + values_total[m][i] = values_sum[m][i]; + count_total[m] = count_sum[m]; + } + normcount = 1; + + } else if (ave == RUNNING) { + for (m = 0; m < nchunk; m++) { + for (i = 0; i < nvalues; i++) + values_total[m][i] += values_sum[m][i]; + count_total[m] += count_sum[m]; + } + normcount++; + + } else if (ave == WINDOW) { + for (m = 0; m < nchunk; m++) { + for (i = 0; i < nvalues; i++) { + values_total[m][i] += values_sum[m][i]; + if (window_limit) values_total[m][i] -= values_list[iwindow][m][i]; + values_list[iwindow][m][i] = values_sum[m][i]; + } + count_total[m] += count_sum[m]; + if (window_limit) count_total[m] -= count_list[iwindow][m]; + count_list[iwindow][m] = count_sum[m]; + } + + iwindow++; + if (iwindow == nwindow) { + iwindow = 0; + window_limit = 1; + } + if (window_limit) normcount = nwindow; + else normcount = iwindow; + } + + // output result to file + + if (fp && me == 0) { + if (overwrite) fseek(fp,filepos,SEEK_SET); + double count = 0.0; + for (m = 0; m < nchunk; m++) count += count_total[m]; + fprintf(fp,BIGINT_FORMAT " %d %g\n",ntimestep,nchunk,count); + + int compress = cchunk->compress; + int *chunkID = cchunk->chunkID; + int ncoord = cchunk->ncoord; + double **coord = cchunk->coord; + + if (!compress) { + if (ncoord == 0) { + for (m = 0; m < nchunk; m++) { + fprintf(fp," %d %g",m+1,count_total[m]/normcount); + for (i = 0; i < nvalues; i++) + fprintf(fp," %g",values_total[m][i]/normcount); + fprintf(fp,"\n"); + } + } else if (ncoord == 1) { + for (m = 0; m < nchunk; m++) { + fprintf(fp," %d %g %g",m+1,coord[m][0], + count_total[m]/normcount); + for (i = 0; i < nvalues; i++) + fprintf(fp," %g",values_total[m][i]/normcount); + fprintf(fp,"\n"); + } + } else if (ncoord == 2) { + for (m = 0; m < nchunk; m++) { + fprintf(fp," %d %g %g %g",m+1,coord[m][0],coord[m][1], + count_total[m]/normcount); + for (i = 0; i < nvalues; i++) + fprintf(fp," %g",values_total[m][i]/normcount); + fprintf(fp,"\n"); + } + } else if (ncoord == 3) { + for (m = 0; m < nchunk; m++) { + fprintf(fp," %d %g %g %g %g",m+1, + coord[m][0],coord[m][1],coord[m][2],count_total[m]/normcount); + for (i = 0; i < nvalues; i++) + fprintf(fp," %g",values_total[m][i]/normcount); + fprintf(fp,"\n"); + } + } + } else { + int j; + if (ncoord == 0) { + for (m = 0; m < nchunk; m++) { + fprintf(fp," %d %d %g",m+1,chunkID[m],count_total[m]/normcount); + for (i = 0; i < nvalues; i++) + fprintf(fp," %g",values_total[m][i]/normcount); + fprintf(fp,"\n"); + } + } else if (ncoord == 1) { + for (m = 0; m < nchunk; m++) { + j = chunkID[m]; + fprintf(fp," %d %d %g %g",m+1,j,coord[j-1][0], + count_total[m]/normcount); + for (i = 0; i < nvalues; i++) + fprintf(fp," %g",values_total[m][i]/normcount); + fprintf(fp,"\n"); + } + } else if (ncoord == 2) { + for (m = 0; m < nchunk; m++) { + j = chunkID[m]; + fprintf(fp," %d %d %g %g %g",m+1,j,coord[j-1][0],coord[j-1][1], + count_total[m]/normcount); + for (i = 0; i < nvalues; i++) + fprintf(fp," %g",values_total[m][i]/normcount); + fprintf(fp,"\n"); + } + } else if (ncoord == 3) { + for (m = 0; m < nchunk; m++) { + j = chunkID[m]; + fprintf(fp," %d %d %g %g %g %g",m+1,j,coord[j-1][0], + coord[j-1][1],coord[j-1][2],count_total[m]/normcount); + for (i = 0; i < nvalues; i++) + fprintf(fp," %g",values_total[m][i]/normcount); + fprintf(fp,"\n"); + } + } + } + + fflush(fp); + if (overwrite) { + long fileend = ftell(fp); + ftruncate(fileno(fp),fileend); + } + } +} + +/* ---------------------------------------------------------------------- + allocate all per-chunk vectors +------------------------------------------------------------------------- */ + +void FixAveChunk::allocate() +{ + size_array_rows = nchunk; + + // reallocate chunk arrays if needed + + if (nchunk > maxchunk) { + maxchunk = nchunk; + memory->grow(count_one,nchunk,"ave/chunk:count_one"); + memory->grow(count_many,nchunk,"ave/chunk:count_many"); + memory->grow(count_sum,nchunk,"ave/chunk:count_sum"); + memory->grow(count_total,nchunk,"ave/chunk:count_total"); + + memory->grow(values_one,nchunk,nvalues,"ave/chunk:values_one"); + memory->grow(values_many,nchunk,nvalues,"ave/chunk:values_many"); + memory->grow(values_sum,nchunk,nvalues,"ave/chunk:values_sum"); + memory->grow(values_total,nchunk,nvalues,"ave/chunk:values_total"); + + // only allocate count and values list for ave = WINDOW + + if (ave == WINDOW) { + memory->create(count_list,nwindow,nchunk,"ave/chunk:count_list"); + memory->create(values_list,nwindow,nchunk,nvalues, + "ave/chunk:values_list"); + } + + // reinitialize regrown count/values total since they accumulate + + int i,m; + for (m = 0; m < nchunk; m++) { + for (i = 0; i < nvalues; i++) values_total[m][i] = 0.0; + count_total[m] = 0.0; + } + } +} + +/* ---------------------------------------------------------------------- + return I,J array value + if I exceeds current nchunks, return 0.0 instead of generating an error + columns 1 to colextra = chunkID + ncoord + next column = count, remaining columns = Nvalues +------------------------------------------------------------------------- */ + +double FixAveChunk::compute_array(int i, int j) +{ + if (values_total == NULL) return 0.0; + if (i >= nchunk) return 0.0; + if (j < colextra) { + if (cchunk->compress) { + if (j == 0) return (double) cchunk->chunkID[i]; + return cchunk->coord[i][j-1]; + } else return cchunk->coord[i][j]; + } + j -= colextra + 1; + if (!normcount) return 0.0; + if (j < 0) return count_total[i]/normcount; + return values_total[i][j]/normcount; +} + +/* ---------------------------------------------------------------------- + calculate nvalid = next step on which end_of_step does something + can be this timestep if multiple of nfreq and nrepeat = 1 + else backup from next multiple of nfreq +------------------------------------------------------------------------- */ + +bigint FixAveChunk::nextvalid() +{ + bigint nvalid = (update->ntimestep/nfreq)*nfreq + nfreq; + if (nvalid-nfreq == update->ntimestep && nrepeat == 1) + nvalid = update->ntimestep; + else + nvalid -= (nrepeat-1)*nevery; + if (nvalid < update->ntimestep) nvalid += nfreq; + return nvalid; +} + +/* ---------------------------------------------------------------------- + memory usage of varatom and bins +------------------------------------------------------------------------- */ + +double FixAveChunk::memory_usage() +{ + double bytes = maxvar * sizeof(double); // varatom + bytes += 4*maxchunk * sizeof(double); // count one,many,sum,total + bytes += nvalues*maxchunk * sizeof(double); // values one,many,sum,total + bytes += nwindow*maxchunk * sizeof(double); // count_list + bytes += nwindow*maxchunk*nvalues * sizeof(double); // values_list + return bytes; +} + +/* ---------------------------------------------------------------------- */ + +void FixAveChunk::reset_timestep(bigint ntimestep) +{ + if (ntimestep > nvalid) error->all(FLERR,"Fix ave/chunk missed timestep"); +} diff --git a/src/fix_ave_chunk.h b/src/fix_ave_chunk.h new file mode 100644 index 0000000000..98d74fbea4 --- /dev/null +++ b/src/fix_ave_chunk.h @@ -0,0 +1,173 @@ +/* -*- c++ -*- ---------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(ave/chunk,FixAveChunk) + +#else + +#ifndef LMP_FIX_AVE_CHUNK_H +#define LMP_FIX_AVE_CHUNK_H + +#include "stdio.h" +#include "fix.h" + +namespace LAMMPS_NS { + +class FixAveChunk : public Fix { + public: + FixAveChunk(class LAMMPS *, int, char **); + ~FixAveChunk(); + int setmask(); + void init(); + void setup(int); + void end_of_step(); + double compute_array(int,int); + double memory_usage(); + void reset_timestep(bigint); + + private: + int me,nvalues; + int nrepeat,nfreq,irepeat; + int normflag,scaleflag,overwrite,biasflag,colextra; + bigint nvalid; + double adof,cdof; + char *tstring,*sstring,*id_bias; + int *which,*argindex,*value2index; + char **ids; + class Compute *tbias; // ptr to additional bias compute + FILE *fp; + + int ave,nwindow; + int normcount,iwindow,window_limit; + + int nchunk,maxchunk; + char *idchunk; + class ComputeChunkAtom *cchunk; + + long filepos; + + int maxvar; + double *varatom; + + // one,many,sum vecs/arrays are used with a single Nfreq epoch + // total,list vecs/arrays are used across epochs + + double *count_one,*count_many,*count_sum; + double **values_one,**values_many,**values_sum; + double *count_total,**count_list; + double **values_total,***values_list; + + void allocate(); + bigint nextvalid(); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Cannot use fix ave/spatial z for 2 dimensional model + +Self-explanatory. + +E: Same dimension twice in fix ave/spatial + +Self-explanatory. + +E: Region ID for fix ave/spatial does not exist + +Self-explanatory. + +E: Cannot open fix ave/spatial file %s + +The specified file cannot be opened. Check that the path and name are +correct. + +E: Compute ID for fix ave/spatial does not exist + +Self-explanatory. + +E: Fix ave/spatial compute does not calculate per-atom values + +A compute used by fix ave/spatial must generate per-atom values. + +E: Fix ave/spatial compute does not calculate a per-atom vector + +A compute used by fix ave/spatial must generate per-atom values. + +E: Fix ave/spatial compute does not calculate a per-atom array + +Self-explanatory. + +E: Fix ave/spatial compute vector is accessed out-of-range + +The index for the vector is out of bounds. + +E: Fix ID for fix ave/spatial does not exist + +Self-explanatory. + +E: Fix ave/spatial fix does not calculate per-atom values + +A fix used by fix ave/spatial must generate per-atom values. + +E: Fix ave/spatial fix does not calculate a per-atom vector + +A fix used by fix ave/spatial must generate per-atom values. + +E: Fix ave/spatial fix does not calculate a per-atom array + +Self-explanatory. + +E: Fix ave/spatial fix vector is accessed out-of-range + +The index for the vector is out of bounds. + +E: Variable name for fix ave/spatial does not exist + +Self-explanatory. + +E: Fix ave/spatial variable is not atom-style variable + +A variable used by fix ave/spatial must generate per-atom values. + +E: Fix ave/spatial for triclinic boxes requires units reduced + +Self-explanatory. + +E: Fix ave/spatial settings invalid with changing box size + +If the box size changes, only the units reduced option can be +used. + +E: Fix for fix ave/spatial not computed at compatible time + +Fixes generate their values on specific timesteps. Fix ave/spatial is +requesting a value on a non-allowed timestep. + +E: Fix ave/spatial missed timestep + +You cannot reset the timestep to a value beyond where the fix +expects to next perform averaging. + +*/ diff --git a/src/fix_ave_spatial.cpp b/src/fix_ave_spatial.cpp index 1fcb3fbd17..9d1d46daf5 100644 --- a/src/fix_ave_spatial.cpp +++ b/src/fix_ave_spatial.cpp @@ -29,6 +29,7 @@ #include "compute.h" #include "input.h" #include "variable.h" +#include "comm.h" #include "memory.h" #include "error.h" @@ -49,6 +50,12 @@ enum{NODISCARD,MIXED,YESDISCARD}; FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { + if (comm->me == 0) + error->warning(FLERR,"The fix ave/spatial command has been replaced " + "by the more flexible fix ave/chunk and compute chunk/atom " + "commands -- fix ave/spatial will be removed in " + "the summer of 2015"); + if (narg < 6) error->all(FLERR,"Illegal fix ave/spatial command"); MPI_Comm_rank(world,&me);