diff --git a/src/Makefile b/src/Makefile index 35c6a4efff..0fad2a5712 100755 --- a/src/Makefile +++ b/src/Makefile @@ -14,7 +14,7 @@ OBJ = $(SRC:.cpp=.o) # Package variables PACKAGE = asphere class2 colloid dipole dpd gpu granular \ - kspace manybody meam molecule opt peri poems reax xtc + kspace manybody meam molecule opt peri poems prd reax xtc PACKUSER = user-ackland user-atc user-cg-cmm user-ewaldn user-smd diff --git a/src/Makefile.package b/src/Makefile.package index e8b9bb6bd2..0c8aa9d850 100644 --- a/src/Makefile.package +++ b/src/Makefile.package @@ -1,9 +1,9 @@ # Settings for libraries used by specific LAMMPS packages # this file is auto-edited when those packages are included/excluded -PKG_INC = -I../../lib/atc -I../../lib/reax -I../../lib/poems -I../../lib/meam -PKG_PATH = -L../../lib/atc -L../../lib/reax -L../../lib/poems -L../../lib/meam -L../../lib/gpu -PKG_LIB = -latc -lreax -lpoems -lmeam -lgpu +PKG_INC = -I../../lib/reax -I../../lib/poems -I../../lib/meam +PKG_PATH = -L../../lib/reax -L../../lib/poems -L../../lib/meam +PKG_LIB = -lreax -lpoems -lmeam -PKG_SYSPATH = $(user-atc_SYSPATH) $(reax_SYSPATH) $(meam_SYSPATH) $(gpu_SYSPATH) -PKG_SYSLIB = $(user-atc_SYSLIB) $(reax_SYSLIB) $(meam_SYSLIB) $(gpu_SYSLIB) +PKG_SYSPATH = $(reax_SYSPATH) $(meam_SYSPATH) +PKG_SYSLIB = $(reax_SYSLIB) $(meam_SYSLIB) diff --git a/src/PRD/Install.csh b/src/PRD/Install.csh new file mode 100644 index 0000000000..3e73cf3c58 --- /dev/null +++ b/src/PRD/Install.csh @@ -0,0 +1,28 @@ +# Install/unInstall package classes in LAMMPS + +if ($1 == 1) then + + cp style_prd.h .. + + cp compute_event_displace.cpp .. + cp fix_event.cpp .. + cp prd.cpp .. + + cp compute_event_displace.h .. + cp fix_event.h .. + cp prd.h .. + +else if ($1 == 0) then + + rm ../style_prd.h + touch ../style_prd.h + + rm ../compute_event_displace.cpp + rm ../fix_event.cpp + rm ../prd.cpp + + rm ../compute_event_displace.h + rm ../fix_event.h + rm ../prd.h + +endif diff --git a/src/PRD/compute_event_displace.cpp b/src/PRD/compute_event_displace.cpp new file mode 100644 index 0000000000..70e14c690b --- /dev/null +++ b/src/PRD/compute_event_displace.cpp @@ -0,0 +1,150 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Mike Brown (SNL) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdlib.h" +#include "string.h" +#include "compute_event_displace.h" +#include "atom.h" +#include "domain.h" +#include "modify.h" +#include "fix.h" +#include "memory.h" +#include "error.h" +#include "update.h" + +using namespace LAMMPS_NS; + +#define INVOKED_SCALAR 1 + +/* ---------------------------------------------------------------------- */ + +ComputeEventDisplace::ComputeEventDisplace(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg != 4) error->all("Illegal compute event/displace command"); + + scalar_flag = 1; + extscalar = 0; + + double displace_dist = atof(arg[3]); + if (displace_dist <= 0.0) + error->all("Distnace must be > 0 for compute event/displace"); + displace_distsq = displace_dist * displace_dist; + + // fix event ID will be set later by PRD + + id_event = NULL; +} + +/* ---------------------------------------------------------------------- */ + +ComputeEventDisplace::~ComputeEventDisplace() +{ + delete [] id_event; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeEventDisplace::init() +{ + // set fix which stores original atom coords + // check if is correct style + + if (id_event == NULL) + error->all("Compute event/displace has not had fix event assigned"); + + int ifix = modify->find_fix(id_event); + if (ifix < 0) error->all("Could not find compute event/displace fix ID"); + fix = modify->fix[ifix]; + + if (strcmp(fix->style,"EVENT") != 0) + error->all("Compute event/displace has invalid fix event assigned"); + + triclinic = domain->triclinic; +} + +/* ---------------------------------------------------------------------- + return non-zero if an atom has moved > displace_dist since last event +------------------------------------------------------------------------- */ + +double ComputeEventDisplace::compute_scalar() +{ + invoked_scalar = update->ntimestep; + + double event = 0.0; + double **xevent = fix->vector_atom; + + double **x = atom->x; + int *mask = atom->mask; + int *image = atom->image; + int nlocal = atom->nlocal; + + double *h = domain->h; + double xprd = domain->xprd; + double yprd = domain->yprd; + double zprd = domain->zprd; + int xbox,ybox,zbox; + double dx,dy,dz,rsq; + + if (triclinic == 0) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + xbox = (image[i] & 1023) - 512; + ybox = (image[i] >> 10 & 1023) - 512; + zbox = (image[i] >> 20) - 512; + dx = x[i][0] + xbox*xprd - xevent[i][0]; + dy = x[i][1] + ybox*yprd - xevent[i][1]; + dz = x[i][2] + zbox*zprd - xevent[i][2]; + rsq = dx*dx + dy*dy + dz*dz; + if (rsq >= displace_distsq) { + event = 1.0; + break; + } + } + + } else { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + xbox = (image[i] & 1023) - 512; + ybox = (image[i] >> 10 & 1023) - 512; + zbox = (image[i] >> 20) - 512; + dx = x[i][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox - xevent[i][0]; + dy = x[i][1] + h[1]*ybox + h[3]*zbox - xevent[i][1]; + dz = x[i][2] + h[2]*zbox - xevent[i][2]; + rsq = dx*dx + dy*dy + dz*dz; + if (rsq >= displace_distsq) { + event = 1.0; + break; + } + } + } + + MPI_Allreduce(&event,&scalar,1,MPI_DOUBLE,MPI_SUM,world); + return scalar; +} + + +/* ---------------------------------------------------------------------- */ + +void ComputeEventDisplace::reset_extra_compute_fix(char *id_new) +{ + delete [] id_event; + int n = strlen(id_new) + 1; + id_event = new char[n]; + strcpy(id_event,id_new); +} diff --git a/src/PRD/compute_event_displace.h b/src/PRD/compute_event_displace.h new file mode 100644 index 0000000000..ece5a262e9 --- /dev/null +++ b/src/PRD/compute_event_displace.h @@ -0,0 +1,38 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef COMPUTE_EVENT_DISPLACE_H +#define COMPUTE_EVENT_DISPLACE_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeEventDisplace : public Compute { + public: + ComputeEventDisplace(class LAMMPS *, int, char **); + ~ComputeEventDisplace(); + void init(); + double compute_scalar(); + void reset_extra_compute_fix(char *); + + private: + int triclinic; + double displace_distsq; + char *id_event; + class Fix *fix; +}; + +} + +#endif diff --git a/src/PRD/fix_event.cpp b/src/PRD/fix_event.cpp new file mode 100644 index 0000000000..263cbabba4 --- /dev/null +++ b/src/PRD/fix_event.cpp @@ -0,0 +1,249 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Mike Brown (SNL) +------------------------------------------------------------------------- */ + +#include "stdlib.h" +#include "string.h" +#include "fix_event.h" +#include "atom.h" +#include "update.h" +#include "domain.h" +#include "neighbor.h" +#include "comm.h" +#include "universe.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +FixEvent::FixEvent(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + if (narg != 3) error->all("Illegal fix event command"); + + restart_global = 1; + + // perform initial allocation of atom-based array + // register with Atom class + + xevent = NULL; + xold = NULL; + imageold = NULL; + grow_arrays(atom->nmax); + atom->add_callback(0); + + event_number = 0; + event_timestep = 0; + clock = 0; +} + +/* ---------------------------------------------------------------------- */ + +FixEvent::~FixEvent() +{ + // unregister callbacks to this fix from Atom class + + atom->delete_callback(id,0); + + // delete locally stored array + + memory->destroy_2d_double_array(xevent); + memory->destroy_2d_double_array(xold); + memory->sfree(imageold); +} + +/* ---------------------------------------------------------------------- */ + +int FixEvent::setmask() +{ + return 0; +} + +/* ---------------------------------------------------------------------- + save current atom coords as an event + called when an event occurs in some replica + set event_timestep = when event occurred in a particular replica + update clock = elapsed time since last event, across all replicas +------------------------------------------------------------------------- */ + +void FixEvent::store_event(int timestep, int delta_clock) +{ + double **x = atom->x; + int *image = atom->image; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + domain->unmap(x[i],image[i],xevent[i]); + + event_timestep = timestep; + clock += delta_clock; + event_number++; +} + +/* ---------------------------------------------------------------------- + store state of all atoms + called before quench and subsequent check for event + so can later restore pre-quench state if no event occurs +------------------------------------------------------------------------- */ + +void FixEvent::store_state() +{ + double **x = atom->x; + double **f = atom->f; + int *image = atom->image; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + xold[i][0] = x[i][0]; + xold[i][1] = x[i][1]; + xold[i][2] = x[i][2]; + imageold[i] = image[i]; + } +} + +/* ---------------------------------------------------------------------- + restore state of all atoms to pre-quench state + called after no event detected so can continue +------------------------------------------------------------------------- */ + +void FixEvent::restore_state() +{ + double **x = atom->x; + int *image = atom->image; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + x[i][0] = xold[i][0]; + x[i][1] = xold[i][1]; + x[i][2] = xold[i][2]; + image[i] = imageold[i]; + } +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based array +------------------------------------------------------------------------- */ + +double FixEvent::memory_usage() +{ + double bytes = 6*atom->nmax * sizeof(double); + bytes += atom->nmax*sizeof(int); + return bytes; +} + +/* ---------------------------------------------------------------------- + allocate atom-based array +------------------------------------------------------------------------- */ + +void FixEvent::grow_arrays(int nmax) +{ + xevent = memory->grow_2d_double_array(xevent,nmax,3,"event:xevent"); + xold = memory->grow_2d_double_array(xold,nmax,3,"event:xold"); + imageold = (int *) + memory->srealloc(imageold,nmax*sizeof(int),"event:imageold"); + + // allow compute event to access stored event coords + + vector_atom = xevent; +} + +/* ---------------------------------------------------------------------- + copy values within local atom-based array +------------------------------------------------------------------------- */ + +void FixEvent::copy_arrays(int i, int j) +{ + xevent[j][0] = xevent[i][0]; + xevent[j][1] = xevent[i][1]; + xevent[j][2] = xevent[i][2]; + xold[j][0] = xold[i][0]; + xold[j][1] = xold[i][1]; + xold[j][2] = xold[i][2]; + imageold[j] = imageold[i]; +} + +/* ---------------------------------------------------------------------- + pack values in local atom-based array for exchange with another proc +------------------------------------------------------------------------- */ + +int FixEvent::pack_exchange(int i, double *buf) +{ + buf[0] = xevent[i][0]; + buf[1] = xevent[i][1]; + buf[2] = xevent[i][2]; + buf[3] = xold[i][0]; + buf[4] = xold[i][1]; + buf[5] = xold[i][2]; + buf[6] = imageold[i]; + + return 7; +} + +/* ---------------------------------------------------------------------- + unpack values in local atom-based array from exchange with another proc +------------------------------------------------------------------------- */ + +int FixEvent::unpack_exchange(int nlocal, double *buf) +{ + xevent[nlocal][0] = buf[0]; + xevent[nlocal][1] = buf[1]; + xevent[nlocal][2] = buf[2]; + xold[nlocal][0] = buf[3]; + xold[nlocal][1] = buf[4]; + xold[nlocal][2] = buf[5]; + imageold[nlocal] = static_cast(buf[6]); + + return 7; +} + +/* ---------------------------------------------------------------------- + pack entire state of Fix into one write +------------------------------------------------------------------------- */ + +void FixEvent::write_restart(FILE *fp) +{ + int n = 0; + double list[5]; + list[n++] = event_number; + list[n++] = event_timestep; + list[n++] = clock; + list[n++] = replica_number; + list[n++] = correlated_event; + + if (comm->me == 0) { + int size = n * sizeof(double); + fwrite(&size,sizeof(int),1,fp); + fwrite(&list,sizeof(double),n,fp); + } +} + +/* ---------------------------------------------------------------------- + use state info from restart file to restart the Fix +------------------------------------------------------------------------- */ + +void FixEvent::restart(char *buf) +{ + int n = 0; + double *list = (double *) buf; + + event_number = 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++]); +} diff --git a/src/PRD/fix_event.h b/src/PRD/fix_event.h new file mode 100644 index 0000000000..ce58d5c6f4 --- /dev/null +++ b/src/PRD/fix_event.h @@ -0,0 +1,55 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef FIX_EVENT_H +#define FIX_EVENT_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixEvent : public Fix { + public: + int event_number; // event counter + int event_timestep; // timestep of last event on any replica + int 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 + + FixEvent(class LAMMPS *, int, char **); + ~FixEvent(); + int setmask(); + + double memory_usage(); + void grow_arrays(int); + void copy_arrays(int, int); + int pack_exchange(int, double *); + int unpack_exchange(int, double *); + void write_restart(FILE *); + void restart(char *); + + // methods specific to FixEvent, invoked by PRD + + void store_event(int, int); + void store_state(); + void restore_state(); + + private: + double **xevent; // atom coords at last event + double **xold; // atom coords for reset/restore + int *imageold; // image flags for reset/restore +}; + +} + +#endif diff --git a/src/PRD/prd.cpp b/src/PRD/prd.cpp new file mode 100644 index 0000000000..83de62c474 --- /dev/null +++ b/src/PRD/prd.cpp @@ -0,0 +1,772 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Mike Brown (SNL) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdlib.h" +#include "string.h" +#include "prd.h" +#include "universe.h" +#include "update.h" +#include "atom.h" +#include "domain.h" +#include "comm.h" +#include "velocity.h" +#include "integrate.h" +#include "min.h" +#include "neighbor.h" +#include "modify.h" +#include "compute.h" +#include "fix.h" +#include "fix_event.h" +#include "force.h" +#include "pair.h" +#include "random_park.h" +#include "random_mars.h" +#include "output.h" +#include "dump.h" +#include "finish.h" +#include "timer.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define MAXINT 0x7FFFFFFF + +/* ---------------------------------------------------------------------- */ + +PRD::PRD(LAMMPS *lmp) : Pointers(lmp) {} + +/* ---------------------------------------------------------------------- + perform PRD +------------------------------------------------------------------------- */ + +void PRD::command(int narg, char **arg) +{ + int i,flag,allflag,ireplica; + + // error checks + + if (domain->box_exist == 0) + error->all("PRD command before simulation box is defined"); + if (universe->nworlds != universe->nprocs && + atom->map_style == 0) + error->all("Cannot use PRD command with multi-proc replicas " + "unless atom map exists"); + if (universe->nworlds == 1 && comm->me == 0) + error->warning("Running PRD with only one replica"); + + if (narg < 7) error->universe_all("Illegal prd command"); + + nsteps = atoi(arg[0]); + t_event = atoi(arg[1]); + n_dephase = atoi(arg[2]); + t_dephase = atoi(arg[3]); + t_corr = atoi(arg[4]); + + char id_compute[strlen(arg[5])+1]; + strcpy(id_compute,arg[5]); + int seed = atoi(arg[6]); + + options(narg-7,&arg[7]); + + // total # of timesteps must be multiple of t_event + + if (t_event == 0) error->universe_all("Invalid t_event in prd command"); + if (nsteps % t_event) + error->universe_all("PRD nsteps must be multiple of t_event"); + if (t_corr % t_event) + error->universe_all("PRD t_corr must be multiple of t_event"); + + // local storage + + int me_universe = universe->me; + int nprocs_universe = universe->nprocs; + int nreplica = universe->nworlds; + int iworld = universe->iworld; + + MPI_Comm_rank(world,&me); + MPI_Comm_size(world,&nprocs); + + // comm_replica = communicator between same proc across replicas + // not used if replicas have unequal number of procs + // equal_size_replicas = 1 if all replicas have same # of procs + + int color = me; + MPI_Comm_split(universe->uworld,color,0,&comm_replica); + + flag = 0; + if (nreplica*nprocs == nprocs_universe) flag = 1; + MPI_Allreduce(&flag,&equal_size_replicas,1,MPI_INT,MPI_MIN,universe->uworld); + + // workspace for inter-replica communication via gathers + + natoms = static_cast (atom->natoms); + + displacements = NULL; + tagall = NULL; + xall = NULL; + imageall = NULL; + + if (nreplica != nprocs_universe) { + displacements = new int[nprocs]; + tagall = (int *) memory->smalloc(natoms*sizeof(int),"prd:tagall"); + xall = memory->create_2d_double_array(natoms,3,"prd:xall"); + imageall = (int *) memory->smalloc(natoms*sizeof(int),"prd:imageall"); + } + + // random_select = same RNG for each replica for multiple event selection + // random_dephase = unique RNG for each replica for dephasing + + random_select = new RanPark(lmp,seed); + random_dephase = new RanMars(lmp,seed+iworld); + + // create ComputeTemp class to monitor temperature + + char **args = new char*[3]; + args[0] = (char *) "prd_temp"; + args[1] = (char *) "all"; + args[2] = (char *) "temp"; + modify->add_compute(3,args); + temperature = modify->compute[modify->ncompute-1]; + + // create Velocity class for velocity creation in dephasing + // pass it temperature compute, loop_setting, dist_setting settings + + atom->check_mass(); + velocity = new Velocity(lmp); + velocity->init_external("all"); + + args[0] = (char *) "temp"; + args[1] = (char *) "prd_temp"; + velocity->options(2,args); + args[0] = (char *) "loop"; + args[1] = (char *) loop_setting; + if (loop_setting) velocity->options(2,args); + args[0] = (char *) "dist"; + args[1] = (char *) dist_setting; + if (dist_setting) velocity->options(2,args); + + // create FixEvent class to store event and pre-quench states + + args[0] = (char *) "prd_event"; + args[1] = (char *) "all"; + args[2] = (char *) "EVENT"; + modify->add_fix(3,args); + fix_event = (FixEvent *) modify->fix[modify->nfix-1]; + + // create Finish for timing output + + finish = new Finish(lmp); + + // string clean-up + + delete [] args; + delete [] loop_setting; + delete [] dist_setting; + + // assign FixEvent to event-detection compute + // necessary so it will know atom coords at last event + + int icompute = modify->find_compute(id_compute); + if (icompute < 0) error->all("Could not find compute ID for PRD"); + compute_event = modify->compute[icompute]; + compute_event->reset_extra_compute_fix("prd_event"); + + // reset reneighboring criteria since will perform minimizations + + neigh_every = neighbor->every; + neigh_delay = neighbor->delay; + neigh_dist_check = neighbor->dist_check; + + if (neigh_every != 1 || neigh_delay != 0 || neigh_dist_check != 1) { + if (me == 0) + error->warning("Resetting reneighboring criteria during PRD"); + } + + neighbor->every = 1; + neighbor->delay = 0; + neighbor->dist_check = 1; + + // initialize PRD as if one long dynamics run + + update->whichflag = 1; + update->nsteps = nsteps; + update->beginstep = update->firststep = update->ntimestep; + update->endstep = update->laststep = update->firststep + nsteps; + update->restrict_output = 1; + + lmp->init(); + + // init minimizer settings and minimizer itself + + update->etol = etol; + update->ftol = ftol; + update->max_eval = maxeval; + + update->minimize->init(); + + // cannot use PRD with time-dependent fixes + + for (int i = 0; i < modify->nfix; i++) + if (modify->fix[i]->time_depend) + error->all("Cannot use PRD with a time-dependent fix defined"); + + // perform PRD simulation + + if (me_universe == 0 && universe->uscreen) + fprintf(universe->uscreen,"Setting up PRD ...\n"); + + if (me_universe == 0) { + if (universe->uscreen) + fprintf(universe->uscreen,"Step Clock Event Correlated Replica\n"); + if (universe->ulogfile) + fprintf(universe->ulogfile,"Step Clock Event Correlated Replica\n"); + } + + // store hot state and quenched event for replica 0 + // use share_event() to copy that info to all replicas + // this insures all start from same place + + // need this line if quench() does only setup_minimal() + //update->minimize->setup(); + + fix_event->store_state(); + quench(); + share_event(0,0); + log_event(); + + // do full init/setup since are starting all replicas after event + // replica 0 bcasts temp to all replicas if user did not set temp_dephase + update->whichflag = 1; + lmp->init(); + update->integrate->setup(); + if (temp_flag == 0) { + if (universe->iworld == 0) + temp_dephase = temperature->compute_scalar(); + MPI_Bcast(&temp_dephase,1,MPI_DOUBLE,universe->root_proc[0], + universe->uworld); + } + + // main loop: look for events until out of time + // (1) dephase independently on each proc after event + // (2) loop: dynamics, store state, quench, check event, restore state + // (3) share and record event + + nbuild = ndanger = 0; + time_dephase = time_dynamics = time_quench = time_comm = time_output = 0.0; + + timer->barrier_start(TIME_LOOP); + double time_start = timer->array[TIME_LOOP]; + + while (update->ntimestep < update->endstep) { + dephase(); + + ireplica = -1; + while (update->ntimestep < update->endstep) { + dynamics(); + fix_event->store_state(); + quench(); + ireplica = check_event(); + if (ireplica >= 0) break; + fix_event->restore_state(); + } + if (ireplica < 0) break; + + // Can potentially be more efficient for correlated events by not + // sharing until correlated check has completed (this will complicate + // the dump (always on replica 0)). + share_event(ireplica,1); + log_event(); + + int restart_flag = 0; + if (output->restart_every && universe->iworld == 0) + if (fix_event->event_number % output->restart_every == 0) + restart_flag = 1; + + // Correlated event loop + // -- We could have other procs doing dephasing during this time + int corr_endstep = update->ntimestep + t_corr; + while (update->ntimestep < corr_endstep) { + if (update->ntimestep == update->endstep) { + restart_flag = 0; + break; + } + dynamics(); + fix_event->store_state(); + quench(); + int corr_event_check = check_event(ireplica); + if (corr_event_check >= 0) { + share_event(ireplica,2); + log_event(); + corr_endstep = update->ntimestep + t_corr; + } else + fix_event->restore_state(); + } + + // do full init/setup since are starting all replicas after event + // -- also, need to get reneighbor before restart + update->whichflag = 1; + lmp->init(); + update->integrate->setup(); + + timer->barrier_start(TIME_LOOP); + if (t_corr > 0) replicate(ireplica); + + // event replica bcasts temp to all replicas if user did not set temp_dephase + if (temp_flag == 0) { + if (ireplica == universe->iworld) + temp_dephase = temperature->compute_scalar(); + MPI_Bcast(&temp_dephase,1,MPI_DOUBLE,universe->root_proc[ireplica], + universe->uworld); + } + timer->barrier_stop(TIME_LOOP); + time_comm += timer->array[TIME_LOOP]; + + // write restart file of hot coords + if (restart_flag) { + timer->barrier_start(TIME_LOOP); + output->write_restart(update->ntimestep); + timer->barrier_stop(TIME_LOOP); + time_output += timer->array[TIME_LOOP]; + } + + } + + // set total timers and counters so Finish() will process them + + timer->array[TIME_LOOP] = time_start; + timer->barrier_stop(TIME_LOOP); + + timer->array[TIME_PAIR] = time_dephase; + timer->array[TIME_BOND] = time_dynamics; + timer->array[TIME_KSPACE] = time_quench; + timer->array[TIME_COMM] = time_comm; + timer->array[TIME_OUTPUT] = time_output; + + neighbor->ncalls = nbuild; + neighbor->ndanger = ndanger; + + if (me_universe == 0) { + if (universe->uscreen) + fprintf(universe->uscreen, + "Loop time of %g on %d procs for %d steps with %.15g atoms\n", + timer->array[TIME_LOOP],nprocs_universe,nsteps,atom->natoms); + if (universe->ulogfile) + fprintf(universe->ulogfile, + "Loop time of %g on %d procs for %d steps with %.15g atoms\n", + timer->array[TIME_LOOP],nprocs_universe,nsteps,atom->natoms); + } + + finish->end(2); + + update->whichflag = 0; + update->firststep = update->laststep = 0; + update->beginstep = update->endstep = 0; + update->restrict_output = 0; + + // reset reneighboring criteria + + neighbor->every = neigh_every; + neighbor->delay = neigh_delay; + neighbor->dist_check = neigh_dist_check; + + // clean up + + delete [] displacements; + memory->sfree(tagall); + memory->destroy_2d_double_array(xall); + memory->sfree(imageall); + + MPI_Comm_free(&comm_replica); + delete random_select; + delete random_dephase; + delete velocity; + delete finish; + modify->delete_compute("prd_temp"); + modify->delete_fix("prd_event"); +} + +/* ---------------------------------------------------------------------- + dephasing = one or more short runs with new random velocities +------------------------------------------------------------------------- */ + +void PRD::dephase() +{ + int ntimestep_hold = update->ntimestep; + + update->whichflag = 1; + update->nsteps = n_dephase*t_dephase; + + timer->barrier_start(TIME_LOOP); + + for (int i = 0; i < n_dephase; i++) { + int seed = static_cast (random_dephase->uniform() * MAXINT); + if (seed == 0) seed = 1; + velocity->create(temp_dephase,seed); + update->integrate->run(t_dephase); + 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 + + update->ntimestep = ntimestep_hold; + for (int i = 0; i < modify->ncompute; i++) + if (modify->compute[i]->timeflag) modify->compute[i]->clearstep(); +} + +/* ---------------------------------------------------------------------- + single short dynamics run +------------------------------------------------------------------------- */ + +void PRD::dynamics() +{ + update->whichflag = 1; + update->nsteps = t_event; + + lmp->init(); + update->integrate->setup(); + //modify->addstep_compute_all(update->ntimestep); + int ncalls = neighbor->ncalls; + + timer->barrier_start(TIME_LOOP); + update->integrate->run(t_event); + timer->barrier_stop(TIME_LOOP); + time_dynamics += timer->array[TIME_LOOP]; + + nbuild += neighbor->ncalls - ncalls; + ndanger += neighbor->ndanger; + + update->integrate->cleanup(); + finish->end(0); +} + +/* ---------------------------------------------------------------------- + quench minimization +------------------------------------------------------------------------- */ + +void PRD::quench() +{ + int ntimestep_hold = update->ntimestep; + + // need to change whichflag so that minimize->setup() calling + // modify->setup() will call fix->min_setup() + + update->whichflag = 2; + update->nsteps = maxiter; + + // these work + lmp->init(); + update->minimize->setup(); + + // these do not work + //modify->addstep_compute_all(update->ntimestep); + //update->minimize->setup_minimal(1); + + int ncalls = neighbor->ncalls; + + timer->barrier_start(TIME_LOOP); + update->minimize->run(maxiter); + timer->barrier_stop(TIME_LOOP); + time_quench += timer->array[TIME_LOOP]; + + if (neighbor->ncalls == ncalls) quench_reneighbor = 0; + else quench_reneighbor = 1; + + update->minimize->cleanup(); + finish->end(0); + + // reset timestep as if dephase did not occur + // clear timestep storage from computes, since now invalid + + update->ntimestep = ntimestep_hold; + for (int i = 0; i < modify->ncompute; i++) + if (modify->compute[i]->timeflag) modify->compute[i]->clearstep(); +} + +/* ---------------------------------------------------------------------- + check for an event in any replica + if replica_num is non-negative only check for event on replica_num + if multiple events, choose one at random + return -1 if no event + else return ireplica = world in which event occured +------------------------------------------------------------------------- */ + +int PRD::check_event(int replica_num) +{ + int worldflag,universeflag,scanflag,replicaflag,ireplica; + + worldflag = 0; + if (compute_event->compute_scalar() > 0.0) worldflag = 1; + if (replica_num >= 0 && replica_num != universe->iworld) worldflag = 0; + + timer->barrier_start(TIME_LOOP); + if (me == 0) MPI_Allreduce(&worldflag,&universeflag,1, + MPI_INT,MPI_SUM,comm_replica); + MPI_Bcast(&universeflag,1,MPI_INT,0,world); + if (!universeflag) + ireplica = -1; + else { + if (universeflag > 1) { + int iwhich = static_cast (universeflag*random_select->uniform()) + 1; + if (me == 0) MPI_Scan(&worldflag,&scanflag,1, + MPI_INT,MPI_SUM,comm_replica); + MPI_Bcast(&scanflag,1,MPI_INT,0,world); + if (scanflag != iwhich) worldflag = 0; + } + + if (worldflag) replicaflag = universe->iworld; + else replicaflag = 0; + if (me == 0) MPI_Allreduce(&replicaflag,&ireplica,1, + MPI_INT,MPI_SUM,comm_replica); + MPI_Bcast(&ireplica,1,MPI_INT,0,world); + } + + timer->barrier_stop(TIME_LOOP); + time_comm += timer->array[TIME_LOOP]; + return ireplica; +} + +/* ---------------------------------------------------------------------- + share quenched and hot coords owned by ireplica with all replicas + all replicas store event in fix_event + replica 0 dumps event snapshot + flag = 0 = called before PRD run + flag = 1 = called during PRD run = not correlated event + flag = 2 = called during PRD run = correlated event +------------------------------------------------------------------------- */ + +void PRD::share_event(int ireplica, int flag) +{ + timer->barrier_start(TIME_LOOP); + + // communicate quenched coords to all replicas and store as event + // decrement event counter if flag = 0 since not really an event + + replicate(ireplica); + timer->barrier_stop(TIME_LOOP); + time_comm += timer->array[TIME_LOOP]; + + // Adjust time for last correlated event check (Not on first event) + int corr_adjust = t_corr; + if (fix_event->event_number<1 || flag == 2) corr_adjust = 0; + // Time since last correlated event check + int delta = update->ntimestep - fix_event->event_timestep - corr_adjust; + // If this is a correlated event, time was only on one partition + if (flag != 2) delta *= universe->nworlds; + delta += corr_adjust; + // Don't change the clock or timestep if this is a restart + if (flag == 0 && fix_event->event_number != 0) + fix_event->store_event(fix_event->event_timestep,0); + else { + fix_event->store_event(update->ntimestep,delta); + fix_event->replica_number = ireplica; + fix_event->correlated_event = 0; + if (flag == 2) fix_event->correlated_event = 1; + } + if (flag == 0) fix_event->event_number--; + + // dump snapshot of quenched coords + // must reneighbor and compute forces before dumping + // since replica 0 possibly has new state from another replica + // addstep_compute_all insures eng/virial are calculated if needed + + if (output->ndump && universe->iworld == 0) { + timer->barrier_start(TIME_LOOP); + modify->addstep_compute_all(update->ntimestep); + update->integrate->setup_minimal(1); + output->write_dump(update->ntimestep); + timer->barrier_stop(TIME_LOOP); + time_output += timer->array[TIME_LOOP]; + } + + // restore and communicate hot coords to all replicas + + fix_event->restore_state(); + timer->barrier_start(TIME_LOOP); + replicate(ireplica); + timer->barrier_stop(TIME_LOOP); + time_comm += timer->array[TIME_LOOP]; +} + +/* ---------------------------------------------------------------------- + universe proc 0 prints event info +------------------------------------------------------------------------- */ + +void PRD::log_event() +{ + if (universe->me == 0) { + if (universe->uscreen) + fprintf(universe->uscreen,"%d %d %d %d %d\n", + fix_event->event_timestep,fix_event->clock, + fix_event->event_number,fix_event->correlated_event, + fix_event->replica_number); + if (universe->ulogfile) + fprintf(universe->ulogfile,"%d %d %d %d %d\n", + fix_event->event_timestep,fix_event->clock, + fix_event->event_number,fix_event->correlated_event, + fix_event->replica_number); + } +} + +/* ---------------------------------------------------------------------- + communicate atom coords and image flags in ireplica to all other replicas + one proc per replica: + direct overwrite via bcast + equal procs per replica and no replica reneighbored: + direct overwrite via bcast + unequal procs per replica or reneighboring occurred: + collect to root proc of event replica + bcast to roots of other replicas + bcast within each replica + each proc extracts info for atoms it owns using atom IDs +------------------------------------------------------------------------- */ + +void PRD::replicate(int ireplica) +{ + int nreplica = universe->nworlds; + int nprocs_universe = universe->nprocs; + int i,m,flag,commflag; + int counts[nprocs]; + + if (nreplica == nprocs_universe) commflag = 0; + else if (equal_size_replicas) { + flag = 0; + if (quench_reneighbor) flag = 1; + MPI_Allreduce(&flag,&commflag,1,MPI_INT,MPI_MAX,universe->uworld); + } else commflag = 1; + + if (commflag == 0) { + MPI_Bcast(atom->image,atom->nlocal,MPI_INT,ireplica,comm_replica); + MPI_Bcast(atom->x[0],3*atom->nlocal,MPI_DOUBLE,ireplica,comm_replica); + } else { + if (universe->iworld == ireplica) { + MPI_Gather(&atom->nlocal,1,MPI_INT,counts,1,MPI_INT,0,world); + displacements[0] = 0; + for (i = 0; i < nprocs-1; i++) + displacements[i+1] = displacements[i] + counts[i]; + MPI_Gatherv(atom->tag,atom->nlocal,MPI_INT, + tagall,counts,displacements,MPI_INT,0,world); + MPI_Gatherv(atom->image,atom->nlocal,MPI_INT, + imageall,counts,displacements,MPI_INT,0,world); + for (i = 0; i < nprocs; i++) counts[i] *= 3; + for (i = 0; i < nprocs-1; i++) + displacements[i+1] = displacements[i] + counts[i]; + MPI_Gatherv(atom->x[0],3*atom->nlocal,MPI_DOUBLE, + xall[0],counts,displacements,MPI_DOUBLE,0,world); + } + + if (me == 0) { + MPI_Bcast(tagall,natoms,MPI_INT,ireplica,comm_replica); + MPI_Bcast(imageall,natoms,MPI_INT,ireplica,comm_replica); + MPI_Bcast(xall[0],3*natoms,MPI_DOUBLE,ireplica,comm_replica); + } + + MPI_Bcast(tagall,natoms,MPI_INT,0,world); + MPI_Bcast(imageall,natoms,MPI_INT,0,world); + MPI_Bcast(xall[0],3*natoms,MPI_DOUBLE,0,world); + + double **x = atom->x; + int nlocal = atom->nlocal; + + for (i = 0; i < natoms; i++) { + m = atom->map(tagall[i]); + if (m >= 0 && m < nlocal) { + x[m][0] = xall[i][0]; + x[m][1] = xall[i][1]; + x[m][2] = xall[i][2]; + atom->image[m] = imageall[i]; + } + } + } +} + +/* ---------------------------------------------------------------------- + parse optional parameters at end of PRD input line +------------------------------------------------------------------------- */ + +void PRD::options(int narg, char **arg) +{ + if (narg < 0) error->all("Illegal prd command"); + + // set defaults + + etol = 0.1; + ftol = 0.1; + maxiter = 40; + maxeval = 50; + temp_flag = 0; + loop_setting = NULL; + dist_setting = NULL; + + int iarg = 0; + while (iarg < narg) { + if (strcmp(arg[iarg],"min") == 0) { + if (iarg+5 > narg) error->all("Illegal prd command"); + etol = atof(arg[iarg+1]); + ftol = atof(arg[iarg+2]); + maxiter = atoi(arg[iarg+3]); + maxeval = atoi(arg[iarg+4]); + if (maxiter < 0) error->all("Illegal prd command"); + iarg += 5; + + } else if (strcmp(arg[iarg],"temp") == 0) { + if (iarg+2 > narg) error->all("Illegal prd command"); + temp_flag = 1; + temp_dephase = atof(arg[iarg+1]); + if (temp_dephase <= 0.0) error->all("Illegal prd command"); + iarg += 2; + + } else if (strcmp(arg[iarg],"vel") == 0) { + if (iarg+3 > narg) error->all("Illegal prd command"); + + if (strcmp(arg[iarg+1],"all") == 0) loop_setting = NULL; + else if (strcmp(arg[iarg+1],"local") == 0) loop_setting = NULL; + else if (strcmp(arg[iarg+1],"geom") == 0) loop_setting = NULL; + else error->all("Illegal prd command"); + int n = strlen(arg[iarg+1]) + 1; + loop_setting = new char[n]; + strcpy(loop_setting,arg[iarg+1]); + + if (strcmp(arg[iarg+2],"uniform") == 0) dist_setting = NULL; + else if (strcmp(arg[iarg+2],"gaussian") == 0) dist_setting = NULL; + else error->all("Illegal prd command"); + n = strlen(arg[iarg+2]) + 1; + dist_setting = new char[n]; + strcpy(dist_setting,arg[iarg+2]); + + iarg += 3; + } else error->all("Illegal prd command"); + } + + // Set defaults + if (loop_setting == NULL) { + loop_setting = new char[5]; + strcpy(loop_setting,"geom"); + } + if (dist_setting == NULL) { + dist_setting = new char[9]; + strcpy(dist_setting,"gaussian"); + } +} diff --git a/src/PRD/prd.h b/src/PRD/prd.h new file mode 100644 index 0000000000..a2d88267ea --- /dev/null +++ b/src/PRD/prd.h @@ -0,0 +1,65 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef PRD_H +#define PRD_H + +#include "pointers.h" + +namespace LAMMPS_NS { + +class PRD : protected Pointers { + public: + PRD(class LAMMPS *); + ~PRD() {} + void command(int, char **); + + private: + int me,nprocs; + int nsteps,t_event,n_dephase,t_dephase,t_corr; + double etol,ftol,temp_dephase; + int maxiter,maxeval,temp_flag; + char *loop_setting,*dist_setting; + + int equal_size_replicas,natoms; + int neigh_every,neigh_delay,neigh_dist_check; + int nbuild,ndanger; + int quench_reneighbor; + + double time_dephase,time_dynamics,time_quench,time_comm,time_output; + + MPI_Comm comm_replica; + int *tagall,*displacements,*imageall; + double **xall; + + class RanPark *random_select; + class RanMars *random_dephase; + class Compute *compute_event; + class FixEvent *fix_event; + class Velocity *velocity; + class Compute *temperature; + class Finish *finish; + + void dephase(); + void dynamics(); + void quench(); + int check_event(int replica = -1); + void share_event(int, int); + void log_event(); + void replicate(int); + void options(int, char **); +}; + +} + +#endif diff --git a/src/PRD/style_prd.h b/src/PRD/style_prd.h new file mode 100644 index 0000000000..7db4a277ee --- /dev/null +++ b/src/PRD/style_prd.h @@ -0,0 +1,36 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef CommandInclude +#include "prd.h" +#endif + +#ifdef CommandClass +CommandStyle(prd,PRD) +#endif + +#ifdef ComputeInclude +#include "compute_event_displace.h" +#endif + +#ifdef ComputeClass +ComputeStyle(event/displace,ComputeEventDisplace) +#endif + +#ifdef FixInclude +#include "fix_event.h" +#endif + +#ifdef FixClass +FixStyle(EVENT,FixEvent) +#endif diff --git a/src/style.h b/src/style.h index 936e8bfbbc..d35ddd3cad 100644 --- a/src/style.h +++ b/src/style.h @@ -385,6 +385,7 @@ RegionStyle(union,RegUnion) #include "style_opt.h" #include "style_peri.h" #include "style_poems.h" +#include "style_prd.h" #include "style_reax.h" #include "style_xtc.h" diff --git a/src/style_user_ackland.h b/src/style_user_ackland.h index 6e7483a9f7..e69de29bb2 100644 --- a/src/style_user_ackland.h +++ b/src/style_user_ackland.h @@ -1,20 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifdef ComputeInclude -#include "compute_ackland_atom.h" -#endif - -#ifdef ComputeClass -ComputeStyle(ackland/atom,ComputeAcklandAtom) -#endif diff --git a/src/style_user_ewaldn.h b/src/style_user_ewaldn.h index 3eafa50744..e69de29bb2 100644 --- a/src/style_user_ewaldn.h +++ b/src/style_user_ewaldn.h @@ -1,30 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifdef KSpaceInclude -#include "ewald_n.h" -#endif - -#ifdef KSpaceClass -KSpaceStyle(ewald/n,EwaldN) -#endif - -#ifdef PairInclude -#include "pair_buck_coul.h" -#include "pair_lj_coul.h" -#endif - -#ifdef PairClass -PairStyle(buck/coul,PairBuckCoul) -PairStyle(lj/coul,PairLJCoul) -#endif