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

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

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];
+    }
+}

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

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);
+    }
+  }
+}

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

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);
+  }
+}