diff --git a/src/compute_msd_chunk.cpp b/src/compute_msd_chunk.cpp index 81603f2e44..9d9ddcb3bd 100644 --- a/src/compute_msd_chunk.cpp +++ b/src/compute_msd_chunk.cpp @@ -14,10 +14,12 @@ #include #include "compute_msd_chunk.h" #include "atom.h" +#include "group.h" #include "update.h" #include "modify.h" #include "compute_chunk_atom.h" #include "domain.h" +#include "fix_store.h" #include "memory.h" #include "error.h" @@ -42,25 +44,51 @@ ComputeMSDChunk::ComputeMSDChunk(LAMMPS *lmp, int narg, char **arg) : idchunk = new char[n]; strcpy(idchunk,arg[3]); + firstflag = 1; init(); massproc = masstotal = NULL; - com = comall = cominit = NULL; + com = comall = NULL; msd = NULL; - firstflag = 1; + // create a new fix STORE style for reference positions + // id = compute-ID + COMPUTE_STORE, fix group = compute group + // do not know size of array at this point, just allocate 1x3 array + // fix creation must be done now so that a restart run can + // potentially re-populate the fix array (and change it to correct size) + // otherwise size reset and init will be done in setup() + + n = strlen(id) + strlen("_COMPUTE_STORE") + 1; + id_fix = new char[n]; + strcpy(id_fix,id); + strcat(id_fix,"_COMPUTE_STORE"); + + char **newarg = new char*[6]; + newarg[0] = id_fix; + newarg[1] = group->names[igroup]; + newarg[2] = (char *) "STORE"; + newarg[3] = (char *) "global"; + newarg[4] = (char *) "1"; + newarg[5] = (char *) "1"; + modify->add_fix(6,newarg); + fix = (FixStore *) modify->fix[modify->nfix-1]; + delete [] newarg; } /* ---------------------------------------------------------------------- */ ComputeMSDChunk::~ComputeMSDChunk() { + // check nfix in case all fixes have already been deleted + + if (modify->nfix) modify->delete_fix(id_fix); + + delete [] id_fix; delete [] idchunk; memory->destroy(massproc); memory->destroy(masstotal); memory->destroy(com); memory->destroy(comall); - memory->destroy(cominit); memory->destroy(msd); } @@ -74,6 +102,15 @@ void ComputeMSDChunk::init() cchunk = (ComputeChunkAtom *) modify->compute[icompute]; if (strcmp(cchunk->style,"chunk/atom") != 0) error->all(FLERR,"Compute msd/chunk does not use chunk/atom compute"); + + // set fix which stores reference atom coords + // if firstflag, will be created in setup() + + if (!firstflag) { + int ifix = modify->find_fix(id_fix); + if (ifix < 0) error->all(FLERR,"Could not find compute msd/chunk fix ID"); + fix = (FixStore *) modify->fix[ifix]; + } } /* ---------------------------------------------------------------------- @@ -87,6 +124,13 @@ void ComputeMSDChunk::setup() compute_array(); firstflag = 0; + // if fix->astore is already correct size, restart file set it up + // otherwise reset its size now and initialize to current COM + + if (fix->nrow == nchunk && fix->ncol == 3) return; + fix->reset_global(nchunk,3); + + double **cominit = fix->astore; for (int i = 0; i < nchunk; i++) { cominit[i][0] = comall[i][0]; cominit[i][1] = comall[i][1]; @@ -163,11 +207,12 @@ void ComputeMSDChunk::compute_array() } // MSD is difference between current and initial COM - // cominit does not yet exist when called first time from setup() + // cominit is initilialized by setup() when firstflag is set if (firstflag) return; double dx,dy,dz; + double **cominit = fix->astore; for (int i = 0; i < nchunk; i++) { dx = comall[i][0] - cominit[i][0]; @@ -246,7 +291,6 @@ void ComputeMSDChunk::allocate() 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; } @@ -258,7 +302,7 @@ void ComputeMSDChunk::allocate() double ComputeMSDChunk::memory_usage() { double bytes = (bigint) nchunk * 2 * sizeof(double); - bytes += (bigint) nchunk * 3*3 * sizeof(double); + bytes += (bigint) nchunk * 2*3 * sizeof(double); bytes += (bigint) nchunk * 4 * sizeof(double); return bytes; } diff --git a/src/compute_msd_chunk.h b/src/compute_msd_chunk.h index f220fde233..6706755e16 100644 --- a/src/compute_msd_chunk.h +++ b/src/compute_msd_chunk.h @@ -44,13 +44,14 @@ class ComputeMSDChunk : public Compute { int nchunk; char *idchunk; class ComputeChunkAtom *cchunk; + char *id_fix; + class FixStore *fix; + int firstflag; double *massproc,*masstotal; - double **com,**comall,**cominit; + double **com,**comall; double **msd; - int firstflag; - void allocate(); }; diff --git a/src/fix_store.cpp b/src/fix_store.cpp index ecd052a295..e3de462dd5 100644 --- a/src/fix_store.cpp +++ b/src/fix_store.cpp @@ -32,15 +32,15 @@ FixStore::FixStore(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) if (narg != 6) error->all(FLERR,"Illegal fix store command"); // 4th arg determines GLOBAL vs PERATOM values - // syntax: id group style peratom 0/1 nvalue - // 0/1 flag = not-store or store peratom values in restart file - // nvalue = # of peratom values, N=1 is vector, N>1 is array // syntax: id group style global nrow ncol // Nrow by Ncol array of global values // Ncol=1 is vector, Nrow>1 is array + // syntax: id group style peratom 0/1 nvalue + // 0/1 flag = not-store or store peratom values in restart file + // nvalue = # of peratom values, N=1 is vector, N>1 is array - if (strcmp(arg[3],"peratom") == 0) flavor = PERATOM; - else if (strcmp(arg[3],"global") == 0) flavor = GLOBAL; + if (strcmp(arg[3],"global") == 0) flavor = GLOBAL; + else if (strcmp(arg[3],"peratom") == 0) flavor = PERATOM; else error->all(FLERR,"Invalid fix store command"); // GLOBAL values are always written to restart file @@ -54,7 +54,8 @@ FixStore::FixStore(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) error->all(FLERR,"Invalid fix store command"); vecflag = 0; if (ncol == 1) vecflag = 1; - } else { + } + if (flavor == PERATOM) { restart_peratom = force->inumeric(FLERR,arg[4]); nvalues = force->inumeric(FLERR,arg[5]); if (restart_peratom < 0 or restart_peratom > 1 || nvalues <= 0) @@ -66,21 +67,23 @@ FixStore::FixStore(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) vstore = NULL; astore = NULL; - // allocate vector or array + // allocate vector or array and restart buffer rbuf // for PERATOM, register with Atom class if (flavor == GLOBAL) { if (vecflag) memory->create(vstore,nrow,"fix/store:vstore"); else memory->create(astore,nrow,ncol,"fix/store:astore"); + memory->create(rbuf,nrow*ncol+2,"fix/store:rbuf"); } if (flavor == PERATOM) { grow_arrays(atom->nmax); atom->add_callback(0); if (restart_peratom) atom->add_callback(1); + rbuf = NULL; } // zero the storage - // PERATOM may be exchanged before filled by caller + // PERATOM may be comm->exchanged before filled by caller if (flavor == GLOBAL) { if (vecflag) @@ -99,6 +102,7 @@ FixStore::FixStore(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) for (int j = 0; j < nvalues; j++) astore[i][j] = 0.0; } + } /* ---------------------------------------------------------------------- */ @@ -114,6 +118,7 @@ FixStore::~FixStore() memory->destroy(vstore); memory->destroy(astore); + memory->destroy(rbuf); } /* ---------------------------------------------------------------------- */ @@ -124,18 +129,47 @@ int FixStore::setmask() return mask; } +/* ---------------------------------------------------------------------- + reset size of global vector/array + invoked by caller if size is unknown at time this fix is instantiated + caller will do subsequent initialization +------------------------------------------------------------------------- */ + +void FixStore::reset_global(int nrow_caller, int ncol_caller) +{ + memory->destroy(vstore); + memory->destroy(astore); + memory->destroy(rbuf); + vstore = NULL; + astore = NULL; + + vecflag = 0; + if (ncol_caller == 1) vecflag = 1; + nrow = nrow_caller; + ncol = ncol_caller; + if (vecflag) memory->create(vstore,nrow,"fix/store:vstore"); + else memory->create(astore,nrow,ncol,"fix/store:astore"); + memory->create(rbuf,nrow*ncol+2,"fix/store:rbuf"); +} + /* ---------------------------------------------------------------------- write global array to restart file ------------------------------------------------------------------------- */ void FixStore::write_restart(FILE *fp) { - int n = nrow*ncol; + // fill rbuf with size and vec/array values + + rbuf[0] = nrow; + rbuf[1] = ncol; + if (vecflag) memcpy(&rbuf[2],vstore,nrow*sizeof(double)); + else memcpy(&rbuf[2],&astore[0][0],nrow*ncol*sizeof(double)); + + int n = nrow*ncol + 2; if (comm->me == 0) { int size = n * sizeof(double); fwrite(&size,sizeof(int),1,fp); - if (vecflag) fwrite(vstore,sizeof(double),n,fp); - else fwrite(&astore[0][0],sizeof(double),n,fp); + fwrite(rbuf,sizeof(double),n,fp); } } @@ -145,12 +179,36 @@ void FixStore::write_restart(FILE *fp) void FixStore::restart(char *buf) { - // HOWTO insure size of buf is the same + // first 2 values in buf are vec/array sizes + + double *dbuf = (double *) buf; + int nrow_restart = dbuf[0]; + int ncol_restart = dbuf[1]; + + // if size of vec/array has changed, + // means the restart file is setting size of vec or array and doing init + // because caller did not know size at time this fix was instantiated + // reallocate vstore or astore accordingly + + if (nrow != nrow_restart || ncol != ncol_restart) { + memory->destroy(vstore); + memory->destroy(astore); + memory->destroy(rbuf); + vstore = NULL; + astore = NULL; + + vecflag = 0; + if (ncol_restart == 1) vecflag = 1; + nrow = nrow_restart; + ncol = ncol_restart; + if (vecflag) memory->create(vstore,nrow,"fix/store:vstore"); + else memory->create(astore,nrow,ncol,"fix/store:astore"); + memory->create(rbuf,nrow*ncol+2,"fix/store:rbuf"); + } int n = nrow*ncol; - double *dbuf = (double *) buf; - if (vecflag) memcpy(vstore,dbuf,n*sizeof(double)); - else memcpy(&astore[0][0],dbuf,n*sizeof(double)); + if (vecflag) memcpy(vstore,&dbuf[2],n*sizeof(double)); + else memcpy(&astore[0][0],&dbuf[2],n*sizeof(double)); } /* ---------------------------------------------------------------------- diff --git a/src/fix_store.h b/src/fix_store.h index 8bfef28c90..11bcd936f3 100644 --- a/src/fix_store.h +++ b/src/fix_store.h @@ -31,6 +31,7 @@ class FixStore : public Fix { int nvalues; // number of per-atom values double *vstore; // vector storage for GLOBAL or PERATOM double **astore; // array storage for GLOBAL or PERATOM + double *rbuf; // restart buffer for GLOBAL vec/array FixStore(class LAMMPS *, int, char **); ~FixStore(); @@ -50,9 +51,12 @@ class FixStore : public Fix { double memory_usage(); + void reset_global(int, int); + private: int flavor; // GLOBAL or PERATOM int vecflag; // 1 if ncol=1 or nvalues=1 + }; }