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

src/comm.cpp

@@ -154,7 +154,7 @@ Comm::~Comm()
    map processors to grid, setup xyz split for a uniform grid
 ------------------------------------------------------------------------- */
 
-void Comm::set_proc_grid()
+void Comm::set_proc_grid(int outflag)
 {
   // recv 3d proc grid of another partition if my 3d grid depends on it
 
@@ -237,7 +237,7 @@ void Comm::set_proc_grid()
 
   // print 3d grid info to screen and logfile
 
-  if (me == 0) {
+  if (outflag && me == 0) {
     if (screen) {
       fprintf(screen,"  %d by %d by %d MPI processor grid\n",
               procgrid[0],procgrid[1],procgrid[2]);

src/comm.h

@@ -42,7 +42,7 @@ class Comm : protected Pointers {
   virtual ~Comm();
 
   virtual void init();
-  virtual void set_proc_grid();               // setup 3d grid of procs
+  virtual void set_proc_grid(int outflag = 1); // setup 3d grid of procs
   virtual void setup();                        // setup 3d comm pattern
   virtual void forward_comm(int dummy = 0);    // forward comm of atom coords
   virtual void reverse_comm();                 // reverse comm of forces

src/fix.h

@@ -145,6 +145,7 @@ class Fix : protected Pointers {
   virtual void deform(int) {}
   virtual void reset_target(double) {}
   virtual void reset_dt() {}
+  virtual void reset_timestep(bigint) {}
 
   virtual void read_data_header(char *) {}
   virtual void read_data_section(char *, int, char *) {}

src/fix_ave_atom.cpp

@@ -38,8 +38,6 @@ FixAveAtom::FixAveAtom(LAMMPS *lmp, int narg, char **arg) :
 {
   if (narg < 7) error->all(FLERR,"Illegal fix ave/atom command");
 
-  time_depend = 1;
-
   nevery = atoi(arg[3]);
   nrepeat = atoi(arg[4]);
   peratom_freq = atoi(arg[5]);
@@ -447,3 +445,10 @@ bigint FixAveAtom::nextvalid()
   if (nvalid < update->ntimestep) nvalid += peratom_freq;
   return nvalid;
 }
+
+/* ---------------------------------------------------------------------- */
+
+void FixAveAtom::reset_timestep(bigint ntimestep)
+{
+  if (ntimestep > nvalid) error->all(FLERR,"Fix ave/atom missed timestep");
+}

src/fix_ave_atom.h

@@ -39,6 +39,7 @@ class FixAveAtom : public Fix {
   void copy_arrays(int, int);
   int pack_exchange(int, double *);
   int unpack_exchange(int, double *);
+  void reset_timestep(bigint);
 
  private:
   int nvalues;

src/fix_ave_correlate.cpp

@@ -53,8 +53,6 @@ FixAveCorrelate::FixAveCorrelate(LAMMPS * lmp, int narg, char **arg):
   nfreq = atoi(arg[5]);
 
   global_freq = nfreq;
-  time_depend = 1;
-
   // parse values until one isn't recognized
 
   which = new int[narg-6];
@@ -592,3 +590,10 @@ bigint FixAveCorrelate::nextvalid()
   if (nvalid % nevery) nvalid = (nvalid/nevery)*nevery + nevery;
   return nvalid;
 }
+
+/* ---------------------------------------------------------------------- */
+
+void FixAveCorrelate::reset_timestep(bigint ntimestep)
+{
+  if (ntimestep > nvalid) error->all(FLERR,"Fix ave/correlate missed timestep");
+}

src/fix_ave_correlate.h

@@ -34,6 +34,7 @@ class FixAveCorrelate : public Fix {
   void setup(int);
   void end_of_step();
   double compute_array(int,int);
+  void reset_timestep(bigint);
 
  private:
   int me,nvalues;

src/fix_ave_histo.cpp

@@ -60,7 +60,6 @@ FixAveHisto::FixAveHisto(LAMMPS *lmp, int narg, char **arg) :
   array_flag = 1;
   size_array_cols = 3;
   extarray = 0;
-  time_depend = 1;
 
   lo = atof(arg[6]);
   hi = atof(arg[7]);
@@ -991,3 +990,10 @@ bigint FixAveHisto::nextvalid()
   if (nvalid < update->ntimestep) nvalid += nfreq;
   return nvalid;
 }
+
+/* ---------------------------------------------------------------------- */
+
+void FixAveHisto::reset_timestep(bigint ntimestep)
+{
+  if (ntimestep > nvalid) error->all(FLERR,"Fix ave/histo missed timestep");
+}

src/fix_ave_histo.h

@@ -35,6 +35,7 @@ class FixAveHisto : public Fix {
   void end_of_step();
   double compute_vector(int);
   double compute_array(int,int);
+  void reset_timestep(bigint);
 
  private:
   int me,nvalues;

src/fix_ave_spatial.cpp

@@ -58,7 +58,6 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) :
 
   global_freq = nfreq;
   no_change_box = 1;
-  time_depend = 1;
 
   ndim = 0;
   int iarg = 6;
@@ -1282,3 +1281,10 @@ double FixAveSpatial::memory_usage()
   bytes += nwindow*nbins*nvalues * sizeof(double);  // values_list
   return bytes;
 }
+
+/* ---------------------------------------------------------------------- */
+
+void FixAveSpatial::reset_timestep(bigint ntimestep)
+{
+  if (ntimestep > nvalid) error->all(FLERR,"Fix ave/spatial missed timestep");
+}

src/fix_ave_spatial.h

@@ -35,6 +35,7 @@ class FixAveSpatial : public Fix {
   void end_of_step();
   double compute_array(int,int);
   double memory_usage();
+  void reset_timestep(bigint);
 
  private:
   int me,nvalues;

src/fix_ave_time.cpp

@@ -52,7 +52,6 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
   nfreq = atoi(arg[5]);
 
   global_freq = nfreq;
-  time_depend = 1;
 
   // scan values to count them
   // then read options so know mode = SCALAR/VECTOR before re-reading values
@@ -181,7 +180,8 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
       if (argindex[i] && modify->compute[icompute]->vector_flag == 0)
         error->all(FLERR,"Fix ave/time compute does not calculate a vector");
       if (argindex[i] && argindex[i] > modify->compute[icompute]->size_vector)
-        error->all(FLERR,"Fix ave/time compute vector is accessed out-of-range");
+        error->all(FLERR,
+                   "Fix ave/time compute vector is accessed out-of-range");
 
     } else if (which[i] == COMPUTE && mode == VECTOR) {
       int icompute = modify->find_compute(ids[i]);
@@ -206,7 +206,8 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
       if (argindex[i] && argindex[i] > modify->fix[ifix]->size_vector)
         error->all(FLERR,"Fix ave/time fix vector is accessed out-of-range");
       if (nevery % modify->fix[ifix]->global_freq)
-        error->all(FLERR,"Fix for fix ave/time not computed at compatible time");
+        error->all(FLERR,
+                   "Fix for fix ave/time not computed at compatible time");
 
     } else if (which[i] == FIX && mode == VECTOR) {
       int ifix = modify->find_fix(ids[i]);
@@ -923,3 +924,10 @@ bigint FixAveTime::nextvalid()
   if (nvalid < update->ntimestep) nvalid += nfreq;
   return nvalid;
 }
+
+/* ---------------------------------------------------------------------- */
+
+void FixAveTime::reset_timestep(bigint ntimestep)
+{
+  if (ntimestep > nvalid) error->all(FLERR,"Fix ave/time missed timestep");
+}

src/fix_ave_time.h

@@ -36,6 +36,7 @@ class FixAveTime : public Fix {
   double compute_scalar();
   double compute_vector(int);
   double compute_array(int,int);
+  void reset_timestep(bigint);
 
  private:
   int me,nvalues;

src/output.cpp

@@ -171,12 +171,15 @@ void Output::setup(int headflag, int memflag)
 
   // perform dump at start of run if current timestep is multiple of every
   //   and last dump was not on this timestep
-  // set next_dump to multiple of every
+  // perform dump on first step ever even if not multiple of every
+  //   if firstflag is set
   // do not write on last step of run unless multiple of every
+  // NOTE: what if dump freq is set by variable, then what happens?
+  // set next_dump to multiple of every or variable value
   // set next_dump_any to smallest next_dump
   // if no dumps, set next_dump_any to last+1 so will not influence next
   // wrap dumps that invoke computes with clear/add
-  // if dump not written now, add_all on future step since clear/add is noop
+  // if dump not written now, add_all of next step since clear/add is no-op
 
   int writeflag;
 
@@ -187,6 +190,7 @@ void Output::setup(int headflag, int memflag)
       if (every_dump[idump] && ntimestep % every_dump[idump] == 0 &&
           last_dump[idump] != ntimestep) writeflag = 1;
       if (last_dump[idump] < 0 && dump[idump]->first_flag == 1) writeflag = 1;
+
       if (writeflag) {
         dump[idump]->write();
         last_dump[idump] = ntimestep;
@@ -424,6 +428,79 @@ void Output::write_restart(bigint ntimestep)
   last_restart = ntimestep;
 }
 
+/* ----------------------------------------------------------------------
+   timestep is being changed
+   reset next timestep values for dumps, restart, thermo output
+   reset to smallest value >= new timestep
+   called by update->reset_timestep()
+------------------------------------------------------------------------- */
+
+void Output::reset_timestep(bigint ntimestep)
+{
+  // worry about clear/add of varialbe?
+
+  for (int idump = 0; idump < ndump; idump++) {
+    if (every_dump[idump]) {
+      next_dump[idump] = (ntimestep/every_dump[idump])*every_dump[idump];
+      if (next_dump[idump] < ntimestep) next_dump[idump] += every_dump[idump];
+    } else {
+      /*
+        bigint nextdump = static_cast<bigint>
+        (input->variable->compute_equal(ivar_dump[idump]));
+        if (nextdump <= ntimestep)
+        error->all(FLERR,"Dump every variable returned a bad timestep");
+        next_dump[idump] = nextdump;
+      */
+    }
+    if (idump) next_dump_any = MIN(next_dump_any,next_dump[idump]);
+    else next_dump_any = next_dump[0];
+  }
+
+  if (restart_flag_single) {
+    if (restart_every_single) {
+      next_restart_single =
+        (ntimestep/restart_every_single)*restart_every_single;
+      if (next_restart_single < ntimestep) 
+        next_restart_single += restart_every_single;
+    } else {
+      /*
+        bigint nextrestart = static_cast<bigint>
+        (input->variable->compute_equal(ivar_restart_single));
+        if (nextrestart <= ntimestep)
+        error->all(FLERR,"Restart variable returned a bad timestep");
+        next_restart_single = nextrestart;
+      */
+    }
+  } else next_restart_single = update->laststep + 1;
+
+  if (restart_flag_double) {
+    if (restart_every_double) {
+      next_restart_double =
+        (ntimestep/restart_every_double)*restart_every_double;
+      if (next_restart_double < ntimestep) 
+        next_restart_double += restart_every_double;
+    } else {
+    /*
+      bigint nextrestart = static_cast<bigint>
+        (input->variable->compute_equal(ivar_restart_double));
+        if (nextrestart <= ntimestep)
+        error->all(FLERR,"Restart variable returned a bad timestep");
+        next_restart_double = nextrestart;
+      */
+    }
+  } else next_restart_double = update->laststep + 1;
+
+  next_restart = MIN(next_restart_single,next_restart_double);
+
+  next_thermo = (ntimestep/thermo_every)*thermo_every;
+  if (next_thermo < ntimestep) next_thermo += thermo_every;
+
+  // worry about thermo output on last step?
+
+  next = MIN(next_dump_any,next_restart);
+  next = MIN(next,next_thermo);
+}
+
 /* ----------------------------------------------------------------------
    add a Dump to list of Dumps
 ------------------------------------------------------------------------- */

src/output.h

@@ -64,6 +64,7 @@ class Output : protected Pointers {
   void write(bigint);                // output for current timestep
   void write_dump(bigint);           // force output of dump snapshots
   void write_restart(bigint);        // force output of a restart file
+  void reset_timestep(bigint);       // reset next timestep for all output
 
   void add_dump(int, char **);       // add a Dump to Dump list
   void modify_dump(int, char **);    // modify a Dump

src/read_dump.cpp

@@ -24,6 +24,8 @@
 #include "atom.h"
 #include "atom_vec.h"
 #include "update.h"
+#include "modify.h"
+#include "fix.h"
 #include "domain.h"
 #include "comm.h"
 #include "irregular.h"
@@ -56,7 +58,6 @@ ReadDump::ReadDump(LAMMPS *lmp) : Pointers(lmp)
   fieldtype = NULL;
   fieldlabel = NULL;
   fields = NULL;
-  uflag = ucflag = ucflag_all = NULL;
 
   reader = NULL;
   fp = NULL;
@@ -73,10 +74,6 @@ ReadDump::~ReadDump()
   delete [] fieldtype;
 
   memory->destroy(fields);
-  memory->destroy(uflag);
-  memory->destroy(ucflag);
-  memory->destroy(ucflag_all);
-
   delete reader;
 }
 
@@ -393,13 +390,6 @@ void ReadDump::atoms()
   // uflag[i] = 1 for each owned atom appearing in dump
   // ucflag = similar flag for each chunk atom, used in process_atoms()
 
-  // NOTE: this logic is sloppy
-
-  memory->destroy(uflag);
-  memory->destroy(ucflag);
-  memory->destroy(ucflag_all);
-  uflag = ucflag = ucflag_all = NULL;
-
   int nlocal = atom->nlocal;
   memory->create(uflag,nlocal,"read_dump:uflag");
   for (int i = 0; i < nlocal; i++) uflag[i] = 0;
@@ -439,6 +429,12 @@ void ReadDump::atoms()
     MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_LMP_BIGINT,MPI_SUM,world);
   }
 
+  // can now delete uflag arrays
+
+  memory->destroy(uflag);
+  memory->destroy(ucflag);
+  memory->destroy(ucflag_all);
+
   // delete atom map if created it above
   // else reinitialize map for current atoms
   // do this before migrating atoms to new procs via Irregular
@@ -453,6 +449,7 @@ void ReadDump::atoms()
   }
 
   // overwrite simulation box with dump snapshot box if requested
+  // reallocate processors to box
 
   if (boxflag) {
     domain->boxlo[0] = xlo;
@@ -473,8 +470,7 @@ void ReadDump::atoms()
 
     domain->set_initial_box();
     domain->set_global_box();
-    // would be OK to do this, except it prints out info
+    comm->set_proc_grid(0);
-    //comm->set_proc_grid();
     domain->set_local_box();
   }
 
@@ -714,6 +710,7 @@ void ReadDump::process_atoms(int n)
 
   MPI_Allreduce(ucflag,ucflag_all,n,MPI_INT,MPI_SUM,world);
 
+  int nlocal_previous = atom->nlocal;
   double lamda[3],one[3];
   double *coord;
   
@@ -789,6 +786,17 @@ void ReadDump::process_atoms(int n)
       image[m] = (xbox << 20) | (ybox << 10) | zbox;
     }
   }
+
+  // invoke set_arrays() for fixes that need initialization of new atoms
+  // same as in CreateAtoms
+
+  nlocal = atom->nlocal;
+  for (m = 0; m < modify->nfix; m++) {
+    Fix *fix = modify->fix[m];
+    if (fix->create_attribute)
+      for (i = nlocal_previous; i < nlocal; i++)
+        fix->set_arrays(i);
+  }
 }
 
 /* ----------------------------------------------------------------------

src/rerun.cpp

@@ -57,9 +57,10 @@ void Rerun::command(int narg, char **arg)
   if (nfile == 0 || nfile == narg) error->all(FLERR,"Illegal rerun command");
 
   // parse optional args up until "dump"
+  // user MAXBIGINT -1 so Output can add 1 to it and still be a big int
 
   bigint first = 0;
-  bigint last = MAXBIGINT;
+  bigint last = MAXBIGINT - 1;
   int nevery = 0;
   int nskip = 1;
   int startflag = 0;
@@ -151,9 +152,11 @@ void Rerun::command(int narg, char **arg)
 
     modify->addstep_compute_all(ntimestep);
     update->integrate->setup_minimal(1);
-    output->setup(firstflag,firstflag);
+    modify->end_of_step();
-    firstflag = 0;
+    if (firstflag) output->setup(firstflag,firstflag);
+    else if (output->next) output->write(ntimestep);
 
+    firstflag = 0;
     ntimestep = rd->next(ntimestep,last,nevery,nskip);
     if (ntimestep < 0) break;
   }

src/update.cpp

@@ -356,9 +356,8 @@ void Update::reset_timestep(int narg, char **arg)
 
 /* ----------------------------------------------------------------------
    reset timestep
-   do not allow dump files or a restart to be defined
+   trigger reset of timestep for output and fixes that require it
    do not allow any timestep-dependent fixes to be defined
-   do not allow any dynamic regions to be defined
    reset eflag/vflag global so nothing will think eng/virial are current
    reset invoked flags of computes,
      so nothing will think they are current between runs
@@ -367,23 +366,18 @@ void Update::reset_timestep(int narg, char **arg)
 
 void Update::reset_timestep(bigint newstep)
 {
-  for (int i = 0; i < output->ndump; i++)
+  ntimestep = newstep;
-    if (output->last_dump[i] >= 0)
+  if (ntimestep < 0) error->all(FLERR,"Timestep must be >= 0");
-      error->all(FLERR,
+  if (ntimestep > MAXBIGINT) error->all(FLERR,"Too big a timestep");
-                 "Cannot reset timestep with dump file already written to");
 
-  if (output->restart && output->last_restart >= 0)
+  output->reset_timestep(ntimestep);
-    error->all(FLERR,
-               "Cannot reset timestep with restart file already written");
 
-  for (int i = 0; i < modify->nfix; i++)
+  for (int i = 0; i < modify->nfix; i++) {
     if (modify->fix[i]->time_depend)
       error->all(FLERR,
                  "Cannot reset timestep with a time-dependent fix defined");
-
+    modify->fix[i]->reset_timestep(newstep);
-  for (int i = 0; i < domain->nregion; i++)
+  }
-    if (domain->regions[i]->dynamic_check())
-      error->all(FLERR,"Cannot reset timestep with a dynamic region defined");
 
   eflag_global = vflag_global = -1;
 
@@ -398,9 +392,11 @@ void Update::reset_timestep(bigint newstep)
   for (int i = 0; i < modify->ncompute; i++)
     if (modify->compute[i]->timeflag) modify->compute[i]->clearstep();
 
-  ntimestep = newstep;
+  // NOTE: 7Jun12, adding rerun command, don't think this is required
-  if (ntimestep < 0) error->all(FLERR,"Timestep must be >= 0");
+
-  if (ntimestep > MAXBIGINT) error->all(FLERR,"Too big a timestep");
+  //for (int i = 0; i < domain->nregion; i++)
+  //  if (domain->regions[i]->dynamic_check())
+  //    error->all(FLERR,"Cannot reset timestep with a dynamic region defined");
 }
 
 /* ----------------------------------------------------------------------