git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14428 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/compute_msd_chunk.cpp

@@ -14,10 +14,12 @@
 #include <string.h>
 #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;
 }

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();
 };
 

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;
+  if (strcmp(arg[3],"global") == 0) flavor = GLOBAL;
-  else 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);
+    fwrite(rbuf,sizeof(double),n,fp);
-    else fwrite(&astore[0][0],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[2],n*sizeof(double));
-  if (vecflag) memcpy(vstore,dbuf,n*sizeof(double));
+  else memcpy(&astore[0][0],&dbuf[2],n*sizeof(double));
-  else memcpy(&astore[0][0],dbuf,n*sizeof(double));
 }
 
 /* ----------------------------------------------------------------------

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
+
 };
 
 }