diff --git a/src/PRD/Install.sh b/src/REPLICA/Install.sh similarity index 64% rename from src/PRD/Install.sh rename to src/REPLICA/Install.sh index 7c6da632ad..bfacd7835c 100644 --- a/src/PRD/Install.sh +++ b/src/REPLICA/Install.sh @@ -4,20 +4,32 @@ if (test $1 = 1) then cp compute_event_displace.cpp .. cp fix_event.cpp .. + cp fix_neb.cpp .. + cp neb.cpp .. cp prd.cpp .. + cp temper.cpp .. cp compute_event_displace.h .. cp fix_event.h .. + cp fix_neb.h .. + cp neb.h .. cp prd.h .. + cp temper.h .. elif (test $1 = 0) then rm ../compute_event_displace.cpp rm ../fix_event.cpp + rm ../fix_neb.cpp + rm ../neb.cpp rm ../prd.cpp + rm ../temper.cpp rm ../compute_event_displace.h rm ../fix_event.h + rm ../fix_neb.h + rm ../neb.h rm ../prd.h + rm ../temper.h fi diff --git a/src/PRD/compute_event_displace.cpp b/src/REPLICA/compute_event_displace.cpp similarity index 100% rename from src/PRD/compute_event_displace.cpp rename to src/REPLICA/compute_event_displace.cpp diff --git a/src/PRD/compute_event_displace.h b/src/REPLICA/compute_event_displace.h similarity index 100% rename from src/PRD/compute_event_displace.h rename to src/REPLICA/compute_event_displace.h diff --git a/src/PRD/fix_event.cpp b/src/REPLICA/fix_event.cpp similarity index 100% rename from src/PRD/fix_event.cpp rename to src/REPLICA/fix_event.cpp diff --git a/src/PRD/fix_event.h b/src/REPLICA/fix_event.h similarity index 100% rename from src/PRD/fix_event.h rename to src/REPLICA/fix_event.h diff --git a/src/REPLICA/fix_neb.cpp b/src/REPLICA/fix_neb.cpp new file mode 100644 index 0000000000..fbe89e188d --- /dev/null +++ b/src/REPLICA/fix_neb.cpp @@ -0,0 +1,299 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "mpi.h" +#include "math.h" +#include "stdlib.h" +#include "string.h" +#include "fix_neb.h" +#include "universe.h" +#include "update.h" +#include "domain.h" +#include "modify.h" +#include "compute.h" +#include "atom.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define MIN(A,B) ((A) < (B)) ? (A) : (B) +#define MAX(A,B) ((A) > (B)) ? (A) : (B) + +/* ---------------------------------------------------------------------- */ + +FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + if (narg != 4) error->all("Illegal fix neb command"); + + kspring = atof(arg[3]); + if (kspring <= 0.0) error->all("Illegal fix neb command"); + + // nreplica = number of partitions + // ireplica = which world I am in universe + // procprev,procnext = root proc in adjacent replicas + + nreplica = universe->nworlds; + ireplica = universe->iworld; + if (nreplica == 1) + error->all("Cannot use fix neb without multiple replicas"); + if (nreplica != universe->nprocs) + error->all("Can only use NEB with 1-processor replicas"); + + if (ireplica > 0) procprev = universe->root_proc[ireplica-1]; + else procprev = -1; + if (ireplica < nreplica-1) procnext = universe->root_proc[ireplica+1]; + else procnext = -1; + uworld = universe->uworld; + + // create a new compute pe style + // id = fix-ID + pe, compute group = all + + int n = strlen(id) + 4; + id_pe = new char[n]; + strcpy(id_pe,id); + strcat(id_pe,"_pe"); + + char **newarg = new char*[3]; + newarg[0] = id_pe; + newarg[1] = (char *) "all"; + newarg[2] = (char *) "pe"; + modify->add_compute(3,newarg); + delete [] newarg; + + xprev = xnext = tangent = NULL; +} + +/* ---------------------------------------------------------------------- */ + +FixNEB::~FixNEB() +{ + modify->delete_compute(id_pe); + delete [] id_pe; + + memory->destroy_2d_double_array(xprev); + memory->destroy_2d_double_array(xnext); + memory->destroy_2d_double_array(tangent); +} + +/* ---------------------------------------------------------------------- */ + +int FixNEB::setmask() +{ + int mask = 0; + mask |= MIN_POST_FORCE; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixNEB::init() +{ + int icompute = modify->find_compute(id_pe); + if (icompute < 0) + error->all("Potential energy ID for fix neb does not exist"); + pe = modify->compute[icompute]; + + // turn off climbing mode, NEB command turns it on after init() + + rclimber = -1; + + // setup xprev and xnext arrays + + memory->destroy_2d_double_array(xprev); + memory->destroy_2d_double_array(xnext); + memory->destroy_2d_double_array(tangent); + natoms = atom->nlocal; + xprev = memory->create_2d_double_array(natoms,3,"neb:xprev"); + xnext = memory->create_2d_double_array(natoms,3,"neb:xnext"); + tangent = memory->create_2d_double_array(natoms,3,"neb:tangent"); +} + +/* ---------------------------------------------------------------------- */ + +void FixNEB::min_setup(int vflag) +{ + min_post_force(vflag); + + // trigger potential energy computation on next timestep + + pe->addstep(update->ntimestep+1); +} + +/* ---------------------------------------------------------------------- */ + +void FixNEB::min_post_force(int vflag) +{ + MPI_Status status; + double vprev,vnext,vmax,vmin; + double delx,dely,delz; + double delta1[3],delta2[3]; + + // veng = PE of this replica + // vprev,vnext = PEs of adjacent replicas + + veng = pe->compute_scalar(); + + if (ireplica < nreplica-1) + MPI_Sendrecv(&veng,1,MPI_DOUBLE,procnext,0, + &vprev,1,MPI_DOUBLE,procprev,0,uworld,&status); + if (ireplica > 0) + MPI_Sendrecv(&veng,1,MPI_DOUBLE,procprev,0, + &vnext,1,MPI_DOUBLE,procnext,0,uworld,&status); + + // xprev,xnext = atom coords of adjacent replicas + // assume order of atoms in all replicas is the same + // check that number of atoms hasn't changed + + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (nlocal != natoms) error->one("Atom count changed in fix neb"); + + 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_Sendrecv(x[0],3*nlocal,MPI_DOUBLE,procprev,0, + xnext[0],3*nlocal,MPI_DOUBLE,procnext,0,uworld,&status); + + // trigger potential energy computation on next timestep + + pe->addstep(update->ntimestep+1); + + // if this is first or last replica, no change to forces, just return + + if (ireplica == 0 || ireplica == nreplica-1) { + plen = nlen = 0.0; + return; + } + + // tangent = unit tangent vector in 3N space + // based on delta vectors between atoms and their images in adjacent replicas + // use one or two delta vecs to compute tangent, + // depending on relative PEs of 3 replicas + // see Henkelman & Jonsson 2000 paper, eqs 8-11 + + if (vnext > veng && veng > vprev) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + tangent[i][0] = xnext[i][0] - x[i][0]; + tangent[i][1] = xnext[i][1] - x[i][1]; + tangent[i][2] = xnext[i][2] - x[i][2]; + domain->minimum_image(tangent[i]); + } + } else if (vnext < veng && veng < vprev) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + tangent[i][0] = x[i][0] - xprev[i][0]; + tangent[i][1] = x[i][1] - xprev[i][1]; + tangent[i][2] = x[i][2] - xprev[i][2]; + domain->minimum_image(tangent[i]); + } + } else { + vmax = MAX(fabs(vnext-veng),fabs(vprev-veng)); + vmin = MIN(fabs(vnext-veng),fabs(vprev-veng)); + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + delta1[0] = xnext[i][0] - x[i][0]; + delta1[1] = xnext[i][1] - x[i][1]; + delta1[2] = xnext[i][2] - x[i][2]; + domain->minimum_image(delta1); + delta2[0] = x[i][0] - xprev[i][0]; + delta2[1] = x[i][1] - xprev[i][1]; + delta2[2] = x[i][2] - xprev[i][2]; + domain->minimum_image(delta2); + if (vnext > vprev) { + tangent[i][0] = vmax*delta1[0] + vmin*delta2[0]; + tangent[i][1] = vmax*delta1[1] + vmin*delta2[1]; + tangent[i][2] = vmax*delta1[2] + vmin*delta2[2]; + } else { + tangent[i][0] = vmin*delta1[0] + vmax*delta2[0]; + tangent[i][1] = vmin*delta1[1] + vmax*delta2[1]; + tangent[i][2] = vmin*delta1[2] + vmax*delta2[2]; + } + } + } + + // tlen,plen,nlen = lengths of tangent, prev, next vectors + + double tlen = 0.0; + plen = 0.0; + nlen = 0.0; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + tlen += tangent[i][0]*tangent[i][0] + tangent[i][1]*tangent[i][1] + + tangent[i][2]*tangent[i][2]; + + delx = x[i][0] - xprev[i][0]; + dely = x[i][1] - xprev[i][1]; + delz = x[i][2] - xprev[i][2]; + domain->minimum_image(delx,dely,delz); + plen += delx*delx + dely*dely + delz*delz; + + delx = xnext[i][0] - x[i][0]; + dely = xnext[i][1] - x[i][1]; + delz = xnext[i][2] - x[i][2]; + domain->minimum_image(delx,dely,delz); + nlen += delx*delx + dely*dely + delz*delz; + } + + tlen = sqrt(tlen); + plen = sqrt(plen); + nlen = sqrt(nlen); + + // normalize tangent vector + + if (tlen > 0.0) { + double tleninv = 1.0/tlen; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + tangent[i][0] *= tleninv; + tangent[i][1] *= tleninv; + tangent[i][2] *= tleninv; + } + } + + // reset force on each atom in this replica + // regular NEB for all replicas except rclimber does hill-climbing NEB + // currently have F = -Grad(V) = -Grad(V)_perp - Grad(V)_parallel + // want F = -Grad(V)_perp + Fspring for regular NEB + // thus Fdelta = Grad(V)_parallel + Fspring for regular NEB + // want F = -Grad(V) + 2 Grad(V)_parallel for hill-climbing NEB + // thus Fdelta = 2 Grad(V)_parallel for hill-climbing NEB + // Grad(V)_parallel = (Grad(V) . utan) * utangent = -(F . utan) * utangent + // Fspring = k (nlen - plen) * utangent + // see Henkelman & Jonsson 2000 paper, eqs 3,4,12 + // see Henkelman & Jonsson 2000a paper, eq 5 + + double **f = atom->f; + + double dot = 0.0; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + dot += f[i][0]*tangent[i][0] + f[i][1]*tangent[i][1] + + f[i][2]*tangent[i][2]; + + double prefactor; + if (ireplica == rclimber) prefactor = -2.0*dot; + else prefactor = -dot + kspring*(nlen-plen); + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + f[i][0] += prefactor*tangent[i][0]; + f[i][1] += prefactor*tangent[i][1]; + f[i][2] += prefactor*tangent[i][2]; + } +} diff --git a/src/REPLICA/fix_neb.h b/src/REPLICA/fix_neb.h new file mode 100644 index 0000000000..7e7d97c635 --- /dev/null +++ b/src/REPLICA/fix_neb.h @@ -0,0 +1,56 @@ +/* ---------------------------------------------------------------------- + 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(neb,FixNEB) + +#else + +#ifndef LMP_FIX_NEB_H +#define LMP_FIX_NEB_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixNEB : public Fix { + public: + double veng,plen,nlen; + int rclimber; + + FixNEB(class LAMMPS *, int, char **); + ~FixNEB(); + int setmask(); + void init(); + void min_setup(int); + void min_post_force(int); + + private: + double kspring; + int ireplica,nreplica; + int procnext,procprev; + MPI_Comm uworld; + + char *id_pe; + class Compute *pe; + + int natoms; + double **xprev,**xnext; + double **tangent; +}; + +} + +#endif +#endif diff --git a/src/REPLICA/neb.cpp b/src/REPLICA/neb.cpp new file mode 100644 index 0000000000..49e8c2b0e7 --- /dev/null +++ b/src/REPLICA/neb.cpp @@ -0,0 +1,389 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "mpi.h" +#include "math.h" +#include "stdlib.h" +#include "string.h" +#include "neb.h" +#include "universe.h" +#include "atom.h" +#include "update.h" +#include "domain.h" +#include "min.h" +#include "modify.h" +#include "fix.h" +#include "fix_neb.h" +#include "output.h" +#include "thermo.h" +#include "finish.h" +#include "timer.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define CHUNK 1000 +#define MAXLINE 256 + +/* ---------------------------------------------------------------------- */ + +NEB::NEB(LAMMPS *lmp) : Pointers(lmp) {} + +/* ---------------------------------------------------------------------- */ + +NEB::~NEB() +{ + MPI_Comm_free(&roots); + memory->destroy_2d_double_array(all); + delete [] rdist; +} + +/* ---------------------------------------------------------------------- + perform NEB on multiple replicas +------------------------------------------------------------------------- */ + +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"); + + double ftol = atof(arg[0]); + int n1steps = atoi(arg[1]); + int n2steps = atoi(arg[2]); + int nevery = atoi(arg[3]); + char *infile = arg[4]; + + // error checks + + if (ftol < 0.0) error->all("Illegal NEB command"); + if (nevery == 0) error->universe_all("Illegal NEB command"); + if (n1steps % nevery || n2steps % nevery) + error->universe_all("Illegal NEB command"); + + // replica info + + nreplica = universe->nworlds; + ireplica = universe->iworld; + me_universe = universe->me; + uworld = universe->uworld; + MPI_Comm_rank(world,&me); + + // create MPI communicator for root proc from each world + + int color; + if (me == 0) color = 0; + else color = 1; + MPI_Comm_split(uworld,color,0,&roots); + + // error checks + + if (nreplica == 1) error->all("Cannot use NEB without multiple replicas"); + if (nreplica != universe->nprocs) + error->all("Can only use NEB with 1-processor replicas"); + if (atom->sortfreq > 0) + error->all("Cannot use NEB with atom_modify sort enabled"); + if (atom->map_style == 0) + error->all("Cannot use NEB unless atom map exists"); + + int ineb,idamp; + for (ineb = 0; ineb < modify->nfix; ineb++) + if (strcmp(modify->fix[ineb]->style,"neb") == 0) break; + 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"); + rdist = new double[nreplica]; + + // initialize LAMMPS + + update->whichflag = 2; + update->etol = 0.0; + update->ftol = ftol; + + lmp->init(); + + 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 + + if (me_universe == 0 && universe->uscreen) + fprintf(universe->uscreen,"Setting up regular NEB ...\n"); + + update->beginstep = update->firststep = update->ntimestep; + update->endstep = update->laststep = update->firststep + n1steps; + update->nsteps = n1steps; + update->max_eval = n1steps; + + update->minimize->setup(); + + if (me_universe == 0) { + if (universe->uscreen) + fprintf(universe->uscreen,"Step RD1 PE1 RD2 PE2 ... RDN PEN\n"); + if (universe->ulogfile) + fprintf(universe->ulogfile,"Step RD1 PE1 RD2 PE2 ... RDN PEN\n"); + } + print_status(); + + // perform regular NEB for n1steps or until all replicas converged + // retrieve PE values from fix NEB and print every nevery iterations + // break induced by MPI_Allreduce() and MPI_Bcast() + // only if all replicas have stop_condition > 0 + + int flag,flagall; + + timer->barrier_start(TIME_LOOP); + + while (update->minimize->niter < n1steps) { + update->minimize->run(nevery,1); + print_status(); + if (me == 0) { + flag = update->minimize->stop_condition; + MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MIN,roots); + } + MPI_Bcast(&flagall,1,MPI_INT,0,world); + if (flagall) break; + } + + timer->barrier_stop(TIME_LOOP); + + update->minimize->cleanup(); + + Finish finish(lmp); + finish.end(1); + + // switch fix NEB to climbing mode + + double vmax = all[0][0]; + int top = 0; + for (int m = 1; m < nreplica; m++) + if (vmax < all[m][0]) { + vmax = all[m][0]; + top = m; + } + + // setup climbing NEB minimization + // must reinitialize minimizer so it re-creates its fix MINIMIZE + + if (me_universe == 0 && universe->uscreen) + fprintf(universe->uscreen,"Setting up climbing NEB ...\n"); + + update->beginstep = update->firststep = update->ntimestep; + update->endstep = update->laststep = update->firststep + n2steps; + update->nsteps = n2steps; + update->max_eval = n2steps; + + update->minimize->init(); + fneb->rclimber = top; + update->minimize->setup(); + + if (me_universe == 0) { + if (universe->uscreen) + fprintf(universe->uscreen,"Step RD1 PE1 RD2 PE2 ... RDN PEN\n"); + if (universe->ulogfile) + fprintf(universe->ulogfile,"Step RD1 PE1 RD2 PE2 ... RDN PEN\n"); + } + print_status(); + + // perform climbing NEB for n2steps or until all replicas converged + // retrieve PE values from fix NEB and print every nevery iterations + // break induced by MPI_Allreduce() and MPI_Bcast() + // only if all replicas have stop_condition > 0 + + timer->barrier_start(TIME_LOOP); + + while (update->minimize->niter < n2steps) { + update->minimize->run(nevery,1); + print_status(); + if (me == 0) { + flag = update->minimize->stop_condition; + MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MIN,roots); + } + MPI_Bcast(&flagall,1,MPI_INT,0,world); + if (flagall) break; + } + + timer->barrier_stop(TIME_LOOP); + + update->minimize->cleanup(); + + finish.end(1); + + update->whichflag = 0; + update->firststep = update->laststep = 0; + update->beginstep = update->endstep = 0; +} + +/* ---------------------------------------------------------------------- + read target coordinates from file, store with appropriate atom + adjust coords of each atom based on ireplica + new coord = replica fraction between current and final state +------------------------------------------------------------------------- */ + +void NEB::readfile(char *file) +{ + if (me_universe == 0) { + if (screen) fprintf(screen,"Reading NEB target file %s ...\n",file); + open(file); + } + + double fraction = ireplica/(nreplica-1.0); + + double **x = atom->x; + int *image = atom->image; + int nlocal = atom->nlocal; + + char *buffer = new char[CHUNK*MAXLINE]; + char *ptr,*next,*bufptr; + int i,m,nlines,tag; + double xx,yy,zz,delx,dely,delz; + + int firstline = 1; + int ncount = 0; + int eof = 0; + + while (!eof) { + if (me_universe == 0) { + m = 0; + for (nlines = 0; nlines < CHUNK; nlines++) { + ptr = fgets(&buffer[m],MAXLINE,fp); + if (ptr == NULL) break; + m += strlen(&buffer[m]); + } + if (ptr == NULL) eof = 1; + buffer[m++] = '\n'; + } + + MPI_Bcast(&eof,1,MPI_INT,0,uworld); + MPI_Bcast(&nlines,1,MPI_INT,0,uworld); + MPI_Bcast(&m,1,MPI_INT,0,uworld); + MPI_Bcast(buffer,m,MPI_CHAR,0,uworld); + + bufptr = buffer; + for (i = 0; i < nlines; i++) { + next = strchr(bufptr,'\n'); + *next = '\0'; + + if (firstline) { + if (atom->count_words(bufptr) == 4) firstline = 0; + else error->all("Incorrect format in NEB target file"); + } + + sscanf(bufptr,"%d %lg %lg %lg",&tag,&xx,&yy,&zz); + + // adjust atom coord based on replica fraction + // ignore image flags of final x + // new x is displacement from old x + // if final x is across periodic boundary: + // new x may be outside box + // will be remapped back into box when simulation starts + // its image flags will be adjusted appropriately + + m = atom->map(tag); + if (m >= 0 && m < nlocal) { + delx = xx - x[m][0]; + dely = yy - x[m][1]; + delz = zz - x[m][2]; + domain->minimum_image(delx,dely,delz); + x[m][0] += fraction*delx; + x[m][1] += fraction*dely; + x[m][2] += fraction*delz; + ncount++; + } + + bufptr = next + 1; + } + } + + // clean up + + delete [] buffer; + + if (me_universe == 0) { + if (compressed) pclose(fp); + else fclose(fp); + } +} + +/* ---------------------------------------------------------------------- + universe proc 0 opens NEB data file + test if gzipped +------------------------------------------------------------------------- */ + +void NEB::open(char *file) +{ + compressed = 0; + char *suffix = file + strlen(file) - 3; + if (suffix > file && strcmp(suffix,".gz") == 0) compressed = 1; + if (!compressed) fp = fopen(file,"r"); + else { +#ifdef LAMMPS_GZIP + char gunzip[128]; + sprintf(gunzip,"gunzip -c %s",file); + fp = popen(gunzip,"r"); +#else + error->one("Cannot open gzipped file"); +#endif + } + + if (fp == NULL) { + char str[128]; + sprintf(str,"Cannot open file %s",file); + error->one(str); + } +} + +/* ---------------------------------------------------------------------- + query fix NEB for PE of each replica + proc 0 prints current NEB status +------------------------------------------------------------------------- */ + +void NEB::print_status() +{ + double one[3]; + one[0] = fneb->veng; + one[1] = fneb->plen; + one[2] = fneb->nlen; + if (output->thermo->normflag) one[0] /= atom->natoms; + if (me == 0) MPI_Allgather(one,3,MPI_DOUBLE,&all[0][0],3,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]; + for (int i = 1; i < nreplica; i++) + rdist[i] /= endpt; + + if (me_universe == 0) { + if (universe->uscreen) { + fprintf(universe->uscreen,"%d ",update->ntimestep); + for (int i = 0; i < nreplica; i++) + fprintf(universe->uscreen,"%g %g ",rdist[i],all[i][0]); + fprintf(universe->uscreen,"\n"); + } + if (universe->ulogfile) { + fprintf(universe->ulogfile,"%d ",update->ntimestep); + for (int i = 0; i < nreplica; i++) + fprintf(universe->ulogfile,"%g %g ",rdist[i],all[i][0]); + fprintf(universe->ulogfile,"\n"); + fflush(universe->ulogfile); + } + } +} diff --git a/src/REPLICA/neb.h b/src/REPLICA/neb.h new file mode 100644 index 0000000000..5d3aa46cc0 --- /dev/null +++ b/src/REPLICA/neb.h @@ -0,0 +1,54 @@ +/* ---------------------------------------------------------------------- + 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(neb,NEB) + +#else + +#ifndef LMP_NEB_H +#define LMP_NEB_H + +#include "stdio.h" +#include "pointers.h" + +namespace LAMMPS_NS { + +class NEB : protected Pointers { + public: + NEB(class LAMMPS *); + ~NEB(); + void command(int, char **); + + private: + int me,me_universe; // my proc ID in world and universe + int ireplica,nreplica; + MPI_Comm uworld; + MPI_Comm roots; // MPI comm with 1 root proc from each world + FILE *fp; + int compressed; + + class FixNEB *fneb; + double **all; // PE,plen,nlen from each replica + double *rdist; // normalize reaction distance, 0 to 1 + + void readfile(char *); + void open(char *); + void print_status(); +}; + +} + +#endif +#endif diff --git a/src/PRD/prd.cpp b/src/REPLICA/prd.cpp similarity index 100% rename from src/PRD/prd.cpp rename to src/REPLICA/prd.cpp diff --git a/src/PRD/prd.h b/src/REPLICA/prd.h similarity index 100% rename from src/PRD/prd.h rename to src/REPLICA/prd.h diff --git a/src/REPLICA/temper.cpp b/src/REPLICA/temper.cpp new file mode 100644 index 0000000000..592e012f3e --- /dev/null +++ b/src/REPLICA/temper.cpp @@ -0,0 +1,357 @@ +/* ---------------------------------------------------------------------- + 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: Mark Sears (SNL) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdlib.h" +#include "string.h" +#include "temper.h" +#include "universe.h" +#include "domain.h" +#include "atom.h" +#include "update.h" +#include "integrate.h" +#include "modify.h" +#include "compute.h" +#include "force.h" +#include "output.h" +#include "thermo.h" +#include "fix.h" +#include "random_park.h" +#include "finish.h" +#include "timer.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +// #define TEMPER_DEBUG 1 + +/* ---------------------------------------------------------------------- */ + +Temper::Temper(LAMMPS *lmp) : Pointers(lmp) {} + +/* ---------------------------------------------------------------------- */ + +Temper::~Temper() +{ + MPI_Comm_free(&roots); + if (ranswap) delete ranswap; + delete ranboltz; + delete [] set_temp; + delete [] temp2world; + delete [] world2temp; + delete [] world2root; +} + +/* ---------------------------------------------------------------------- + perform tempering with inter-world swaps +------------------------------------------------------------------------- */ + +void Temper::command(int narg, char **arg) +{ + if (universe->nworlds == 1) + error->all("Must have more than one processor partition to temper"); + if (domain->box_exist == 0) + error->all("Temper command before simulation box is defined"); + if (narg != 6 && narg != 7) error->universe_all("Illegal temper command"); + + int nsteps = atoi(arg[0]); + nevery = atoi(arg[1]); + double temp = atof(arg[2]); + + for (whichfix = 0; whichfix < modify->nfix; whichfix++) + if (strcmp(arg[3],modify->fix[whichfix]->id) == 0) break; + if (whichfix == modify->nfix) + error->universe_all("Tempering fix ID is not defined"); + + seed_swap = atoi(arg[4]); + seed_boltz = atoi(arg[5]); + + my_set_temp = universe->iworld; + if (narg == 7) my_set_temp = atoi(arg[6]); + + // swap frequency must evenly divide total # of timesteps + + if (nevery == 0) error->universe_all("Invalid frequency in temper command"); + nswaps = nsteps/nevery; + if (nswaps*nevery != nsteps) + error->universe_all("Non integer # of swaps in temper command"); + + // fix style must be appropriate for temperature control + + if ((strcmp(modify->fix[whichfix]->style,"nvt") != 0) && + (strcmp(modify->fix[whichfix]->style,"langevin") != 0) && + (strcmp(modify->fix[whichfix]->style,"temp/berendsen") != 0) && + (strcmp(modify->fix[whichfix]->style,"temp/rescale") != 0)) + error->universe_all("Tempering temperature fix is not valid"); + + // setup for long tempering run + + update->whichflag = 1; + update->nsteps = nsteps; + update->beginstep = update->firststep = update->ntimestep; + update->endstep = update->laststep = update->firststep + nsteps; + + lmp->init(); + + // local storage + + me_universe = universe->me; + MPI_Comm_rank(world,&me); + nworlds = universe->nworlds; + iworld = universe->iworld; + boltz = force->boltz; + + // pe_compute = ptr to thermo_pe compute + // notify compute it will be called at first swap + + int id = modify->find_compute("thermo_pe"); + if (id < 0) error->all("Tempering could not find thermo_pe compute"); + Compute *pe_compute = modify->compute[id]; + pe_compute->addstep(update->ntimestep + nevery); + + // create MPI communicator for root proc from each world + + int color; + if (me == 0) color = 0; + else color = 1; + MPI_Comm_split(universe->uworld,color,0,&roots); + + // RNGs for swaps and Boltzmann test + // warm up Boltzmann RNG + + if (seed_swap) ranswap = new RanPark(lmp,seed_swap); + else ranswap = NULL; + ranboltz = new RanPark(lmp,seed_boltz + me_universe); + for (int i = 0; i < 100; i++) ranboltz->uniform(); + + // world2root[i] = global proc that is root proc of world i + + world2root = new int[nworlds]; + if (me == 0) + MPI_Allgather(&me_universe,1,MPI_INT,world2root,1,MPI_INT,roots); + MPI_Bcast(world2root,nworlds,MPI_INT,0,world); + + // create static list of set temperatures + // allgather tempering arg "temp" across root procs + // bcast from each root to other procs in world + + set_temp = new double[nworlds]; + if (me == 0) MPI_Allgather(&temp,1,MPI_DOUBLE,set_temp,1,MPI_DOUBLE,roots); + MPI_Bcast(set_temp,nworlds,MPI_DOUBLE,0,world); + + // create world2temp only on root procs from my_set_temp + // create temp2world on root procs from world2temp, + // then bcast to all procs within world + + world2temp = new int[nworlds]; + temp2world = new int[nworlds]; + if (me == 0) { + MPI_Allgather(&my_set_temp,1,MPI_INT,world2temp,1,MPI_INT,roots); + for (int i = 0; i < nworlds; i++) temp2world[world2temp[i]] = i; + } + MPI_Bcast(temp2world,nworlds,MPI_INT,0,world); + + // if restarting tempering, reset temp target of Fix to current my_set_temp + + if (narg == 7) { + double new_temp = set_temp[my_set_temp]; + modify->fix[whichfix]->reset_target(new_temp); + } + + // setup tempering runs + + int i,which,partner,swap,partner_set_temp,partner_world; + double pe,pe_partner,boltz_factor,new_temp; + MPI_Status status; + + if (me_universe == 0 && universe->uscreen) + fprintf(universe->uscreen,"Setting up tempering ...\n"); + + update->integrate->setup(); + + if (me_universe == 0) { + if (universe->uscreen) { + fprintf(universe->uscreen,"Step"); + for (int i = 0; i < nworlds; i++) + fprintf(universe->uscreen," T%d",i); + fprintf(universe->uscreen,"\n"); + } + if (universe->ulogfile) { + fprintf(universe->ulogfile,"Step"); + for (int i = 0; i < nworlds; i++) + fprintf(universe->ulogfile," T%d",i); + fprintf(universe->ulogfile,"\n"); + } + print_status(); + } + + timer->barrier_start(TIME_LOOP); + + for (int iswap = 0; iswap < nswaps; iswap++) { + + // run for nevery timesteps + + update->integrate->run(nevery); + + // compute PE + // notify compute it will be called at next swap + + pe = pe_compute->compute_scalar(); + pe_compute->addstep(update->ntimestep + nevery); + + // which = which of 2 kinds of swaps to do (0,1) + + if (!ranswap) which = iswap % 2; + else if (ranswap->uniform() < 0.5) which = 0; + else which = 1; + + // partner_set_temp = which set temp I am partnering with for this swap + + if (which == 0) { + if (my_set_temp % 2 == 0) partner_set_temp = my_set_temp + 1; + else partner_set_temp = my_set_temp - 1; + } else { + if (my_set_temp % 2 == 1) partner_set_temp = my_set_temp + 1; + else partner_set_temp = my_set_temp - 1; + } + + // partner = proc ID to swap with + // if partner = -1, then I am not a proc that swaps + + partner = -1; + if (me == 0 && partner_set_temp >= 0 && partner_set_temp < nworlds) { + partner_world = temp2world[partner_set_temp]; + partner = world2root[partner_world]; + } + + // swap with a partner, only root procs in each world participate + // hi proc sends PE to low proc + // lo proc make Boltzmann decision on whether to swap + // lo proc communicates decision back to hi proc + + swap = 0; + if (partner != -1) { + if (me_universe > partner) + MPI_Send(&pe,1,MPI_DOUBLE,partner,0,universe->uworld); + else + MPI_Recv(&pe_partner,1,MPI_DOUBLE,partner,0,universe->uworld,&status); + + if (me_universe < partner) { + boltz_factor = (pe - pe_partner) * + (1.0/(boltz*set_temp[my_set_temp]) - + 1.0/(boltz*set_temp[partner_set_temp])); + if (boltz_factor >= 0.0) swap = 1; + else if (ranboltz->uniform() < exp(boltz_factor)) swap = 1; + } + + if (me_universe < partner) + MPI_Send(&swap,1,MPI_INT,partner,0,universe->uworld); + else + MPI_Recv(&swap,1,MPI_INT,partner,0,universe->uworld,&status); + +#ifdef TEMPER_DEBUG + if (me_universe < partner) + printf("SWAP %d & %d: yes = %d,Ts = %d %d, PEs = %g %g, Bz = %g %g\n", + me_universe,partner,swap,my_set_temp,partner_set_temp, + pe,pe_partner,boltz_factor,exp(boltz_factor)); +#endif + + } + + // bcast swap result to other procs in my world + + MPI_Bcast(&swap,1,MPI_INT,0,world); + + // rescale kinetic energy via velocities if move is accepted + + if (swap) scale_velocities(partner_set_temp,my_set_temp); + + // if my world swapped, all procs in world reset temp target of Fix + + if (swap) { + new_temp = set_temp[partner_set_temp]; + modify->fix[whichfix]->reset_target(new_temp); + } + + // update my_set_temp and temp2world on every proc + // root procs update their value if swap took place + // allgather across root procs + // bcast within my world + + if (swap) my_set_temp = partner_set_temp; + if (me == 0) { + MPI_Allgather(&my_set_temp,1,MPI_INT,world2temp,1,MPI_INT,roots); + for (i = 0; i < nworlds; i++) temp2world[world2temp[i]] = i; + } + MPI_Bcast(temp2world,nworlds,MPI_INT,0,world); + + // print out current swap status + + if (me_universe == 0) print_status(); + } + + timer->barrier_stop(TIME_LOOP); + + update->integrate->cleanup(); + + Finish finish(lmp); + finish.end(1); + + update->whichflag = 0; + update->firststep = update->laststep = 0; + update->beginstep = update->endstep = 0; +} + +/* ---------------------------------------------------------------------- + scale kinetic energy via velocities a la Sugita +------------------------------------------------------------------------- */ + +void Temper::scale_velocities(int t_partner, int t_me) +{ + double sfactor = sqrt(set_temp[t_partner]/set_temp[t_me]); + + double **v = atom->v; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + v[i][0] = v[i][0]*sfactor; + v[i][1] = v[i][1]*sfactor; + v[i][2] = v[i][2]*sfactor; + } +} + +/* ---------------------------------------------------------------------- + proc 0 prints current tempering status +------------------------------------------------------------------------- */ + +void Temper::print_status() +{ + if (universe->uscreen) { + fprintf(universe->uscreen,"%d ",update->ntimestep); + for (int i = 0; i < nworlds; i++) + fprintf(universe->uscreen,"%d ",world2temp[i]); + fprintf(universe->uscreen,"\n"); + } + if (universe->ulogfile) { + fprintf(universe->ulogfile,"%d ",update->ntimestep); + for (int i = 0; i < nworlds; i++) + fprintf(universe->ulogfile,"%d ",world2temp[i]); + fprintf(universe->ulogfile,"\n"); + fflush(universe->ulogfile); + } +}