From c69f3aea57560b41e4c90581376a41c693014244 Mon Sep 17 00:00:00 2001 From: athomps Date: Wed, 5 Jan 2011 00:23:34 +0000 Subject: [PATCH] Added tad example git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@5475 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/REPLICA/Install.sh | 12 + src/REPLICA/compute_event_displace.cpp | 11 +- src/REPLICA/compute_event_displace.h | 2 +- src/REPLICA/fix_event.cpp | 108 ++- src/REPLICA/fix_event.h | 25 +- src/REPLICA/fix_event_prd.cpp | 107 +++ src/REPLICA/fix_event_prd.h | 53 ++ src/REPLICA/fix_event_tad.cpp | 100 +++ src/REPLICA/fix_event_tad.h | 51 ++ src/REPLICA/fix_neb.cpp | 31 +- src/REPLICA/neb.cpp | 124 +++- src/REPLICA/neb.h | 14 +- src/REPLICA/prd.cpp | 10 +- src/REPLICA/prd.h | 2 +- src/REPLICA/tad.cpp | 987 +++++++++++++++++++++++++ src/REPLICA/tad.h | 88 +++ 16 files changed, 1597 insertions(+), 128 deletions(-) create mode 100644 src/REPLICA/fix_event_prd.cpp create mode 100644 src/REPLICA/fix_event_prd.h create mode 100644 src/REPLICA/fix_event_tad.cpp create mode 100644 src/REPLICA/fix_event_tad.h create mode 100644 src/REPLICA/tad.cpp create mode 100644 src/REPLICA/tad.h diff --git a/src/REPLICA/Install.sh b/src/REPLICA/Install.sh index bfacd7835c..407f6d5876 100644 --- a/src/REPLICA/Install.sh +++ b/src/REPLICA/Install.sh @@ -4,32 +4,44 @@ if (test $1 = 1) then cp compute_event_displace.cpp .. cp fix_event.cpp .. + cp fix_event_prd.cpp .. + cp fix_event_tad.cpp .. cp fix_neb.cpp .. cp neb.cpp .. cp prd.cpp .. + cp tad.cpp .. cp temper.cpp .. cp compute_event_displace.h .. cp fix_event.h .. + cp fix_event_prd.h .. + cp fix_event_tad.h .. cp fix_neb.h .. cp neb.h .. cp prd.h .. + cp tad.h .. cp temper.h .. elif (test $1 = 0) then rm ../compute_event_displace.cpp rm ../fix_event.cpp + rm ../fix_event_prd.cpp + rm ../fix_event_tad.cpp rm ../fix_neb.cpp rm ../neb.cpp rm ../prd.cpp + rm ../tad.cpp rm ../temper.cpp rm ../compute_event_displace.h rm ../fix_event.h + rm ../fix_event_prd.h + rm ../fix_event_tad.h rm ../fix_neb.h rm ../neb.h rm ../prd.h + rm ../tad.h rm ../temper.h fi diff --git a/src/REPLICA/compute_event_displace.cpp b/src/REPLICA/compute_event_displace.cpp index 9a1f3ba009..53e9d9f8f9 100644 --- a/src/REPLICA/compute_event_displace.cpp +++ b/src/REPLICA/compute_event_displace.cpp @@ -23,7 +23,7 @@ #include "atom.h" #include "domain.h" #include "modify.h" -#include "fix.h" +#include "fix_event.h" #include "memory.h" #include "error.h" #include "update.h" @@ -70,9 +70,10 @@ void ComputeEventDisplace::init() if (id_event != NULL) { int ifix = modify->find_fix(id_event); if (ifix < 0) error->all("Could not find compute event/displace fix ID"); - fix = modify->fix[ifix]; + fix_event = (FixEvent*) modify->fix[ifix]; - if (strcmp(fix->style,"EVENT") != 0) + if (strcmp(fix_event->style,"EVENT/PRD") != 0 && + strcmp(fix_event->style,"EVENT/TAD") != 0) error->all("Compute event/displace has invalid fix event assigned"); } @@ -90,7 +91,7 @@ double ComputeEventDisplace::compute_scalar() if (id_event == NULL) return 0.0; double event = 0.0; - double **xevent = fix->array_atom; + double **xevent = fix_event->array_atom; double **x = atom->x; int *mask = atom->mask; @@ -119,7 +120,6 @@ double ComputeEventDisplace::compute_scalar() break; } } - } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { @@ -138,6 +138,7 @@ double ComputeEventDisplace::compute_scalar() } MPI_Allreduce(&event,&scalar,1,MPI_DOUBLE,MPI_SUM,world); + return scalar; } diff --git a/src/REPLICA/compute_event_displace.h b/src/REPLICA/compute_event_displace.h index 393e4d1bec..21cf8d09f7 100644 --- a/src/REPLICA/compute_event_displace.h +++ b/src/REPLICA/compute_event_displace.h @@ -36,7 +36,7 @@ class ComputeEventDisplace : public Compute { int triclinic; double displace_distsq; char *id_event; - class Fix *fix; + class FixEvent *fix_event; }; } diff --git a/src/REPLICA/fix_event.cpp b/src/REPLICA/fix_event.cpp index 8b3105f310..dfea7a4cec 100644 --- a/src/REPLICA/fix_event.cpp +++ b/src/REPLICA/fix_event.cpp @@ -12,7 +12,7 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing author: Mike Brown (SNL) + Contributing author: Mike Brown (SNL), Aidan Thompson (SNL) ------------------------------------------------------------------------- */ #include "stdlib.h" @@ -43,13 +43,10 @@ FixEvent::FixEvent(LAMMPS *lmp, int narg, char **arg) : xevent = NULL; xold = NULL; + vold = NULL; imageold = NULL; grow_arrays(atom->nmax); atom->add_callback(0); - - event_number = 0; - event_timestep = update->ntimestep; - clock = 0; } /* ---------------------------------------------------------------------- */ @@ -64,6 +61,7 @@ FixEvent::~FixEvent() memory->destroy_2d_double_array(xevent); memory->destroy_2d_double_array(xold); + memory->destroy_2d_double_array(vold); memory->sfree(imageold); } @@ -76,12 +74,10 @@ int FixEvent::setmask() /* ---------------------------------------------------------------------- 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 + called when an event occurs ------------------------------------------------------------------------- */ -void FixEvent::store_event(int timestep, int delta_clock) +void FixEvent::store_event() { double **x = atom->x; int *image = atom->image; @@ -90,9 +86,41 @@ void FixEvent::store_event(int timestep, int delta_clock) for (int i = 0; i < nlocal; i++) domain->unmap(x[i],image[i],xevent[i]); - event_timestep = timestep; - clock += delta_clock; - event_number++; +// printf("store_event %g %d %g %d \n", +// x[8][1],image[8],xevent[8][1],0); + +} + +/* ---------------------------------------------------------------------- + restore atom coords to quenched initial state + called prior to NEB calculation +------------------------------------------------------------------------- */ + +void FixEvent::restore_event() +{ + double **x = atom->x; + int *image = atom->image; + int nlocal = atom->nlocal; + +// printf("restore_event1 %g %d %g %d \n", +// x[8][1],image[8],xevent[8][1],0); + + for (int i = 0; i < nlocal; i++) { + x[i][0] = xevent[i][0]; + x[i][1] = xevent[i][1]; + x[i][2] = xevent[i][2]; + + // Since xevent is unwrapped coordinate, need to + // adjust image flags when remapping + + image[i] = (512 << 20) | (512 << 10) | 512; + domain->remap(x[i],image[i]); + // domain->remap(x[i]); + } + +// printf("restore_event2 %g %d %g %d \n", +// x[8][1],image[8],xevent[8][1],0); + } /* ---------------------------------------------------------------------- @@ -104,14 +132,20 @@ void FixEvent::store_event(int timestep, int delta_clock) void FixEvent::store_state() { double **x = atom->x; - double **f = atom->f; + double **v = atom->v; int *image = atom->image; int nlocal = atom->nlocal; +// printf("store_state %g %d %g %d \n", +// xold[8][1],imageold[8],x[8][1],image[8]); + 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]; + vold[i][0] = v[i][0]; + vold[i][1] = v[i][1]; + vold[i][2] = v[i][2]; imageold[i] = image[i]; } } @@ -124,13 +158,20 @@ void FixEvent::store_state() void FixEvent::restore_state() { double **x = atom->x; + double **v = atom->v; int *image = atom->image; int nlocal = atom->nlocal; +// printf("restore_state %g %d %g %d \n", +// xold[8][1],imageold[8],x[8][1],image[8]); + 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]; + v[i][0] = vold[i][0]; + v[i][1] = vold[i][1]; + v[i][2] = vold[i][2]; image[i] = imageold[i]; } } @@ -154,6 +195,7 @@ 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"); + vold = memory->grow_2d_double_array(vold,nmax,3,"event:vold"); imageold = (int *) memory->srealloc(imageold,nmax*sizeof(int),"event:imageold"); @@ -174,6 +216,9 @@ void FixEvent::copy_arrays(int i, int j) xold[j][0] = xold[i][0]; xold[j][1] = xold[i][1]; xold[j][2] = xold[i][2]; + vold[j][0] = vold[i][0]; + vold[j][1] = vold[i][1]; + vold[j][2] = vold[i][2]; imageold[j] = imageold[i]; } @@ -189,9 +234,12 @@ int FixEvent::pack_exchange(int i, double *buf) buf[3] = xold[i][0]; buf[4] = xold[i][1]; buf[5] = xold[i][2]; - buf[6] = imageold[i]; + buf[6] = vold[i][0]; + buf[7] = vold[i][1]; + buf[8] = vold[i][2]; + buf[9] = imageold[i]; - return 7; + return 10; } /* ---------------------------------------------------------------------- @@ -206,9 +254,12 @@ int FixEvent::unpack_exchange(int nlocal, double *buf) xold[nlocal][0] = buf[3]; xold[nlocal][1] = buf[4]; xold[nlocal][2] = buf[5]; - imageold[nlocal] = static_cast(buf[6]); + vold[nlocal][0] = buf[6]; + vold[nlocal][1] = buf[7]; + vold[nlocal][2] = buf[8]; + imageold[nlocal] = static_cast(buf[9]); - return 7; + return 10; } /* ---------------------------------------------------------------------- @@ -217,20 +268,6 @@ int FixEvent::unpack_exchange(int nlocal, double *buf) 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; - list[n++] = ncoincident; - - if (comm->me == 0) { - int size = n * sizeof(double); - fwrite(&size,sizeof(int),1,fp); - fwrite(list,sizeof(double),n,fp); - } } /* ---------------------------------------------------------------------- @@ -239,13 +276,4 @@ void FixEvent::write_restart(FILE *fp) 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++]); - ncoincident = static_cast (list[n++]); } diff --git a/src/REPLICA/fix_event.h b/src/REPLICA/fix_event.h index 2db729bc80..63ae46a06f 100644 --- a/src/REPLICA/fix_event.h +++ b/src/REPLICA/fix_event.h @@ -11,12 +11,6 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#ifdef FIX_CLASS - -FixStyle(EVENT,FixEvent) - -#else - #ifndef LMP_FIX_EVENT_H #define LMP_FIX_EVENT_H @@ -26,15 +20,9 @@ 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 - int ncoincident; // # of simultaneous events on different replicas FixEvent(class LAMMPS *, int, char **); - ~FixEvent(); + virtual ~FixEvent()=0; // Use destructor to make base class virtual int setmask(); double memory_usage(); @@ -42,22 +30,23 @@ class FixEvent : public Fix { void copy_arrays(int, int); int pack_exchange(int, double *); int unpack_exchange(int, double *); - void write_restart(FILE *); - void restart(char *); + virtual void write_restart(FILE *); + virtual void restart(char *); - // methods specific to FixEvent, invoked by PRD + // methods specific to FixEvent - void store_event(int, int); + virtual void store_event(); // base class stores quenched atoms + void restore_event(); // restore quenched atoms void store_state(); void restore_state(); 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 }; } #endif -#endif diff --git a/src/REPLICA/fix_event_prd.cpp b/src/REPLICA/fix_event_prd.cpp new file mode 100644 index 0000000000..d6d0e063b1 --- /dev/null +++ b/src/REPLICA/fix_event_prd.cpp @@ -0,0 +1,107 @@ +/* ---------------------------------------------------------------------- + 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_prd.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; + +/* ---------------------------------------------------------------------- */ + +FixEventPRD::FixEventPRD(LAMMPS *lmp, int narg, char **arg) : + FixEvent(lmp, narg, arg) +{ + if (narg != 3) error->all("Illegal fix event command"); + + restart_global = 1; + + event_number = 0; + event_timestep = update->ntimestep; + clock = 0; +} + +/* ---------------------------------------------------------------------- */ + +FixEventPRD::~FixEventPRD() +{ +} + +/* ---------------------------------------------------------------------- + save current atom coords as an event (via call to base class) + 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 FixEventPRD::store_event(int timestep, int delta_clock) +{ + FixEvent::store_event(); + event_timestep = timestep; + clock += delta_clock; + event_number++; +} + +/* ---------------------------------------------------------------------- + pack entire state of Fix into one write +------------------------------------------------------------------------- */ + +void FixEventPRD::write_restart(FILE *fp) +{ + int n = 0; + double list[6]; + list[n++] = event_number; + list[n++] = event_timestep; + list[n++] = clock; + list[n++] = replica_number; + list[n++] = correlated_event; + list[n++] = ncoincident; + + 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 FixEventPRD::restart(char *buf) +{ + int n = 0; + double *list = (double *) buf; + + event_number = static_cast (list[n++]); + printf("Event number restart = %d\n",event_number); + 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 new file mode 100644 index 0000000000..03b0cb6bb1 --- /dev/null +++ b/src/REPLICA/fix_event_prd.h @@ -0,0 +1,53 @@ +/* ---------------------------------------------------------------------- + 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 FIX_CLASS + +FixStyle(EVENT/PRD,FixEventPRD) + +#else + +#ifndef LMP_FIX_EVENT_PRD_H +#define LMP_FIX_EVENT_PRD_H + +#include "fix_event.h" + +namespace LAMMPS_NS { + +class FixEventPRD : public FixEvent { + 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 + int ncoincident; // # of simultaneous events on different replicas + + FixEventPRD(class LAMMPS *, int, char **); + ~FixEventPRD(); + + void write_restart(FILE *); + void restart(char *); + + // methods specific to FixEventPRD, invoked by PRD + + void store_event(int, int); + + private: + +}; + +} + +#endif +#endif diff --git a/src/REPLICA/fix_event_tad.cpp b/src/REPLICA/fix_event_tad.cpp new file mode 100644 index 0000000000..c027f4957c --- /dev/null +++ b/src/REPLICA/fix_event_tad.cpp @@ -0,0 +1,100 @@ +/* ---------------------------------------------------------------------- + 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_tad.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; + +/* ---------------------------------------------------------------------- */ + +FixEventTAD::FixEventTAD(LAMMPS *lmp, int narg, char **arg) : + FixEvent(lmp, narg, arg) +{ + if (narg != 3) error->all("Illegal fix event command"); + + restart_global = 1; + + event_number = 0; + event_timestep = update->ntimestep; + tlo = 0.0; + ebarrier = 0.0; +} + +/* ---------------------------------------------------------------------- */ + +FixEventTAD::~FixEventTAD() +{ +} + +/* ---------------------------------------------------------------------- + save current atom coords as an event (via call to base class) + called when an event occurs in some replica + set event_timestep = when event occurred +------------------------------------------------------------------------- */ + +void FixEventTAD::store_event(int timestep) +{ + FixEvent::store_event(); + event_timestep = timestep; +} + +/* ---------------------------------------------------------------------- + pack entire state of Fix into one write +------------------------------------------------------------------------- */ + +void FixEventTAD::write_restart(FILE *fp) +{ + int n = 0; + double list[4]; + list[n++] = event_number; + list[n++] = event_timestep; + list[n++] = tlo; + list[n++] = ebarrier; + + 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 FixEventTAD::restart(char *buf) +{ + int n = 0; + double *list = (double *) buf; + + event_number = static_cast (list[n++]); + event_timestep = static_cast (list[n++]); + tlo = static_cast (list[n++]); + ebarrier = list[n++]; +} + + diff --git a/src/REPLICA/fix_event_tad.h b/src/REPLICA/fix_event_tad.h new file mode 100644 index 0000000000..129b97b874 --- /dev/null +++ b/src/REPLICA/fix_event_tad.h @@ -0,0 +1,51 @@ +/* ---------------------------------------------------------------------- + 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 FIX_CLASS + +FixStyle(EVENT/TAD,FixEventTAD) + +#else + +#ifndef LMP_FIX_EVENT_TAD_H +#define LMP_FIX_EVENT_TAD_H + +#include "fix_event.h" + +namespace LAMMPS_NS { + +class FixEventTAD : public FixEvent { + public: + int event_number; // event counter + int event_timestep; // timestep of last event + double tlo; // event time at low temperature + double ebarrier; // energy barrier for this event + + FixEventTAD(class LAMMPS *, int, char **); + ~FixEventTAD(); + + void write_restart(FILE *); + void restart(char *); + + // methods specific to FixEventTAD, invoked by TAD + + void store_event(int); + + private: + +}; + +} + +#endif +#endif diff --git a/src/REPLICA/fix_neb.cpp b/src/REPLICA/fix_neb.cpp index 8e7cb00497..6d9074964f 100644 --- a/src/REPLICA/fix_neb.cpp +++ b/src/REPLICA/fix_neb.cpp @@ -132,7 +132,6 @@ void FixNEB::min_setup(int vflag) void FixNEB::min_post_force(int vflag) { MPI_Status status; - MPI_Request request; double vprev,vnext,vmax,vmin; double delx,dely,delz; double delta1[3],delta2[3]; @@ -142,21 +141,10 @@ void FixNEB::min_post_force(int vflag) veng = pe->compute_scalar(); - if (ireplica == 0) - MPI_Send(&veng,1,MPI_DOUBLE,procnext,0,uworld); - else if (ireplica == nreplica-1) { - MPI_Irecv(&vprev,1,MPI_DOUBLE,procprev,0,uworld,&request); - MPI_Wait(&request,&status); - } else + if (ireplica < nreplica-1) MPI_Sendrecv(&veng,1,MPI_DOUBLE,procnext,0, &vprev,1,MPI_DOUBLE,procprev,0,uworld,&status); - - if (ireplica == 0) { - MPI_Irecv(&vnext,1,MPI_DOUBLE,procnext,0,uworld,&request); - MPI_Wait(&request,&status); - } else if (ireplica == nreplica-1) - MPI_Send(&veng,1,MPI_DOUBLE,procprev,0,uworld); - else + if (ireplica > 0) MPI_Sendrecv(&veng,1,MPI_DOUBLE,procprev,0, &vnext,1,MPI_DOUBLE,procnext,0,uworld,&status); @@ -169,21 +157,10 @@ void FixNEB::min_post_force(int vflag) int nlocal = atom->nlocal; if (nlocal != nebatoms) error->one("Atom count changed in fix neb"); - if (ireplica == 0) - MPI_Send(x[0],3*nlocal,MPI_DOUBLE,procnext,0,uworld); - else if (ireplica == nreplica-1) { - MPI_Irecv(xprev[0],3*nlocal,MPI_DOUBLE,procprev,0,uworld,&request); - MPI_Wait(&request,&status); - } else + if (ireplica < nreplica-1) MPI_Sendrecv(x[0],3*nlocal,MPI_DOUBLE,procnext,0, xprev[0],3*nlocal,MPI_DOUBLE,procprev,0,uworld,&status); - - if (ireplica == 0) { - MPI_Irecv(xnext[0],3*nlocal,MPI_DOUBLE,procnext,0,uworld,&request); - MPI_Wait(&request,&status); - } else if (ireplica == nreplica-1) - MPI_Send(x[0],3*nlocal,MPI_DOUBLE,procprev,0,uworld); - else + if (ireplica > 0) MPI_Sendrecv(x[0],3*nlocal,MPI_DOUBLE,procprev,0, xnext[0],3*nlocal,MPI_DOUBLE,procnext,0,uworld,&status); diff --git a/src/REPLICA/neb.cpp b/src/REPLICA/neb.cpp index cc488e72fb..e44c95f7f0 100644 --- a/src/REPLICA/neb.cpp +++ b/src/REPLICA/neb.cpp @@ -40,6 +40,50 @@ using namespace LAMMPS_NS; NEB::NEB(LAMMPS *lmp) : Pointers(lmp) {} +/* ---------------------------------------------------------------------- + internal NEB constructor +------------------------------------------------------------------------- */ + +NEB::NEB(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in, + int n2steps_in, int nevery_in, double *buf_init, double *buf_final) + : Pointers(lmp) +{ + double delx,dely,delz; + + etol = etol_in; + ftol = ftol_in; + n1steps = n1steps_in; + n2steps = n2steps_in; + nevery = nevery_in; + + // replica info + + nreplica = universe->nworlds; + ireplica = universe->iworld; + me_universe = universe->me; + uworld = universe->uworld; + MPI_Comm_rank(world,&me); + + // generate linear interpolate replica + + double fraction = ireplica/(nreplica-1.0); + + double **x = atom->x; + int nlocal = atom->nlocal; + + int ii = 0; + for (int i = 0; i < nlocal; i++) { + delx = buf_final[ii] - buf_init[ii]; + dely = buf_final[ii+1] - buf_init[ii+1]; + delz = buf_final[ii+2] - buf_init[ii+2]; + domain->minimum_image(delx,dely,delz); + x[i][0] = buf_init[ii] + fraction*delx; + x[i][1] = buf_init[ii+1] + fraction*dely; + x[i][2] = buf_init[ii+2] + fraction*delz; + ii += 3; + } +} + /* ---------------------------------------------------------------------- */ NEB::~NEB() @@ -58,16 +102,18 @@ void NEB::command(int narg, char **arg) if (domain->box_exist == 0) error->all("NEB command before simulation box is defined"); - if (narg != 5) error->universe_all("Illegal NEB command"); + if (narg != 6) error->universe_all("Illegal NEB command"); - double ftol = atof(arg[0]); - int n1steps = atoi(arg[1]); - int n2steps = atoi(arg[2]); - int nevery = atoi(arg[3]); - char *infile = arg[4]; + etol = atof(arg[0]); + ftol = atof(arg[1]); + n1steps = atoi(arg[2]); + n2steps = atoi(arg[3]); + nevery = atoi(arg[4]); + infile = arg[5]; // error checks + if (etol < 0.0) error->all("Illegal NEB command"); if (ftol < 0.0) error->all("Illegal NEB command"); if (nevery == 0) error->universe_all("Illegal NEB command"); if (n1steps % nevery || n2steps % nevery) @@ -81,6 +127,22 @@ void NEB::command(int narg, char **arg) uworld = universe->uworld; MPI_Comm_rank(world,&me); + // read in file of final state atom coords and reset my coords + + readfile(infile); + + // run the NEB calculation + + run(); + +} + +/* ---------------------------------------------------------------------- + run NEB on multiple replicas +------------------------------------------------------------------------- */ + +void NEB::run() +{ // create MPI communicator for root proc from each world int color; @@ -104,13 +166,14 @@ void NEB::command(int narg, char **arg) if (ineb == modify->nfix) error->all("NEB requires use of fix neb"); fneb = (FixNEB *) modify->fix[ineb]; - all = memory->create_2d_double_array(nreplica,3,"neb:all"); + nall = 4; + all = memory->create_2d_double_array(nreplica,nall,"neb:all"); rdist = new double[nreplica]; // initialize LAMMPS update->whichflag = 2; - update->etol = 0.0; + update->etol = etol; update->ftol = ftol; update->multireplica = 1; @@ -119,11 +182,7 @@ void NEB::command(int narg, char **arg) if (update->minimize->searchflag) error->all("NEB requires damped dynamics minimizer"); - // read in file of final state atom coords and reset my coords - - readfile(infile); - - // setup regular NEB minimizaiton + // setup regular NEB minimization if (me_universe == 0 && universe->uscreen) fprintf(universe->uscreen,"Setting up regular NEB ...\n"); @@ -139,10 +198,12 @@ void NEB::command(int narg, char **arg) if (universe->uscreen) fprintf(universe->uscreen,"Step MaxReplicaForce MaxAtomForce " "GradV0 GradV1 GradVc " + "EBF EBR RDT " "RD1 PE1 RD2 PE2 ... RDN PEN\n"); if (universe->ulogfile) fprintf(universe->ulogfile,"Step MaxReplicaForce MaxAtomForce " "GradV0 GradV1 GradVc " + "EBF EBR RDT " "RD1 PE1 RD2 PE2 ... RDN PEN\n"); } print_status(); @@ -206,10 +267,12 @@ void NEB::command(int narg, char **arg) if (universe->uscreen) fprintf(universe->uscreen,"Step MaxReplicaForce MaxAtomForce " "GradV0 GradV1 GradVc " + "EBF EBR RDT " "RD1 PE1 RD2 PE2 ... RDN PEN\n"); if (universe->ulogfile) fprintf(universe->ulogfile,"Step MaxReplicaForce MaxAtomForce " "GradV0 GradV1 GradVc " + "EBF EBR RDT " "RD1 PE1 RD2 PE2 ... RDN PEN\n"); } print_status(); @@ -220,7 +283,6 @@ void NEB::command(int narg, char **arg) // damped dynamic min styles insure all replicas converge together timer->barrier_start(TIME_LOOP); - while (update->minimize->niter < n2steps) { update->minimize->run(nevery); print_status(); @@ -372,38 +434,39 @@ void NEB::print_status() double fmaxatom; MPI_Allreduce(&fnorminf,&fmaxatom,1,MPI_DOUBLE,MPI_MAX,roots); - double one[3]; + double one[4]; one[0] = fneb->veng; one[1] = fneb->plen; one[2] = fneb->nlen; + one[nall-1] = fneb->gradvnorm; + if (output->thermo->normflag) one[0] /= atom->natoms; - if (me == 0) MPI_Allgather(one,3,MPI_DOUBLE,&all[0][0],3,MPI_DOUBLE,roots); + if (me == 0) + MPI_Allgather(one,nall,MPI_DOUBLE,&all[0][0],nall,MPI_DOUBLE,roots); rdist[0] = 0.0; for (int i = 1; i < nreplica; i++) rdist[i] = rdist[i-1] + all[i][1]; double endpt = rdist[nreplica-1] = rdist[nreplica-2] + all[nreplica-2][2]; - if (endpt > 0.0) - for (int i = 1; i < nreplica; i++) - rdist[i] /= endpt; + for (int i = 1; i < nreplica; i++) + rdist[i] /= endpt; // look up GradV for the initial, final, and climbing replicas - // these should be identical to fnorm2 - // but to be safe take them straight from fix neb + // these are identical to fnorm2, but to be safe we + // take them straight from fix_neb double gradvnorm0, gradvnorm1, gradvnormc; int irep; irep = 0; - if (me_universe == irep) gradvnorm0 = fneb->gradvnorm; - MPI_Bcast(&gradvnorm0,1,MPI_DOUBLE,irep,uworld); + gradvnorm0 = all[irep][3]; irep = nreplica-1; - if (me_universe == irep) gradvnorm1 = fneb->gradvnorm; - MPI_Bcast(&gradvnorm1,1,MPI_DOUBLE,irep,uworld); + gradvnorm1 = all[irep][3]; irep = fneb->rclimber; if (irep > -1) { - if (me_universe == irep) gradvnormc = fneb->gradvnorm; - MPI_Bcast(&gradvnormc,1,MPI_DOUBLE,irep,uworld); + gradvnormc = all[irep][3]; + ebf = all[irep][0]-all[0][0]; + ebr = all[irep][0]-all[nreplica-1][0]; } else { double vmax = all[0][0]; int top = 0; @@ -413,8 +476,9 @@ void NEB::print_status() top = m; } irep = top; - if (me_universe == irep) gradvnormc = fneb->gradvnorm; - MPI_Bcast(&gradvnormc,1,MPI_DOUBLE,irep,uworld); + gradvnormc = all[irep][3]; + ebf = all[irep][0]-all[0][0]; + ebr = all[irep][0]-all[nreplica-1][0]; } if (me_universe == 0) { @@ -423,6 +487,7 @@ void NEB::print_status() fmaxreplica,fmaxatom); fprintf(universe->uscreen,"%g %g %g ", gradvnorm0,gradvnorm1,gradvnormc); + fprintf(universe->uscreen,"%g %g %g ",ebf,ebr,endpt); for (int i = 0; i < nreplica; i++) fprintf(universe->uscreen,"%g %g ",rdist[i],all[i][0]); fprintf(universe->uscreen,"\n"); @@ -432,6 +497,7 @@ void NEB::print_status() fmaxreplica,fmaxatom); fprintf(universe->ulogfile,"%g %g %g ", gradvnorm0,gradvnorm1,gradvnormc); + fprintf(universe->ulogfile,"%g %g %g ",ebf,ebr,endpt); for (int i = 0; i < nreplica; i++) fprintf(universe->ulogfile,"%g %g ",rdist[i],all[i][0]); fprintf(universe->ulogfile,"\n"); diff --git a/src/REPLICA/neb.h b/src/REPLICA/neb.h index 5d3aa46cc0..363396764c 100644 --- a/src/REPLICA/neb.h +++ b/src/REPLICA/neb.h @@ -28,8 +28,12 @@ namespace LAMMPS_NS { class NEB : protected Pointers { public: NEB(class LAMMPS *); + NEB(class LAMMPS *, double, double, int, int, int, double *, double *); ~NEB(); - void command(int, char **); + void command(int, char **); // process neb command + void run(); // run NEB + + double ebf,ebr; // forward and reverse energy barriers private: int me,me_universe; // my proc ID in world and universe @@ -38,9 +42,15 @@ class NEB : protected Pointers { MPI_Comm roots; // MPI comm with 1 root proc from each world FILE *fp; int compressed; + double etol; // energy tolerance convergence criterion + double ftol; // force tolerance convergence criterion + int n1steps, n2steps; // number of steps in stage 1 and 2 + int nevery; // output interval + char *infile; // name of file containing final state class FixNEB *fneb; - double **all; // PE,plen,nlen from each replica + int nall; // per-replica dimension of array all + double **all; // PE,plen,nlen,gradvnorm from each replica double *rdist; // normalize reaction distance, 0 to 1 void readfile(char *); diff --git a/src/REPLICA/prd.cpp b/src/REPLICA/prd.cpp index 7864eb28a9..de7f1be790 100644 --- a/src/REPLICA/prd.cpp +++ b/src/REPLICA/prd.cpp @@ -34,7 +34,7 @@ #include "modify.h" #include "compute.h" #include "fix.h" -#include "fix_event.h" +#include "fix_event_prd.h" #include "force.h" #include "pair.h" #include "random_park.h" @@ -162,13 +162,13 @@ void PRD::command(int narg, char **arg) args[1] = (char *) dist_setting; if (dist_setting) velocity->options(2,args); - // create FixEvent class to store event and pre-quench states + // create FixEventPRD class to store event and pre-quench states args[0] = (char *) "prd_event"; args[1] = (char *) "all"; - args[2] = (char *) "EVENT"; + args[2] = (char *) "EVENT/PRD"; modify->add_fix(3,args); - fix_event = (FixEvent *) modify->fix[modify->nfix-1]; + fix_event = (FixEventPRD *) modify->fix[modify->nfix-1]; // create Finish for timing output @@ -180,7 +180,7 @@ void PRD::command(int narg, char **arg) delete [] loop_setting; delete [] dist_setting; - // assign FixEvent to event-detection compute + // assign FixEventPRD to event-detection compute // necessary so it will know atom coords at last event int icompute = modify->find_compute(id_compute); diff --git a/src/REPLICA/prd.h b/src/REPLICA/prd.h index 916403d4ee..c060718de6 100644 --- a/src/REPLICA/prd.h +++ b/src/REPLICA/prd.h @@ -54,7 +54,7 @@ class PRD : protected Pointers { class RanPark *random_select; class RanMars *random_dephase; class Compute *compute_event; - class FixEvent *fix_event; + class FixEventPRD *fix_event; class Velocity *velocity; class Compute *temperature; class Finish *finish; diff --git a/src/REPLICA/tad.cpp b/src/REPLICA/tad.cpp new file mode 100644 index 0000000000..5f35ebab08 --- /dev/null +++ b/src/REPLICA/tad.cpp @@ -0,0 +1,987 @@ +// To do: +// Mysterious problem with // if (universe->iworld == 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Aidan Thompson (SNL) +------------------------------------------------------------------------- */ + +#include "mpi.h" +#include "math.h" +#include "stdlib.h" +#include "string.h" +#include "tad.h" +#include "universe.h" +#include "update.h" +#include "atom.h" +#include "domain.h" +#include "region.h" +#include "comm.h" +#include "velocity.h" +#include "integrate.h" +#include "min.h" +#include "neighbor.h" +#include "modify.h" +#include "neb.h" +#include "compute.h" +#include "fix.h" +#include "fix_event_tad.h" +#include "fix_store_state.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; + +/* ---------------------------------------------------------------------- */ + +TAD::TAD(LAMMPS *lmp) : Pointers(lmp) {} + +/* ---------------------------------------------------------------------- */ + +TAD::~TAD() { + memory->sfree(fix_event_list); + if (neb_logfilename != NULL) delete [] neb_logfilename; + delete [] min_style; + delete [] min_style_neb; +} + +/* ---------------------------------------------------------------------- + perform TAD simulation on root proc + other procs only used for NEB calcs +------------------------------------------------------------------------- */ + +void TAD::command(int narg, char **arg) +{ + + fix_event_list = NULL; + n_event_list = 0; + nmax_event_list = 0; + nmin_event_list = 10; + + // error checks + + if (domain->box_exist == 0) + error->all("tad command before simulation box is defined"); + if (universe->nworlds == 1) error->all("Cannot use TAD with a single replica for NEB"); + if (universe->nworlds != universe->nprocs) + error->all("Can only use TAD with 1-processor replicas for NEB"); + if (atom->sortfreq > 0) + error->all("Cannot use TAD with atom_modify sort enabled for NEB"); + if (atom->map_style == 0) + error->all("Cannot use TAD unless atom map exists for NEB"); + + if (narg < 7) error->universe_all("Illegal tad command"); + + nsteps = atoi(arg[0]); + t_event = atoi(arg[1]); + templo = atof(arg[2]); + temphi = atof(arg[3]); + delta_conf = atof(arg[4]); + tmax = atof(arg[5]); + + char id_compute[strlen(arg[6])+1]; + strcpy(id_compute,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 tad command"); + if (nsteps % t_event) + error->universe_all("TAD nsteps must be multiple of t_event"); + + if (delta_conf <= 0.0 || delta_conf >= 1.0) + error->universe_all("Invalid delta_conf in tad command"); + + if (tmax <= 0.0) + error->universe_all("Invalid tmax in tad command"); + + // deltconf = (ln(1/delta))/freq_min (timestep units) + + deltconf = -log(delta_conf)*tmax/update->dt; + + // local storage + + int me_universe = universe->me; + int nprocs_universe = universe->nprocs; + + MPI_Comm_rank(world,&me); + MPI_Comm_size(world,&nprocs); + + delta_beta = (1.0/templo - 1.0/temphi) / force->boltz; + ratio_beta = templo/temphi; + + // create FixEventTAD object to store last event + + int narg2 = 3; + char **args = new char*[narg2]; + args[0] = (char *) "tad_event"; + args[1] = (char *) "all"; + args[2] = (char *) "EVENT/TAD"; + modify->add_fix(narg2,args); + fix_event = (FixEventTAD *) modify->fix[modify->nfix-1]; + delete [] args; + + // create FixStoreState object to store revert state + + narg2 = 13; + args = new char*[narg2]; + args[0] = (char *) "tad_revert"; + args[1] = (char *) "all"; + args[2] = (char *) "store/state"; + args[3] = (char *) "0"; + args[4] = (char *) "x"; + args[5] = (char *) "y"; + args[6] = (char *) "z"; + args[7] = (char *) "ix"; + args[8] = (char *) "iy"; + args[9] = (char *) "iz"; + args[10] = (char *) "vx"; + args[11] = (char *) "vy"; + args[12] = (char *) "vz"; + modify->add_fix(narg2,args); + fix_revert = (FixStoreState *) modify->fix[modify->nfix-1]; + delete [] args; + + // create Finish for timing output + + finish = new Finish(lmp); + + // assign FixEventTAD 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 TAD"); + compute_event = modify->compute[icompute]; + compute_event->reset_extra_compute_fix("tad_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_universe == 0) + error->warning("Resetting reneighboring criteria during TAD"); + } + + neighbor->every = 1; + neighbor->delay = 0; + neighbor->dist_check = 1; + + // initialize TAD 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(); + + // set minimize style for quench + + narg2 = 1; + args = new char*[narg2]; + args[0] = min_style; + + update->create_minimize(narg2,args); + + delete [] args; + + // init minimizer settings and minimizer itself + + update->etol = etol; + update->ftol = ftol; + update->max_eval = maxeval; + + update->minimize->init(); + + // perform TAD simulation + + if (me_universe == 0 && universe->uscreen) + fprintf(universe->uscreen,"Setting up TAD ...\n"); + + if (me_universe == 0) { + if (universe->uscreen) + fprintf(universe->uscreen,"Step CPU Clock Event " + "\n"); + if (universe->ulogfile) + fprintf(universe->ulogfile,"Step CPU Clock Event " + "\n"); + } + + ulogfile_lammps = universe->ulogfile; + uscreen_lammps = universe->uscreen; + ulogfile_neb = NULL; + uscreen_neb = NULL; + if (me_universe == 0 && neb_logfilename) + ulogfile_neb = fopen(neb_logfilename,"w"); + + // store hot state and quenched event, only on replica 0 + + // need this line if quench() does only setup_minimal() + // update->minimize->setup(); + + // This should work with if uncommented, but does not + // if (universe->iworld == 0) { + fix_event->store_state(); + quench(); + timer->barrier_start(TIME_LOOP); + time_start = timer->array[TIME_LOOP]; + fix_event->store_event(update->ntimestep); + log_event(); + fix_event->restore_state(); + + // do full init/setup + + update->whichflag = 1; + lmp->init(); + update->integrate->setup(); + // } + + // main loop: look for events until out of time + // (1) dynamics, store state, quench, check event, restore state + // (2) if event, perform NEB, record in fix_event_list + // (3) if confident, pick earliest event + + nbuild = ndanger = 0; + time_neb = time_dynamics = time_quench = time_comm = time_output = 0.0; + + timer->barrier_start(TIME_LOOP); + time_start = timer->array[TIME_LOOP]; + + int confident_flag, event_flag; + + if (universe->iworld == 0) { + while (update->ntimestep < update->endstep) { + + // initialize list of possible events + + initialize_event_list(); + confident_flag = 0; + + while (update->ntimestep < update->endstep) { + event_flag = 0; + while (update->ntimestep < update->endstep) { + + dynamics(); + + + fix_event->store_state(); + + + quench(); + + + event_flag = check_event(); + MPI_Bcast(&event_flag,1,MPI_INT,0,universe->uworld); + + if (event_flag) break; + + // restore hot state + + fix_event->restore_state(); + + // store hot state in revert + + fix_revert->end_of_step(); + } + if (!event_flag) break; + + add_event(); + + perform_neb(n_event_list-1); + compute_tlo(n_event_list-1); + confident_flag = check_confidence(); + MPI_Bcast(&confident_flag,1,MPI_INT,0,universe->uworld); + if (confident_flag) break; + if (universe->iworld == 0) revert(); + } + if (!confident_flag) break; + + perform_event(event_first); + + // need to sync timestep with TAD + + MPI_Bcast(&(update->ntimestep),1,MPI_INT,0,universe->uworld); + + int restart_flag = 0; + if (output->restart_every && universe->iworld == 0) + if (fix_event->event_number % output->restart_every == 0) + restart_flag = 1; + + // full init/setup since are starting after event + + update->whichflag = 1; + lmp->init(); + update->integrate->setup(); + + // 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]; + } + } + + } else { + + while (update->ntimestep < update->endstep) { + confident_flag = 0; + while (update->ntimestep < update->endstep) { + event_flag = 0; + while (update->ntimestep < update->endstep) { + update->ntimestep += t_event; + MPI_Bcast(&event_flag,1,MPI_INT,0,universe->uworld); + + if (event_flag) break; + } + if (!event_flag) break; + perform_neb(-1); + MPI_Bcast(&confident_flag,1,MPI_INT,0,universe->uworld); + if (confident_flag) break; + } + if (!confident_flag) break; + + // need to sync timestep with TAD + + MPI_Bcast(&(update->ntimestep),1,MPI_INT,0,universe->uworld); + } + } + + // 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_neb; + 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); + } + + if (me_universe == 0) fclose(ulogfile_neb); + + finish->end(3); + + 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; + + delete finish; + modify->delete_fix("tad_event"); + modify->delete_fix("tad_revert"); + delete_event_list(); + + compute_event->reset_extra_compute_fix(NULL); +} + +/* ---------------------------------------------------------------------- + single short dynamics run +------------------------------------------------------------------------- */ + +void TAD::dynamics() +{ + update->whichflag = 1; + update->nsteps = t_event; + + lmp->init(); + update->integrate->setup(); + // this may be needed if don't do full init + //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 TAD::quench() +{ + int ntimestep_hold = update->ntimestep; + int endstep_hold = update->endstep; + + // need to change whichflag so that minimize->setup() calling + // modify->setup() will call fix->min_setup() + + update->whichflag = 2; + update->nsteps = maxiter; + update->endstep = update->laststep = update->firststep + maxiter; + + // full init works + + lmp->init(); + update->minimize->setup(); + + // partial init does 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(1); + + // reset timestep as if quench did not occur + // clear timestep storage from computes, since now invalid + + update->ntimestep = ntimestep_hold; + update->endstep = update->laststep = endstep_hold; + for (int i = 0; i < modify->ncompute; i++) + if (modify->compute[i]->timeflag) modify->compute[i]->clearstep(); +} + +/* ---------------------------------------------------------------------- + check for an event + return 0 if no event + return 1 if event +------------------------------------------------------------------------- */ + +int TAD::check_event() +{ + int flag; + + flag = 0; + if (compute_event->compute_scalar() > 0.0) flag = 1; + + return flag; +} + +/* ---------------------------------------------------------------------- + universe proc 0 prints event info +------------------------------------------------------------------------- */ + +void TAD::log_event() +{ + timer->array[TIME_LOOP] = time_start; + if (universe->me == 0) { + if (universe->uscreen) + fprintf(universe->uscreen,"%d %.3f %.3f %d\n", + fix_event->event_timestep, + timer->elapsed(TIME_LOOP), + fix_event->tlo, + fix_event->event_number); + if (universe->ulogfile) + fprintf(universe->ulogfile,"%d %.3f %.3f %d\n", + fix_event->event_timestep, + timer->elapsed(TIME_LOOP), + fix_event->tlo, + fix_event->event_number); + } + + // dump snapshot of quenched coords + // must reneighbor and compute forces before dumping + // 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]; + } + +} + +/* ---------------------------------------------------------------------- + parse optional parameters at end of TAD input line +------------------------------------------------------------------------- */ + +void TAD::options(int narg, char **arg) +{ + if (narg < 0) error->all("Illegal tad command"); + + // set defaults + + etol = 0.1; + ftol = 0.1; + maxiter = 40; + maxeval = 50; + + etol_neb = 0.01; + ftol_neb = 0.01; + n1steps_neb = 100; + n2steps_neb = 100; + nevery_neb = 10; + + min_style = new char[3]; + strcpy(min_style,"cg"); + min_style_neb = new char[9]; + strcpy(min_style_neb,"quickmin"); + neb_logfilename = NULL; + + int iarg = 0; + while (iarg < narg) { + if (strcmp(arg[iarg],"min") == 0) { + if (iarg+5 > narg) error->all("Illegal tad command"); + etol = atof(arg[iarg+1]); + ftol = atof(arg[iarg+2]); + maxiter = atoi(arg[iarg+3]); + maxeval = atoi(arg[iarg+4]); + if (maxiter < 0 || maxeval < 0 || + etol < 0.0 || ftol < 0.0 ) + error->all("Illegal tad command"); + iarg += 5; + + } else if (strcmp(arg[iarg],"neb") == 0) { + if (iarg+6 > narg) error->all("Illegal tad command"); + etol_neb = atof(arg[iarg+1]); + ftol_neb = atof(arg[iarg+2]); + n1steps_neb = atoi(arg[iarg+3]); + n2steps_neb = atoi(arg[iarg+4]); + nevery_neb = atoi(arg[iarg+5]); + if (etol_neb < 0.0 || ftol_neb < 0.0 || + n1steps_neb < 0 || n2steps_neb < 0 || + nevery_neb < 0) error->all("Illegal tad command"); + iarg += 6; + + } else if (strcmp(arg[iarg],"min_style") == 0) { + if (iarg+2 > narg) error->all("Illegal tad command"); + int n = strlen(arg[iarg+1]) + 1; + delete [] min_style; + min_style = new char[n]; + strcpy(min_style,arg[iarg+1]); + iarg += 2; + + } else if (strcmp(arg[iarg],"neb_style") == 0) { + if (iarg+2 > narg) error->all("Illegal tad command"); + int n = strlen(arg[iarg+1]) + 1; + delete [] min_style_neb; + min_style_neb = new char[n]; + strcpy(min_style_neb,arg[iarg+1]); + iarg += 2; + + } else if (strcmp(arg[iarg],"neb_log") == 0) { + delete [] neb_logfilename; + if (iarg+2 > narg) error->all("Illegal tad command"); + if (strcmp(arg[iarg+1],"none") == 0) neb_logfilename = NULL; + else { + int n = strlen(arg[iarg+1]) + 1; + neb_logfilename = new char[n]; + strcpy(neb_logfilename,arg[iarg+1]); + } + iarg += 2; + } else error->all("Illegal tad command"); + } +} + +/* ---------------------------------------------------------------------- + perform NEB calculation +------------------------------------------------------------------------- */ + +void TAD::perform_neb(int ievent) +{ + + double **x = atom->x; + int nlocal = atom->nlocal; + + double *buf_final = (double *) + memory->smalloc(3*nlocal*sizeof(double),"tad:buffinal"); + + // set system to quenched state of event ievent + + if (universe->iworld == 0) { + + fix_event_list[ievent]->restore_event(); + + int ii = 0; + for (int i = 0; i < nlocal; i++) { + buf_final[ii++] = x[i][0]; + buf_final[ii++] = x[i][1]; + buf_final[ii++] = x[i][2]; + } + } + + MPI_Bcast(buf_final,3*nlocal,MPI_DOUBLE,universe->root_proc[0], + universe->uworld); + + double *buf_init = (double *) + memory->smalloc(3*nlocal*sizeof(double),"tad:bufinit"); + + // set system to quenched state of fix_event + + if (universe->iworld == 0) { + + fix_event->restore_event(); + + int ii = 0; + for (int i = 0; i < nlocal; i++) { + buf_init[ii++] = x[i][0]; + buf_init[ii++] = x[i][1]; + buf_init[ii++] = x[i][2]; + } + } + + MPI_Bcast(buf_init,3*nlocal,MPI_DOUBLE,universe->root_proc[0], + universe->uworld); + + // create FixNEB object to support NEB + + int narg2 = 4; + char **args = new char*[narg2]; + args[0] = (char *) "neb"; + args[1] = (char *) "all"; + args[2] = (char *) "neb"; + char str[128]; + args[3] = str; + double kspring = 1.0; + sprintf(args[3],"%f",kspring); + modify->add_fix(narg2,args); + fix_neb = (Fix *) modify->fix[modify->nfix-1]; + delete [] args; + + // switch minimize style to quickmin for NEB + + narg2 = 1; + args = new char*[narg2]; + args[0] = min_style_neb; + + update->create_minimize(narg2,args); + + delete [] args; + + // create NEB object + + neb = new NEB(lmp,etol_neb,ftol_neb,n1steps_neb, + n2steps_neb,nevery_neb,buf_init,buf_final); + + // free up temporary arrays + + memory->sfree(buf_init); + memory->sfree(buf_final); + + // Run NEB + + double beginstep_hold = update->beginstep; + double endstep_hold = update->endstep; + double ntimestep_hold = update->ntimestep; + double nsteps_hold = update->nsteps; + + if (universe->me == 0) { + universe->ulogfile = ulogfile_neb; + universe->uscreen = uscreen_neb; + } + + // Had to bypass timer interface + // because timer->array is reset + // inside neb->run() + +// timer->barrier_start(TIME_LOOP); +// neb->run(); +// timer->barrier_stop(TIME_LOOP); +// time_neb += timer->array[TIME_LOOP]; + + MPI_Barrier(world); + double time_tmp = MPI_Wtime(); + neb->run(); + MPI_Barrier(world); + time_neb += MPI_Wtime() - time_tmp; + + if (universe->me == 0) { + universe->ulogfile = ulogfile_lammps; + universe->uscreen = uscreen_lammps; + } + + // Extract barrier energy from NEB + + if (universe->iworld == 0) + fix_event_list[ievent]->ebarrier = neb->ebf; + + update->beginstep = update->firststep = beginstep_hold; + update->endstep = update->laststep = endstep_hold; + update->ntimestep = ntimestep_hold; + update->nsteps = nsteps_hold; + + // switch minimize style back for quench + + narg2 = 1; + args = new char*[narg2]; + args[0] = min_style; + + update->create_minimize(narg2,args); + + update->etol = etol; + update->ftol = ftol; + + delete [] args; + + // Clean up + + modify->delete_fix("neb"); + delete neb; +} + +/* ---------------------------------------------------------------------- + check if confidence criterion for tstop is satisfied + return 0 if not satisfied + return 1 if satisfied +------------------------------------------------------------------------- */ + +int TAD::check_confidence() +{ + int flag; + + // update stopping time + + deltstop = deltconf*pow(deltfirst/deltconf, ratio_beta); + + flag = 0; + if (deltstop < update->ntimestep - fix_event->event_timestep) flag = 1; + + return flag; +} + +/* ---------------------------------------------------------------------- + reflect back in to starting state +------------------------------------------------------------------------- */ + +void TAD::revert() +{ + double **x = atom->x; + double **v = atom->v; + int *image = atom->image; + int nlocal = atom->nlocal; + + double **array_atom = fix_revert->array_atom; + + for (int i = 0; i < nlocal; i++) { + x[i][0] = array_atom[i][0]; + x[i][1] = array_atom[i][1]; + x[i][2] = array_atom[i][2]; + image[i] = ((int(array_atom[i][5]) + 512 & 1023) << 20) | + ((int(array_atom[i][4]) + 512 & 1023) << 10) | + (int(array_atom[i][3]) + 512 & 1023); + v[i][0] = -array_atom[i][6]; + v[i][1] = -array_atom[i][7]; + v[i][2] = -array_atom[i][8]; + } +} + +/* ---------------------------------------------------------------------- + Initialize list of possible events +------------------------------------------------------------------------- */ + +void TAD::initialize_event_list() { + + // First delete old events, if any + + delete_event_list(); + + // Create new list + + n_event_list = 0; + grow_event_list(nmin_event_list); +} + +/* ---------------------------------------------------------------------- + Delete list of possible events +------------------------------------------------------------------------- */ + +void TAD::delete_event_list() { + + for (int i = 0; i < n_event_list; i++) { + char str[128]; + sprintf(str,"tad_event_%d",i); + modify->delete_fix(str); + } + memory->sfree(fix_event_list); + fix_event_list = NULL; + n_event_list = 0; + nmax_event_list = 0; + +} + +/* ---------------------------------------------------------------------- + add event +------------------------------------------------------------------------- */ + +void TAD::add_event() +{ + + // create FixEventTAD object to store possible event + + int narg = 3; + char **args = new char*[narg]; + + char str[128]; + sprintf(str,"tad_event_%d",n_event_list); + + args[0] = str; + args[1] = (char *) "all"; + args[2] = (char *) "EVENT/TAD"; + modify->add_fix(narg,args); + + if (n_event_list == nmax_event_list) + grow_event_list(nmax_event_list+nmin_event_list); + n_event_list += 1; + int ievent = n_event_list-1; + fix_event_list[ievent] = (FixEventTAD *) modify->fix[modify->nfix-1]; + + // store quenched state for new event + + fix_event_list[ievent]->store_event(update->ntimestep); + + // store hot state for new event + + fix_event->restore_state(); + fix_event_list[ievent]->store_state(); + + // string clean-up + + delete [] args; + +} + +/* ---------------------------------------------------------------------- + compute cold time for event ievent +------------------------------------------------------------------------- */ + +void TAD::compute_tlo(int ievent) +{ + double deltlo,delthi,ebarrier; + + ebarrier = fix_event_list[ievent]->ebarrier; + delthi = fix_event_list[ievent]->event_timestep + - fix_event->event_timestep; + deltlo = delthi*exp(ebarrier*delta_beta); + fix_event_list[ievent]->tlo = fix_event->tlo + deltlo; + + // first-replica output about each event + + if (universe->me == 0) { + char str[128]; + double tfrac = 0.0; + if (ievent > 0) tfrac = delthi/deltstop; +// sprintf(str, +// "New event: ievent = %d eb = %g t_lo = %g t_hi = %g t_hi/t_stop = %g \n", +// ievent,ebarrier,deltlo,delthi,tfrac); +// error->warning(str); + + if (screen) + fprintf(screen, + "New event: t_hi = %d ievent = %d eb = %g dt_lo = %g dt_hi/t_stop = %g \n", + fix_event_list[ievent]->event_timestep, + ievent,ebarrier,deltlo,tfrac); + if (logfile) + fprintf(logfile, + "New event: t_hi = %d ievent = %d eb = %g dt_lo = %g dt_hi/t_stop = %g \n", + fix_event_list[ievent]->event_timestep, + ievent,ebarrier,deltlo,tfrac); + } + + // update first event + + if (ievent == 0) { + deltfirst = deltlo; + event_first = ievent; + } else if (deltlo < deltfirst) { + deltfirst = deltlo; + event_first = ievent; + } + +} + +/* ---------------------------------------------------------------------- + perform event +------------------------------------------------------------------------- */ + +void TAD::perform_event(int ievent) +{ + // reset timestep to that of event + + update->ntimestep = fix_event_list[ievent]->event_timestep; + + // Copy event to current event + + fix_event->tlo = fix_event_list[ievent]->tlo; + fix_event->event_number++; + fix_event->event_timestep = update->ntimestep; + fix_event_list[ievent]->restore_event(); + fix_event->store_event(fix_event_list[ievent]->event_timestep); + + // output stats and dump for quenched state + + log_event(); + + // load and store hot state + + fix_event_list[ievent]->restore_state(); + fix_event->store_state(); +} + +/* ---------------------------------------------------------------------- + Allocate list of pointers to events +------------------------------------------------------------------------- */ + +void TAD::grow_event_list(int nmax) { + if (nmax_event_list > nmax) return; + fix_event_list = (FixEventTAD **) memory->srealloc(fix_event_list,nmax*sizeof(FixEventTAD *),"tad:eventlist"); + nmax_event_list = nmax; +} + diff --git a/src/REPLICA/tad.h b/src/REPLICA/tad.h new file mode 100644 index 0000000000..8d7f2c3f8e --- /dev/null +++ b/src/REPLICA/tad.h @@ -0,0 +1,88 @@ +/* ---------------------------------------------------------------------- + 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 COMMAND_CLASS + +CommandStyle(tad,TAD) + +#else + +#ifndef LMP_TAD_H +#define LMP_TAD_H + +#include "pointers.h" + +namespace LAMMPS_NS { + +class TAD : protected Pointers { + public: + TAD(class LAMMPS *); + ~TAD(); + void command(int, char **); + + private: + int me,nprocs; + int nsteps,t_event; + double templo,temphi,delta_conf,tmax; + double etol,ftol,etol_neb,ftol_neb; + int maxiter,maxeval,n1steps_neb,n2steps_neb,nevery_neb; + char *min_style, *min_style_neb; + double delta_beta,ratio_beta; + double deltconf,deltstop,deltfirst; // Times since last event + int event_first; + + int neigh_every,neigh_delay,neigh_dist_check; + int nbuild,ndanger; + int quench_reneighbor; + + double time_dynamics,time_quench,time_neb,time_comm,time_output; + double time_start; + + class NEB *neb; // NEB object + class Fix *fix_neb; // FixNEB object + class Compute *compute_event; // compute to detect event + class FixEventTAD *fix_event; // current event/state + class FixStoreState *fix_revert; // revert state + FixEventTAD **fix_event_list; // list of possible events + int n_event_list; // number of events + int nmax_event_list; // allocated events + int nmin_event_list; // minimum allocation + + char *neb_logfilename; // filename for ulogfile_neb + FILE *uscreen_neb; // neb universe screen output + FILE *ulogfile_neb; // neb universe logfile + FILE *uscreen_lammps; // lammps universe screen output + FILE *ulogfile_lammps; // lammps universe logfile + + class Finish *finish; + + void dynamics(); + void quench(); + int check_event(); + int check_confidence(); + void perform_neb(int); + void log_event(); + void options(int, char **); + void revert(); + void add_event(); + void perform_event(int); + void compute_tlo(int); + void grow_event_list(int); + void initialize_event_list(); + void delete_event_list(); +}; + +} + +#endif +#endif