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

src/REPLICA/fix_event.cpp

@@ -46,6 +46,10 @@ FixEvent::FixEvent(LAMMPS *lmp, int narg, char **arg) :
   xold = NULL;
   vold = NULL;
   imageold = NULL;
+  xorig = NULL;
+  vorig = NULL;
+  imageorig = NULL;
+
   grow_arrays(atom->nmax);
   atom->add_callback(0);
 }
@@ -64,6 +68,9 @@ FixEvent::~FixEvent()
   memory->destroy(xold);
   memory->destroy(vold);
   memory->destroy(imageold);
+  memory->destroy(xorig);
+  memory->destroy(vorig);
+  memory->destroy(imageorig);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -121,7 +128,7 @@ void FixEvent::restore_event()
    so can later restore pre-quench state if no event occurs
 ------------------------------------------------------------------------- */
 
-void FixEvent::store_state()
+void FixEvent::store_state_quench()
 {
   double **x = atom->x;
   double **v = atom->v;
@@ -144,7 +151,7 @@ void FixEvent::store_state()
    called after no event detected so can continue
 ------------------------------------------------------------------------- */
 
-void FixEvent::restore_state()
+void FixEvent::restore_state_quench()
 {
   double **x = atom->x;
   double **v = atom->v;
@@ -162,13 +169,57 @@ void FixEvent::restore_state()
   }
 }
 
+/* ----------------------------------------------------------------------
+   store original state of all atoms
+------------------------------------------------------------------------- */
+
+void FixEvent::store_state_dephase()
+{
+  double **x = atom->x;
+  double **v = atom->v;
+  imageint *image = atom->image;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++) {
+    xorig[i][0] = x[i][0];
+    xorig[i][1] = x[i][1];
+    xorig[i][2] = x[i][2];
+    vorig[i][0] = v[i][0];
+    vorig[i][1] = v[i][1];
+    vorig[i][2] = v[i][2];
+    imageorig[i] = image[i];
+  }
+}
+
+/* ----------------------------------------------------------------------
+   restore state of all atoms to original state
+------------------------------------------------------------------------- */
+
+void FixEvent::restore_state_dephase()
+{
+  double **x = atom->x;
+  double **v = atom->v;
+  imageint *image = atom->image;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++) {
+    x[i][0] = xorig[i][0];
+    x[i][1] = xorig[i][1];
+    x[i][2] = xorig[i][2];
+    v[i][0] = vorig[i][0];
+    v[i][1] = vorig[i][1];
+    v[i][2] = vorig[i][2];
+    image[i] = imageorig[i];
+  }
+}
+
 /* ----------------------------------------------------------------------
    memory usage of local atom-based array
 ------------------------------------------------------------------------- */
 
 double FixEvent::memory_usage()
 {
-  double bytes = 6*atom->nmax * sizeof(double);
+  double bytes = 12*atom->nmax * sizeof(double);
   bytes += atom->nmax*sizeof(int);
   return bytes;
 }
@@ -183,6 +234,9 @@ void FixEvent::grow_arrays(int nmax)
   memory->grow(xold,nmax,3,"event:xold");
   memory->grow(vold,nmax,3,"event:vold");
   memory->grow(imageold,nmax,"event:imageold");
+  memory->grow(xorig,nmax,3,"event:xorig");
+  memory->grow(vorig,nmax,3,"event:vorig");
+  memory->grow(imageorig,nmax,"event:imageorig");
 
   // allow compute event to access stored event coords
 
@@ -205,6 +259,13 @@ void FixEvent::copy_arrays(int i, int j, int delflag)
   vold[j][1] = vold[i][1];
   vold[j][2] = vold[i][2];
   imageold[j] = imageold[i];
+  xorig[j][0] = xorig[i][0];
+  xorig[j][1] = xorig[i][1];
+  xorig[j][2] = xorig[i][2];
+  vorig[j][0] = vorig[i][0];
+  vorig[j][1] = vorig[i][1];
+  vorig[j][2] = vorig[i][2];
+  imageorig[j] = imageorig[i];
 }
 
 /* ----------------------------------------------------------------------
@@ -223,8 +284,15 @@ int FixEvent::pack_exchange(int i, double *buf)
   buf[7] = vold[i][1];
   buf[8] = vold[i][2];
   buf[9] = imageold[i];
+  buf[10] = xorig[i][0];
+  buf[11] = xorig[i][1];
+  buf[12] = xorig[i][2];
+  buf[13] = vorig[i][0];
+  buf[14] = vorig[i][1];
+  buf[15] = vorig[i][2];
+  buf[16] = imageorig[i];
 
-  return 10;
+  return 17;
 }
 
 /* ----------------------------------------------------------------------
@@ -242,23 +310,14 @@ int FixEvent::unpack_exchange(int nlocal, double *buf)
   vold[nlocal][0] = buf[6];
   vold[nlocal][1] = buf[7];
   vold[nlocal][2] = buf[8];
-  imageold[nlocal] = static_cast<int>(buf[9]);
+  imageold[nlocal] = static_cast<imageint>(buf[9]);
+  xorig[nlocal][0] = buf[10];
+  xorig[nlocal][1] = buf[11];
+  xorig[nlocal][2] = buf[12];
+  vorig[nlocal][0] = buf[13];
+  vorig[nlocal][1] = buf[14];
+  vorig[nlocal][2] = buf[15];
+  imageorig[nlocal] = static_cast<imageint>(buf[16]);
 
-  return 10;
+  return 17;
-}
-
-/* ----------------------------------------------------------------------
-   pack entire state of Fix into one write
-------------------------------------------------------------------------- */
-
-void FixEvent::write_restart(FILE *fp)
-{
-}
-
-/* ----------------------------------------------------------------------
-   use state info from restart file to restart the Fix
-------------------------------------------------------------------------- */
-
-void FixEvent::restart(char *buf)
-{
 }

src/REPLICA/fix_event.h

@@ -21,7 +21,7 @@ namespace LAMMPS_NS {
 class FixEvent : public Fix {
  public:
   FixEvent(class LAMMPS *, int, char **);
-  virtual ~FixEvent()=0;    // Use destructor to make base class virtual
+  virtual ~FixEvent()=0;      // use destructor to make base class virtual
   int setmask();
 
   double memory_usage();
@@ -29,21 +29,26 @@ class FixEvent : public Fix {
   void copy_arrays(int, int, int);
   int pack_exchange(int, double *);
   int unpack_exchange(int, double *);
-  virtual void write_restart(FILE *);
+  virtual void write_restart(FILE *) {}
-  virtual void restart(char *);
+  virtual void restart(char *) {}
 
   // methods specific to FixEvent
 
-  void store_event();         // store quenched atoms
+  void store_event();              // store quenched atoms
-  void restore_event();       // restore quenched atoms
+  void restore_event();            // restore quenched atoms
-  void store_state();         // store hot atoms
+  void store_state_quench();       // store hot atoms prior to quench
-  void restore_state();       // restore hot atoms
+  void restore_state_quench();     // restore hot atoms after quench
+  void store_state_dephase();      // store atoms before dephase iteration
+  void restore_state_dephase();    // restore atoms if dephase had event
 
  private:
   double **xevent;       // atom coords at last event
   double **xold;         // atom coords for reset/restore
   double **vold;         // atom vels for reset/restore
-  int *imageold;         // image flags for reset/restore
+  imageint *imageold;    // image flags for reset/restore
+  double **xorig;        // original atom coords for reset/restore
+  double **vorig;        // original atom vels for reset/restore
+  imageint *imageorig;   // original image flags for reset/restore
 };
 
 }

src/REPLICA/fix_event_prd.cpp

@@ -91,8 +91,8 @@ void FixEventPRD::restart(char *buf)
   double *list = (double *) buf;
 
   event_number = static_cast<int> (list[n++]);
-  event_timestep = static_cast<int> (list[n++]);
+  event_timestep = static_cast<bigint> (list[n++]);
-  clock = static_cast<int> (list[n++]);
+  clock = static_cast<bigint> (list[n++]);
   replica_number = static_cast<int> (list[n++]);
   correlated_event = static_cast<int> (list[n++]);
   ncoincident = static_cast<int> (list[n++]);

src/REPLICA/fix_event_prd.h

@@ -28,7 +28,7 @@ class FixEventPRD : public FixEvent {
  public:
   int event_number;      // event counter
   bigint event_timestep; // timestep of last event on any replica
-  int clock;             // total elapsed timesteps across all replicas
+  bigint clock;          // total elapsed timesteps across all replicas
   int replica_number;    // replica where last event occured
   int correlated_event;  // 1 if last event was correlated, 0 otherwise
   int ncoincident;       // # of simultaneous events on different replicas

src/REPLICA/prd.cpp

@@ -73,7 +73,9 @@ void PRD::command(int narg, char **arg)
 
   if (narg < 7) error->universe_all(FLERR,"Illegal prd command");
 
-  nsteps = force->inumeric(FLERR,arg[0]);
+  // read as double so can cast to bigint
+
+  int nsteps = force->inumeric(FLERR,arg[0]);
   t_event = force->inumeric(FLERR,arg[1]);
   n_dephase = force->inumeric(FLERR,arg[2]);
   t_dephase = force->inumeric(FLERR,arg[3]);
@@ -132,10 +134,12 @@ void PRD::command(int narg, char **arg)
     memory->create(imageall,natoms,"prd:imageall");
   }
 
-  // random_select = same RNG for each replica for multiple event selection
+  // random_select = same RNG for each replica, for multiple event selection
-  // random_dephase = unique RNG for each replica for dephasing
+  // random_clock = same RNG for each replica, for clock updates
+  // random_dephase = unique RNG for each replica, for dephasing
 
   random_select = new RanPark(lmp,seed);
+  random_clock = new RanPark(lmp,seed+1000);
   random_dephase = new RanMars(lmp,seed+iworld);
 
   // create ComputeTemp class to monitor temperature
@@ -264,10 +268,10 @@ void PRD::command(int narg, char **arg)
   // need this line if quench() does only setup_minimal()
   // update->minimize->setup();
 
-  fix_event->store_state();
+  fix_event->store_state_quench();
   quench();
   ncoincident = 0;
-  share_event(0,0);
+  share_event(0,0,0);
 
   timer->init();
   timer->barrier_start(TIME_LOOP);
@@ -295,29 +299,49 @@ void PRD::command(int narg, char **arg)
 
   nbuild = ndanger = 0;
   time_dephase = time_dynamics = time_quench = time_comm = time_output = 0.0;
+  bigint clock = 0;
 
   timer->barrier_start(TIME_LOOP);
   time_start = timer->array[TIME_LOOP];
 
-  while (update->ntimestep < update->endstep) {
+  int istep = 0;
+
+  while (istep < nsteps) {
     dephase();
 
+    if (stepmode == 0) istep = update->ntimestep - update->beginstep;
+    else istep = clock;
+
     ireplica = -1;
-    while (update->ntimestep < update->endstep) {
+    while (istep < nsteps) {
-      dynamics();
+      dynamics(t_event,time_dynamics);
-      fix_event->store_state();
+      fix_event->store_state_quench();
       quench();
+      clock = clock + t_event*universe->nworlds;
       ireplica = check_event();
       if (ireplica >= 0) break;
-      fix_event->restore_state();
+      fix_event->restore_state_quench();
+      if (stepmode == 0) istep = update->ntimestep - update->beginstep;
+      else istep = clock;
     }
     if (ireplica < 0) break;
 
-    // potentially more efficient for correlated events if don't
+    // decrement clock by random time at which 1 or more events occurred
-    // share until correlated check has completed
+
+    int frac_t_event = t_event;
+    for (int i = 0; i <= fix_event->ncoincident; i++) {
+      int frac_rand = static_cast<int> (random_clock->uniform() * t_event);
+      frac_t_event = MIN(frac_t_event,frac_rand);
+    }
+    int decrement = frac_t_event*universe->nworlds;
+    clock -= decrement;
+
+    // share event across replicas
+    // NOTE: would be potentially more efficient for correlated events
+    //   if don't share until correlated check below has completed
     // this will complicate the dump (always on replica 0)
 
-    share_event(ireplica,1);
+    share_event(ireplica,1,decrement);
     log_event();
 
     int restart_flag = 0;
@@ -339,15 +363,16 @@ void PRD::command(int narg, char **arg)
         restart_flag = 0;
         break;
       }
-      dynamics();
+      dynamics(t_event,time_dynamics);
-      fix_event->store_state();
+      fix_event->store_state_quench();
       quench();
+      clock += t_event;
       int corr_event_check = check_event(ireplica);
       if (corr_event_check >= 0) {
-        share_event(ireplica,2);
+        share_event(ireplica,2,0);
         log_event();
         corr_endstep = update->ntimestep + t_corr;
-      } else fix_event->restore_state();
+      } else fix_event->restore_state_quench();
     }
 
     // full init/setup since are starting all replicas after event
@@ -378,8 +403,13 @@ void PRD::command(int narg, char **arg)
       timer->barrier_stop(TIME_LOOP);
       time_output += timer->array[TIME_LOOP];
     }
+
+    if (stepmode == 0) istep = update->ntimestep - update->beginstep;
+    else istep = clock;
   }
 
+  if (stepmode) nsteps = update->ntimestep - update->beginstep;
+
   // set total timers and counters so Finish() will process them
 
   timer->array[TIME_LOOP] = time_start;
@@ -407,6 +437,11 @@ void PRD::command(int narg, char **arg)
               timer->array[TIME_LOOP],nprocs_universe,nsteps,atom->natoms);
   }
 
+  if (me == 0) {
+    if (screen) fprintf(screen,"\nPRD done\n");
+    if (logfile) fprintf(logfile,"\nPRD done\n");
+  }
+
   finish->end(2);
 
   update->whichflag = 0;
@@ -430,6 +465,7 @@ void PRD::command(int narg, char **arg)
   delete [] id_compute;
   MPI_Comm_free(&comm_replica);
   delete random_select;
+  delete random_clock;
   delete random_dephase;
   delete velocity;
   delete finish;
@@ -447,25 +483,36 @@ void PRD::dephase()
 {
   bigint ntimestep_hold = update->ntimestep;
 
-  update->whichflag = 1;
+  // n_dephase iterations of dephasing, each of t_dephase steps
-  update->nsteps = n_dephase*t_dephase;
-
-  timer->barrier_start(TIME_LOOP);
 
   for (int i = 0; i < n_dephase; i++) {
-    int seed = static_cast<int> (random_dephase->uniform() * MAXSMALLINT);
+
-    if (seed == 0) seed = 1;
+    fix_event->store_state_dephase();
-    velocity->create(temp_dephase,seed);
+
-    update->integrate->run(t_dephase);
+    // do not proceed to next iteration until an event-free run occurs
-    if (temp_flag == 0) temp_dephase = temperature->compute_scalar();
+
+    int done = 0;
+    while (!done) {
+      int seed = static_cast<int> (random_dephase->uniform() * MAXSMALLINT);
+      if (seed == 0) seed = 1;
+      velocity->create(temp_dephase,seed);
+
+      dynamics(t_dephase,time_dephase);
+      fix_event->store_state_quench();
+      quench();
+
+      if (compute_event->compute_scalar() > 0.0) {
+        fix_event->restore_state_dephase();
+        update->ntimestep -= t_dephase;
+      } else {
+        fix_event->restore_state_quench();
+        done = 1;
+      }
+
+      if (temp_flag == 0) temp_dephase = temperature->compute_scalar();
+    }
   }
 
-  timer->barrier_stop(TIME_LOOP);
-  time_dephase += timer->array[TIME_LOOP];
-
-  update->integrate->cleanup();
-  finish->end(0);
-
   // reset timestep as if dephase did not occur
   // clear timestep storage from computes, since now invalid
 
@@ -475,13 +522,13 @@ void PRD::dephase()
 }
 
 /* ----------------------------------------------------------------------
-   single short dynamics run
+   short dynamics run: for event search, decorrelation, or dephasing
 ------------------------------------------------------------------------- */
 
-void PRD::dynamics()
+void PRD::dynamics(int nsteps, double &time_category)
 {
   update->whichflag = 1;
-  update->nsteps = t_event;
+  update->nsteps = nsteps;
 
   lmp->init();
   update->integrate->setup();
@@ -490,9 +537,9 @@ void PRD::dynamics()
   bigint ncalls = neighbor->ncalls;
 
   timer->barrier_start(TIME_LOOP);
-  update->integrate->run(t_event);
+  update->integrate->run(nsteps);
   timer->barrier_stop(TIME_LOOP);
-  time_dynamics += timer->array[TIME_LOOP];
+  time_category += timer->array[TIME_LOOP];
 
   nbuild += neighbor->ncalls - ncalls;
   ndanger += neighbor->ndanger;
@@ -542,7 +589,7 @@ void PRD::quench()
   update->minimize->cleanup();
   finish->end(0);
 
-  // reset timestep as if dephase did not occur
+  // reset timestep as if quench did not occur
   // clear timestep storage from computes, since now invalid
 
   update->ntimestep = ntimestep_hold;
@@ -618,7 +665,7 @@ int PRD::check_event(int replica_num)
    flag = 2 = called during PRD run = correlated event
 ------------------------------------------------------------------------- */
 
-void PRD::share_event(int ireplica, int flag)
+void PRD::share_event(int ireplica, int flag, int decrement)
 {
   timer->barrier_start(TIME_LOOP);
 
@@ -641,8 +688,11 @@ void PRD::share_event(int ireplica, int flag)
   // if this is a correlated event, time elapsed only on one partition
 
   if (flag != 2) delta *= universe->nworlds;
+  if (delta > 0 && flag != 2) delta -= decrement;
   delta += corr_adjust;
-
+  
+  // delta passed to store_event_prd() should make its clock update
+  //   be consistent with clock in main PRD loop
   // don't change the clock or timestep if this is a restart
 
   if (flag == 0 && fix_event->event_number != 0)
@@ -672,7 +722,7 @@ void PRD::share_event(int ireplica, int flag)
 
   // restore and communicate hot coords to all replicas
 
-  fix_event->restore_state();
+  fix_event->restore_state_quench();
   timer->barrier_start(TIME_LOOP);
   replicate(ireplica);
   timer->barrier_stop(TIME_LOOP);
@@ -689,7 +739,7 @@ void PRD::log_event()
   if (universe->me == 0) {
     if (universe->uscreen)
       fprintf(universe->uscreen,
-              BIGINT_FORMAT " %.3f %d %d %d %d %d\n",
+              BIGINT_FORMAT " %.3f " BIGINT_FORMAT " %d %d %d %d\n",
               fix_event->event_timestep,
               timer->elapsed(TIME_LOOP),
               fix_event->clock,
@@ -698,7 +748,7 @@ void PRD::log_event()
               fix_event->replica_number);
     if (universe->ulogfile)
       fprintf(universe->ulogfile,
-              BIGINT_FORMAT " %.3f %d %d %d %d %d\n",
+              BIGINT_FORMAT " %.3f " BIGINT_FORMAT " %d %d %d %d\n",
               fix_event->event_timestep,
               timer->elapsed(TIME_LOOP),
               fix_event->clock,
@@ -790,6 +840,7 @@ void PRD::options(int narg, char **arg)
   maxiter = 40;
   maxeval = 50;
   temp_flag = 0;
+  stepmode = 0;
 
   char *str = (char *) "geom";
   int n = strlen(str) + 1;
@@ -840,6 +891,14 @@ void PRD::options(int narg, char **arg)
       strcpy(dist_setting,arg[iarg+2]);
 
       iarg += 3;
+
+    } else if (strcmp(arg[iarg],"time") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal prd command");
+      if (strcmp(arg[iarg+1],"steps") == 0) stepmode = 0;
+      else if (strcmp(arg[iarg+1],"clock") == 0) stepmode = 1;
+      else error->all(FLERR,"Illegal prd command");
+      iarg += 2;
+
     } else error->all(FLERR,"Illegal prd command");
   }
 }

src/REPLICA/prd.h

@@ -32,9 +32,9 @@ class PRD : protected Pointers {
 
  private:
   int me,nprocs;
-  int nsteps,t_event,n_dephase,t_dephase,t_corr;
+  int t_event,n_dephase,t_dephase,t_corr;
   double etol,ftol,temp_dephase;
-  int maxiter,maxeval,temp_flag;
+  int maxiter,maxeval,temp_flag,stepmode;
   char *loop_setting,*dist_setting;
 
   int equal_size_replicas,natoms;
@@ -52,7 +52,7 @@ class PRD : protected Pointers {
 
   int ncoincident;
 
-  class RanPark *random_select;
+  class RanPark *random_select,*random_clock;
   class RanMars *random_dephase;
   class Compute *compute_event;
   class FixEventPRD *fix_event;
@@ -61,10 +61,10 @@ class PRD : protected Pointers {
   class Finish *finish;
 
   void dephase();
-  void dynamics();
+  void dynamics(int, double &);
   void quench();
   int check_event(int replica = -1);
-  void share_event(int, int);
+  void share_event(int, int, int);
   void log_event();
   void replicate(int);
   void options(int, char **);

src/REPLICA/tad.cpp

@@ -247,7 +247,7 @@ void TAD::command(int narg, char **arg)
   // This should work with if uncommented, but does not
   // if (universe->iworld == 0) {
 
-  fix_event->store_state();
+  fix_event->store_state_quench();
   quench();
 
   timer->init();
@@ -255,7 +255,7 @@ void TAD::command(int narg, char **arg)
   time_start = timer->array[TIME_LOOP];
   fix_event->store_event_tad(update->ntimestep);
   log_event(0);
-  fix_event->restore_state();
+  fix_event->restore_state_quench();
 
   // do full init/setup
 
@@ -289,7 +289,7 @@ void TAD::command(int narg, char **arg)
         while (update->ntimestep < update->endstep) {
 
           dynamics();
-          fix_event->store_state();
+          fix_event->store_state_quench();
           quench();
 
           event_flag = check_event();
@@ -299,7 +299,7 @@ void TAD::command(int narg, char **arg)
 
           // restore hot state
 
-          fix_event->restore_state();
+          fix_event->restore_state_quench();
 
           // store hot state in revert
 
@@ -404,6 +404,11 @@ void TAD::command(int narg, char **arg)
 
   if (me_universe == 0) fclose(ulogfile_neb);
 
+  if (me == 0) {
+    if (screen) fprintf(screen,"\nTAD done\n");
+    if (logfile) fprintf(logfile,"\nTAD done\n");
+  }
+
   finish->end(3);
 
   update->whichflag = 0;
@@ -929,8 +934,8 @@ void TAD::add_event()
 
   // store hot state for new event
 
-  fix_event->restore_state();
+  fix_event->restore_state_quench();
-  fix_event_list[ievent]->store_state();
+  fix_event_list[ievent]->store_state_quench();
 
   // string clean-up
 
@@ -1019,8 +1024,8 @@ void TAD::perform_event(int ievent)
 
   // load and store hot state
 
-  fix_event_list[ievent]->restore_state();
+  fix_event_list[ievent]->restore_state_quench();
-  fix_event->store_state();
+  fix_event->store_state_quench();
 }
 
 /* ----------------------------------------------------------------------