From 4833f14be89102644c760c8ad9a5f1e08fd2ea62 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 14 Feb 2014 16:31:23 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11565 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/REPLICA/fix_event.cpp | 103 ++++++++++++++++++------ src/REPLICA/fix_event.h | 21 +++-- src/REPLICA/fix_event_prd.cpp | 4 +- src/REPLICA/fix_event_prd.h | 2 +- src/REPLICA/prd.cpp | 145 ++++++++++++++++++++++++---------- src/REPLICA/prd.h | 10 +-- src/REPLICA/tad.cpp | 21 +++-- 7 files changed, 217 insertions(+), 89 deletions(-) diff --git a/src/REPLICA/fix_event.cpp b/src/REPLICA/fix_event.cpp index 458d5c5818..a1eb7f7894 100644 --- a/src/REPLICA/fix_event.cpp +++ b/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(buf[9]); + imageold[nlocal] = static_cast(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(buf[16]); - return 10; -} - -/* ---------------------------------------------------------------------- - 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) -{ + return 17; } diff --git a/src/REPLICA/fix_event.h b/src/REPLICA/fix_event.h index f8fe5d2d81..b607098063 100644 --- a/src/REPLICA/fix_event.h +++ b/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 restart(char *); + virtual void write_restart(FILE *) {} + virtual void restart(char *) {} // methods specific to FixEvent - void store_event(); // store quenched atoms - void restore_event(); // restore quenched atoms - void store_state(); // store hot atoms - void restore_state(); // restore hot atoms + void store_event(); // store quenched atoms + void restore_event(); // restore quenched atoms + void store_state_quench(); // store hot atoms prior to quench + 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 }; } diff --git a/src/REPLICA/fix_event_prd.cpp b/src/REPLICA/fix_event_prd.cpp index 338c6a2f76..57db1e17b7 100644 --- a/src/REPLICA/fix_event_prd.cpp +++ b/src/REPLICA/fix_event_prd.cpp @@ -91,8 +91,8 @@ void FixEventPRD::restart(char *buf) double *list = (double *) buf; event_number = static_cast (list[n++]); - event_timestep = static_cast (list[n++]); - clock = static_cast (list[n++]); + event_timestep = static_cast (list[n++]); + clock = static_cast (list[n++]); replica_number = static_cast (list[n++]); correlated_event = static_cast (list[n++]); ncoincident = static_cast (list[n++]); diff --git a/src/REPLICA/fix_event_prd.h b/src/REPLICA/fix_event_prd.h index 8cf7ef7846..d4442925d8 100644 --- a/src/REPLICA/fix_event_prd.h +++ b/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 diff --git a/src/REPLICA/prd.cpp b/src/REPLICA/prd.cpp index 145435e1a1..6e008840a8 100644 --- a/src/REPLICA/prd.cpp +++ b/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_dephase = unique RNG for each replica for dephasing + // random_select = same RNG for each replica, for multiple event selection + // 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) { - dynamics(); - fix_event->store_state(); + while (istep < nsteps) { + dynamics(t_event,time_dynamics); + 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 - // share until correlated check has completed + // decrement clock by random time at which 1 or more events occurred + + int frac_t_event = t_event; + for (int i = 0; i <= fix_event->ncoincident; i++) { + int frac_rand = static_cast (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(); - fix_event->store_state(); + dynamics(t_event,time_dynamics); + 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; - update->nsteps = n_dephase*t_dephase; - - timer->barrier_start(TIME_LOOP); + // n_dephase iterations of dephasing, each of t_dephase steps for (int i = 0; i < n_dephase; i++) { - int seed = static_cast (random_dephase->uniform() * MAXSMALLINT); - if (seed == 0) seed = 1; - velocity->create(temp_dephase,seed); - update->integrate->run(t_dephase); - if (temp_flag == 0) temp_dephase = temperature->compute_scalar(); + + fix_event->store_state_dephase(); + + // do not proceed to next iteration until an event-free run occurs + + int done = 0; + while (!done) { + int seed = static_cast (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"); } } diff --git a/src/REPLICA/prd.h b/src/REPLICA/prd.h index 74333da021..d7a58166e9 100644 --- a/src/REPLICA/prd.h +++ b/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 **); diff --git a/src/REPLICA/tad.cpp b/src/REPLICA/tad.cpp index 87e3ba0051..8b5b7ab861 100644 --- a/src/REPLICA/tad.cpp +++ b/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_list[ievent]->store_state(); + fix_event->restore_state_quench(); + 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->store_state(); + fix_event_list[ievent]->restore_state_quench(); + fix_event->store_state_quench(); } /* ----------------------------------------------------------------------