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

src/MISC/fix_deposit.cpp

@@ -0,0 +1,508 @@
+/* ----------------------------------------------------------------------
+   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 "math.h"
+#include "stdlib.h"
+#include "string.h"
+#include "fix_deposit.h"
+#include "atom.h"
+#include "atom_vec.h"
+#include "force.h"
+#include "update.h"
+#include "modify.h"
+#include "fix.h"
+#include "comm.h"
+#include "domain.h"
+#include "lattice.h"
+#include "region.h"
+#include "random_park.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+/* ---------------------------------------------------------------------- */
+
+FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg)
+{
+  if (narg < 7) error->all(FLERR,"Illegal fix deposit command");
+
+  restart_global = 1;
+  time_depend = 1;
+
+  // required args
+
+  ninsert = force->inumeric(FLERR,arg[3]);
+  ntype = force->inumeric(FLERR,arg[4]);
+  nfreq = force->inumeric(FLERR,arg[5]);
+  seed = force->inumeric(FLERR,arg[6]);
+
+  if (seed <= 0) error->all(FLERR,"Illegal fix deposit command");
+
+  // set defaults
+
+  iregion = -1;
+  idregion = NULL;
+  idnext = 0;
+  globalflag = localflag = 0;
+  lo = hi = deltasq = 0.0;
+  nearsq = 0.0;
+  maxattempt = 10;
+  rateflag = 0;
+  vxlo = vxhi = vylo = vyhi = vzlo = vzhi = 0.0;
+  scaleflag = 1;
+  targetflag = 0;
+
+  // read options from end of input line
+
+  options(narg-7,&arg[7]);
+
+  // error checks on region and its extent being inside simulation box
+
+  if (iregion == -1) error->all(FLERR,"Must specify a region in fix deposit");
+  if (domain->regions[iregion]->bboxflag == 0)
+    error->all(FLERR,"Fix deposit region does not support a bounding box");
+  if (domain->regions[iregion]->dynamic_check())
+    error->all(FLERR,"Fix deposit region cannot be dynamic");
+
+  xlo = domain->regions[iregion]->extent_xlo;
+  xhi = domain->regions[iregion]->extent_xhi;
+  ylo = domain->regions[iregion]->extent_ylo;
+  yhi = domain->regions[iregion]->extent_yhi;
+  zlo = domain->regions[iregion]->extent_zlo;
+  zhi = domain->regions[iregion]->extent_zhi;
+
+  if (domain->triclinic == 0) {
+    if (xlo < domain->boxlo[0] || xhi > domain->boxhi[0] ||
+        ylo < domain->boxlo[1] || yhi > domain->boxhi[1] ||
+        zlo < domain->boxlo[2] || zhi > domain->boxhi[2])
+      error->all(FLERR,"Deposition region extends outside simulation box");
+  } else {
+    if (xlo < domain->boxlo_bound[0] || xhi > domain->boxhi_bound[0] ||
+        ylo < domain->boxlo_bound[1] || yhi > domain->boxhi_bound[1] ||
+        zlo < domain->boxlo_bound[2] || zhi > domain->boxhi_bound[2])
+      error->all(FLERR,"Deposition region extends outside simulation box");
+  }
+
+  // setup scaling
+
+  double xscale,yscale,zscale;
+  if (scaleflag) {
+    xscale = domain->lattice->xlattice;
+    yscale = domain->lattice->ylattice;
+    zscale = domain->lattice->zlattice;
+  }
+  else xscale = yscale = zscale = 1.0;
+
+  // apply scaling to all input parameters with dist/vel units
+
+  if (domain->dimension == 2) {
+    lo *= yscale;
+    hi *= yscale;
+    rate *= yscale;
+  } else {
+    lo *= zscale;
+    hi *= zscale;
+    rate *= zscale;
+  }
+  deltasq *= xscale*xscale;
+  nearsq *= xscale*xscale;
+  vxlo *= xscale;
+  vxhi *= xscale;
+  vylo *= yscale;
+  vyhi *= yscale;
+  vzlo *= zscale;
+  vzhi *= zscale;
+  tx *= xscale;
+  ty *= yscale;
+  tz *= zscale;
+
+  // maxtag_all = current max tag for all atoms
+
+  if (idnext) {
+    int *tag = atom->tag;
+    int nlocal = atom->nlocal;
+
+    int maxtag = 0;
+    for (int i = 0; i < nlocal; i++) maxtag = MAX(maxtag,tag[i]);
+    MPI_Allreduce(&maxtag,&maxtag_all,1,MPI_INT,MPI_MAX,world);
+  }
+
+  // random number generator, same for all procs
+
+  random = new RanPark(lmp,seed);
+
+  // set up reneighboring
+
+  force_reneighbor = 1;
+  next_reneighbor = update->ntimestep + 1;
+  nfirst = next_reneighbor;
+  ninserted = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixDeposit::~FixDeposit()
+{
+  delete random;
+  delete [] idregion;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixDeposit::setmask()
+{
+  int mask = 0;
+  mask |= PRE_EXCHANGE;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixDeposit::init()
+{
+  // set index and check validity of region
+
+  iregion = domain->find_region(idregion);
+  if (iregion == -1)
+    error->all(FLERR,"Region ID for fix deposit does not exist");
+}
+
+/* ----------------------------------------------------------------------
+   perform particle insertion
+------------------------------------------------------------------------- */
+
+void FixDeposit::pre_exchange()
+{
+  int i,j;
+  int flag,flagall;
+  double coord[3],lamda[3],delx,dely,delz,rsq;
+  double *newcoord;
+
+  // just return if should not be called on this timestep
+
+  if (next_reneighbor != update->ntimestep) return;
+
+  // compute current offset = bottom of insertion volume
+
+  double offset = 0.0;
+  if (rateflag) offset = (update->ntimestep - nfirst) * update->dt * rate;
+
+  double *sublo,*subhi;
+  if (domain->triclinic == 0) {
+    sublo = domain->sublo;
+    subhi = domain->subhi;
+  } else {
+    sublo = domain->sublo_lamda;
+    subhi = domain->subhi_lamda;
+  }
+
+  // attempt an insertion until successful
+
+  int nfix = modify->nfix;
+  Fix **fix = modify->fix;
+
+  int success = 0;
+  int attempt = 0;
+  while (attempt < maxattempt) {
+    attempt++;
+
+    // choose random position for new atom within region
+
+    coord[0] = xlo + random->uniform() * (xhi-xlo);
+    coord[1] = ylo + random->uniform() * (yhi-ylo);
+    coord[2] = zlo + random->uniform() * (zhi-zlo);
+    while (domain->regions[iregion]->match(coord[0],coord[1],coord[2]) == 0) {
+      coord[0] = xlo + random->uniform() * (xhi-xlo);
+      coord[1] = ylo + random->uniform() * (yhi-ylo);
+      coord[2] = zlo + random->uniform() * (zhi-zlo);
+    }
+
+    // adjust vertical coord by offset
+
+    if (domain->dimension == 2) coord[1] += offset;
+    else coord[2] += offset;
+
+    // if global, reset vertical coord to be lo-hi above highest atom
+    // if local, reset vertical coord to be lo-hi above highest "nearby" atom
+    // local computation computes lateral distance between 2 particles w/ PBC
+
+    if (globalflag || localflag) {
+      int dim;
+      double max,maxall,delx,dely,delz,rsq;
+
+      if (domain->dimension == 2) {
+        dim = 1;
+        max = domain->boxlo[1];
+      } else {
+        dim = 2;
+        max = domain->boxlo[2];
+      }
+
+      double **x = atom->x;
+      int nlocal = atom->nlocal;
+      for (i = 0; i < nlocal; i++) {
+        if (localflag) {
+          delx = coord[0] - x[i][0];
+          dely = coord[1] - x[i][1];
+          delz = 0.0;
+          domain->minimum_image(delx,dely,delz);
+          if (domain->dimension == 2) rsq = delx*delx;
+          else rsq = delx*delx + dely*dely;
+          if (rsq > deltasq) continue;
+        }
+        if (x[i][dim] > max) max = x[i][dim];
+      }
+
+      MPI_Allreduce(&max,&maxall,1,MPI_DOUBLE,MPI_MAX,world);
+      if (domain->dimension == 2)
+        coord[1] = maxall + lo + random->uniform()*(hi-lo);
+      else
+        coord[2] = maxall + lo + random->uniform()*(hi-lo);
+    }
+
+    // now have final coord
+    // if distance to any atom is less than near, try again
+
+    double **x = atom->x;
+    int nlocal = atom->nlocal;
+
+    flag = 0;
+    for (i = 0; i < nlocal; i++) {
+      delx = coord[0] - x[i][0];
+      dely = coord[1] - x[i][1];
+      delz = coord[2] - x[i][2];
+      domain->minimum_image(delx,dely,delz);
+      rsq = delx*delx + dely*dely + delz*delz;
+      if (rsq < nearsq) flag = 1;
+    }
+    MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MAX,world);
+    if (flagall) continue;
+
+    // insertion will proceed
+    // choose random velocity for new atom
+
+    double vxtmp = vxlo + random->uniform() * (vxhi-vxlo);
+    double vytmp = vylo + random->uniform() * (vyhi-vylo);
+    double vztmp = vzlo + random->uniform() * (vzhi-vzlo);
+
+    // if we have a sputter target change velocity vector accordingly
+    if (targetflag) {
+      double vel = sqrt(vxtmp*vxtmp + vytmp*vytmp + vztmp*vztmp);
+      delx = tx - coord[0];
+      dely = ty - coord[1];
+      delz = tz - coord[2];
+      double rsq = delx*delx + dely*dely + delz*delz;
+      if (rsq > 0.0) {
+        double rinv = sqrt(1.0/rsq);
+        vxtmp = delx*rinv*vel;
+        vytmp = dely*rinv*vel;
+        vztmp = delz*rinv*vel;
+      }
+    }
+
+    // check if new atom is in my sub-box or above it if I'm highest proc
+    // if so, add to my list via create_atom()
+    // initialize info about the atoms
+    // set group mask to "all" plus fix group
+
+    if (domain->triclinic) {
+      domain->x2lamda(coord,lamda);
+      newcoord = lamda;
+    } else newcoord = coord;
+
+    flag = 0;
+    if (newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] &&
+        newcoord[1] >= sublo[1] && newcoord[1] < subhi[1] &&
+        newcoord[2] >= sublo[2] && newcoord[2] < subhi[2]) flag = 1;
+    else if (domain->dimension == 3 && newcoord[2] >= domain->boxhi[2] &&
+             comm->myloc[2] == comm->procgrid[2]-1 &&
+             newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] &&
+             newcoord[1] >= sublo[1] && newcoord[1] < subhi[1]) flag = 1;
+    else if (domain->dimension == 2 && newcoord[1] >= domain->boxhi[1] &&
+             comm->myloc[1] == comm->procgrid[1]-1 &&
+             newcoord[0] >= sublo[0] && newcoord[0] < subhi[0]) flag = 1;
+
+    if (flag) {
+      atom->avec->create_atom(ntype,coord);
+      int m = atom->nlocal - 1;
+      atom->type[m] = ntype;
+      atom->mask[m] = 1 | groupbit;
+      atom->v[m][0] = vxtmp;
+      atom->v[m][1] = vytmp;
+      atom->v[m][2] = vztmp;
+      for (j = 0; j < nfix; j++)
+        if (fix[j]->create_attribute) fix[j]->set_arrays(m);
+    }
+    MPI_Allreduce(&flag,&success,1,MPI_INT,MPI_MAX,world);
+    break;
+  }
+
+  // warn if not successful b/c too many attempts or no proc owned particle
+
+  if (!success && comm->me == 0)
+    error->warning(FLERR,"Particle deposition was unsuccessful",0);
+
+  // reset global natoms
+  // if idnext, set new atom ID to incremented maxtag_all
+  // else set new atom ID to value beyond all current atoms
+  // if global map exists, reset it now instead of waiting for comm
+  // since adding an atom messes up ghosts
+
+  if (success) {
+    atom->natoms += 1;
+    if (atom->tag_enable) {
+      if (idnext) {
+        maxtag_all++;
+        if (atom->nlocal && atom->tag[atom->nlocal-1] == 0) 
+          atom->tag[atom->nlocal-1] = maxtag_all;
+      } else atom->tag_extend();
+      if (atom->map_style) {
+        atom->nghost = 0;
+        atom->map_init();
+        atom->map_set();
+      }
+    }
+  }
+
+  // next timestep to insert
+  // next_reneighbor = 0 if done
+
+  if (success) ninserted++;
+  if (ninserted < ninsert) next_reneighbor += nfreq;
+  else next_reneighbor = 0;
+}
+
+/* ----------------------------------------------------------------------
+   parse optional parameters at end of input line
+------------------------------------------------------------------------- */
+
+void FixDeposit::options(int narg, char **arg)
+{
+  if (narg < 0) error->all(FLERR,"Illegal fix indent command");
+
+  int iarg = 0;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"region") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");
+      iregion = domain->find_region(arg[iarg+1]);
+      if (iregion == -1)
+        error->all(FLERR,"Region ID for fix deposit does not exist");
+      int n = strlen(arg[iarg+1]) + 1;
+      idregion = new char[n];
+      strcpy(idregion,arg[iarg+1]);
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"id") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");
+      if (strcmp(arg[iarg+1],"max") == 0) idnext = 0;
+      else if (strcmp(arg[iarg+1],"next") == 0) idnext = 1;
+      else error->all(FLERR,"Illegal fix deposit command");
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"global") == 0) {
+      if (iarg+3 > narg) error->all(FLERR,"Illegal fix deposit command");
+      globalflag = 1;
+      localflag = 0;
+      lo = force->numeric(FLERR,arg[iarg+1]);
+      hi = force->numeric(FLERR,arg[iarg+2]);
+      iarg += 3;
+    } else if (strcmp(arg[iarg],"local") == 0) {
+      if (iarg+4 > narg) error->all(FLERR,"Illegal fix deposit command");
+      localflag = 1;
+      globalflag = 0;
+      lo = force->numeric(FLERR,arg[iarg+1]);
+      hi = force->numeric(FLERR,arg[iarg+2]);
+      deltasq = force->numeric(FLERR,arg[iarg+3])*force->numeric(FLERR,arg[iarg+3]);
+      iarg += 4;
+    } else if (strcmp(arg[iarg],"near") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");
+      nearsq = force->numeric(FLERR,arg[iarg+1])*force->numeric(FLERR,arg[iarg+1]);
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"attempt") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");
+      maxattempt = force->inumeric(FLERR,arg[iarg+1]);
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"rate") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");
+      rateflag = 1;
+      rate = force->numeric(FLERR,arg[iarg+1]);
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"vx") == 0) {
+      if (iarg+3 > narg) error->all(FLERR,"Illegal fix deposit command");
+      vxlo = force->numeric(FLERR,arg[iarg+1]);
+      vxhi = force->numeric(FLERR,arg[iarg+2]);
+      iarg += 3;
+    } else if (strcmp(arg[iarg],"vy") == 0) {
+      if (iarg+3 > narg) error->all(FLERR,"Illegal fix deposit command");
+      vylo = force->numeric(FLERR,arg[iarg+1]);
+      vyhi = force->numeric(FLERR,arg[iarg+2]);
+      iarg += 3;
+    } else if (strcmp(arg[iarg],"vz") == 0) {
+      if (iarg+3 > narg) error->all(FLERR,"Illegal fix deposit command");
+      vzlo = force->numeric(FLERR,arg[iarg+1]);
+      vzhi = force->numeric(FLERR,arg[iarg+2]);
+      iarg += 3;
+    } else if (strcmp(arg[iarg],"units") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");
+      if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
+      else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
+      else error->all(FLERR,"Illegal fix deposit command");
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"target") == 0) {
+      if (iarg+4 > narg) error->all(FLERR,"Illegal fix deposit command");
+      tx = force->numeric(FLERR,arg[iarg+1]);
+      ty = force->numeric(FLERR,arg[iarg+2]);
+      tz = force->numeric(FLERR,arg[iarg+3]);
+      targetflag = 1;
+      iarg += 4;
+    } else error->all(FLERR,"Illegal fix deposit command");
+  }
+}
+
+/* ----------------------------------------------------------------------
+   pack entire state of Fix into one write
+------------------------------------------------------------------------- */
+
+void FixDeposit::write_restart(FILE *fp)
+{
+  int n = 0;
+  double list[4];
+  list[n++] = random->state();
+  list[n++] = ninserted;
+  list[n++] = nfirst;
+  list[n++] = next_reneighbor;
+
+  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 FixDeposit::restart(char *buf)
+{
+  int n = 0;
+  double *list = (double *) buf;
+
+  seed = static_cast<int> (list[n++]);
+  ninserted = static_cast<int> (list[n++]);
+  nfirst = static_cast<int> (list[n++]);
+  next_reneighbor = static_cast<int> (list[n++]);
+
+  random->reset(seed);
+}

src/MISC/fix_deposit.h

@@ -0,0 +1,98 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   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(deposit,FixDeposit)
+
+#else
+
+#ifndef LMP_FIX_DEPOSIT_H
+#define LMP_FIX_DEPOSIT_H
+
+#include "stdio.h"
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixDeposit : public Fix {
+ public:
+  FixDeposit(class LAMMPS *, int, char **);
+  ~FixDeposit();
+  int setmask();
+  void init();
+  void pre_exchange();
+  void write_restart(FILE *);
+  void restart(char *);
+
+ private:
+  int ninsert,ntype,nfreq,seed;
+  int iregion,globalflag,localflag,maxattempt,rateflag,scaleflag,targetflag;
+  char *idregion;
+  double lo,hi,deltasq,nearsq,rate;
+  double vxlo,vxhi,vylo,vyhi,vzlo,vzhi;
+  double xlo,xhi,ylo,yhi,zlo,zhi;
+  double tx,ty,tz;
+  int nfirst,ninserted;
+  int idnext,maxtag_all;
+  class RanPark *random;
+
+  void options(int, char **);
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Must specify a region in fix deposit
+
+The region keyword must be specified with this fix.
+
+E: Fix deposit region does not support a bounding box
+
+Not all regions represent bounded volumes.  You cannot use
+such a region with the fix deposit command.
+
+E: Fix deposit region cannot be dynamic
+
+Only static regions can be used with fix deposit.
+
+E: Deposition region extends outside simulation box
+
+Self-explanatory.
+
+E: Region ID for fix deposit does not exist
+
+Self-explanatory.
+
+W: Particle deposition was unsuccessful
+
+The fix deposit command was not able to insert as many atoms as
+needed.  The requested volume fraction may be too high, or other atoms
+may be in the insertion region.
+
+U: Use of fix deposit with undefined lattice
+
+Must use lattice command with compute fix deposit command if units
+option is set to lattice.
+
+*/

src/MISC/fix_efield.cpp

@@ -0,0 +1,406 @@
+/* ----------------------------------------------------------------------
+   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: Christina Payne (Vanderbilt U)
+                        Stan Moore (Sandia) for dipole terms
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "string.h"
+#include "stdlib.h"
+#include "fix_efield.h"
+#include "atom.h"
+#include "update.h"
+#include "domain.h"
+#include "comm.h"
+#include "modify.h"
+#include "force.h"
+#include "respa.h"
+#include "input.h"
+#include "variable.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+enum{NONE,CONSTANT,EQUAL,ATOM};
+
+/* ---------------------------------------------------------------------- */
+
+FixEfield::FixEfield(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg)
+{
+  if (narg != 6) error->all(FLERR,"Illegal fix efield command");
+
+  vector_flag = 1;
+  scalar_flag = 1;
+  size_vector = 3;
+  global_freq = 1;
+  extvector = 1;
+  extscalar = 1;
+
+  qe2f = force->qe2f;
+  xstr = ystr = zstr = NULL;
+
+  if (strstr(arg[3],"v_") == arg[3]) {
+    int n = strlen(&arg[3][2]) + 1;
+    xstr = new char[n];
+    strcpy(xstr,&arg[3][2]);
+  } else {
+    ex = qe2f * force->numeric(FLERR,arg[3]);
+    xstyle = CONSTANT;
+  }
+
+  if (strstr(arg[4],"v_") == arg[4]) {
+    int n = strlen(&arg[4][2]) + 1;
+    ystr = new char[n];
+    strcpy(ystr,&arg[4][2]);
+  } else {
+    ey = qe2f * force->numeric(FLERR,arg[4]);
+    ystyle = CONSTANT;
+  }
+
+  if (strstr(arg[5],"v_") == arg[5]) {
+    int n = strlen(&arg[5][2]) + 1;
+    zstr = new char[n];
+    strcpy(zstr,&arg[5][2]);
+  } else {
+    ez = qe2f * force->numeric(FLERR,arg[5]);
+    zstyle = CONSTANT;
+  }
+
+  // optional args
+
+  estr = NULL;
+
+  int iarg = 6;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"energy") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix efield command");
+      if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
+        int n = strlen(&arg[iarg+1][2]) + 1;
+        estr = new char[n];
+        strcpy(estr,&arg[iarg+1][2]);
+      } else error->all(FLERR,"Illegal fix efield command");
+      iarg += 2;
+    } else error->all(FLERR,"Illegal fix efield command");
+  }
+
+  force_flag = 0;
+  fsum[0] = fsum[1] = fsum[2] = fsum[3] = 0.0;
+
+  maxatom = 0;
+  efield = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixEfield::~FixEfield()
+{
+  delete [] xstr;
+  delete [] ystr;
+  delete [] zstr;
+  delete [] estr;
+  memory->destroy(efield);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixEfield::setmask()
+{
+  int mask = 0;
+  mask |= THERMO_ENERGY;
+  mask |= POST_FORCE;
+  mask |= POST_FORCE_RESPA;
+  mask |= MIN_POST_FORCE;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixEfield::init()
+{
+  qflag = muflag = 0;
+  if (atom->q_flag) qflag = 1;
+  if (atom->mu_flag && atom->torque_flag) muflag = 1;
+  if (!qflag && !muflag) 
+    error->all(FLERR,"Fix efield requires atom attribute q or mu");
+
+  // check variables
+
+  if (xstr) {
+    xvar = input->variable->find(xstr);
+    if (xvar < 0) 
+      error->all(FLERR,"Variable name for fix efield does not exist");
+    if (input->variable->equalstyle(xvar)) xstyle = EQUAL;
+    else if (input->variable->atomstyle(xvar)) xstyle = ATOM;
+    else error->all(FLERR,"Variable for fix efield is invalid style");
+  }
+  if (ystr) {
+    yvar = input->variable->find(ystr);
+    if (yvar < 0) 
+      error->all(FLERR,"Variable name for fix efield does not exist");
+    if (input->variable->equalstyle(yvar)) ystyle = EQUAL;
+    else if (input->variable->atomstyle(yvar)) ystyle = ATOM;
+    else error->all(FLERR,"Variable for fix efield is invalid style");
+  }
+  if (zstr) {
+    zvar = input->variable->find(zstr);
+    if (zvar < 0) 
+      error->all(FLERR,"Variable name for fix efield does not exist");
+    if (input->variable->equalstyle(zvar)) zstyle = EQUAL;
+    else if (input->variable->atomstyle(zvar)) zstyle = ATOM;
+    else error->all(FLERR,"Variable for fix efield is invalid style");
+  }
+  if (estr) {
+    evar = input->variable->find(estr);
+    if (evar < 0) 
+      error->all(FLERR,"Variable name for fix efield does not exist");
+    if (input->variable->atomstyle(evar)) estyle = ATOM;
+    else error->all(FLERR,"Variable for fix efield is invalid style");
+  } else estyle = NONE;
+
+  if (xstyle == ATOM || ystyle == ATOM || zstyle == ATOM)
+    varflag = ATOM;
+  else if (xstyle == EQUAL || ystyle == EQUAL || zstyle == EQUAL)
+    varflag = EQUAL;
+  else varflag = CONSTANT;
+
+  if (muflag && varflag == ATOM)
+    error->all(FLERR,"Fix efield with dipoles cannot use atom-style variables");
+
+  if (muflag && update->whichflag == 2 && comm->me == 0)
+    error->warning(FLERR,
+                   "The minimizer does not re-orient dipoles "
+                   "when using fix efield");
+
+  if (varflag == CONSTANT && estyle != NONE)
+    error->all(FLERR,"Cannot use variable energy with "
+               "constant efield in fix efield");
+  if ((varflag == EQUAL || varflag == ATOM) &&
+      update->whichflag == 2 && estyle == NONE)
+    error->all(FLERR,"Must use variable energy with fix efield");
+
+  if (strstr(update->integrate_style,"respa"))
+    nlevels_respa = ((Respa *) update->integrate)->nlevels;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixEfield::setup(int vflag)
+{
+  if (strstr(update->integrate_style,"verlet"))
+    post_force(vflag);
+  else {
+    ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
+    post_force_respa(vflag,nlevels_respa-1,0);
+    ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixEfield::min_setup(int vflag)
+{
+  post_force(vflag);
+}
+
+/* ----------------------------------------------------------------------
+   apply F = qE
+------------------------------------------------------------------------- */
+
+void FixEfield::post_force(int vflag)
+{
+  double **f = atom->f;
+  double *q = atom->q;
+  int *mask = atom->mask;
+  tagint *image = atom->image;
+  int nlocal = atom->nlocal;
+
+  // reallocate efield array if necessary
+
+  if (varflag == ATOM && nlocal > maxatom) {
+    maxatom = atom->nmax;
+    memory->destroy(efield);
+    memory->create(efield,maxatom,4,"efield:efield");
+  }
+
+  // fsum[0] = "potential energy" for added force
+  // fsum[123] = extra force added to atoms
+
+  fsum[0] = fsum[1] = fsum[2] = fsum[3] = 0.0;
+  force_flag = 0;
+
+  double **x = atom->x;
+  double fx,fy,fz;
+
+  // constant efield
+
+  if (varflag == CONSTANT) {
+    double unwrap[3];
+
+    // charge interactions
+    // force = qE, potential energy = F dot x in unwrapped coords
+
+    if (qflag) {
+      for (int i = 0; i < nlocal; i++)
+        if (mask[i] & groupbit) {
+          fx = q[i]*ex;
+          fy = q[i]*ey;
+          fz = q[i]*ez;
+          f[i][0] += fx;
+          f[i][1] += fy;
+          f[i][2] += fz;
+          
+          domain->unmap(x[i],image[i],unwrap);
+          fsum[0] -= fx*unwrap[0]+fy*unwrap[1]+fz*unwrap[2];
+          fsum[1] += fx;
+          fsum[2] += fy;
+          fsum[3] += fz;
+        }
+    }
+
+    // dipole interactions
+    // no force, torque = mu cross E, potential energy = -mu dot E
+
+    if (muflag) {
+      double **mu = atom->mu;
+      double **t = atom->torque;
+      double tx,ty,tz;
+      for (int i = 0; i < nlocal; i++)
+        if (mask[i] & groupbit) {
+          tx = ez*mu[i][1] - ey*mu[i][2];
+          ty = ex*mu[i][2] - ez*mu[i][0];
+          tz = ey*mu[i][0] - ex*mu[i][1];
+          t[i][0] += tx;
+          t[i][1] += ty;
+          t[i][2] += tz;
+          fsum[0] -= mu[i][0]*ex + mu[i][1]*ey + mu[i][2]*ez;
+        }
+    }
+
+  // variable efield, wrap with clear/add
+  // potential energy = evar if defined, else 0.0
+
+  } else {
+
+    modify->clearstep_compute();
+
+    if (xstyle == EQUAL) ex = qe2f * input->variable->compute_equal(xvar);
+    else if (xstyle == ATOM && efield)
+      input->variable->compute_atom(xvar,igroup,&efield[0][0],3,0);
+    if (ystyle == EQUAL) ey = qe2f * input->variable->compute_equal(yvar);
+    else if (ystyle == ATOM && efield)
+      input->variable->compute_atom(yvar,igroup,&efield[0][1],3,0);
+    if (zstyle == EQUAL) ez = qe2f * input->variable->compute_equal(zvar);
+    else if (zstyle == ATOM && efield)
+      input->variable->compute_atom(zvar,igroup,&efield[0][2],3,0);
+    if (estyle == ATOM && efield)
+      input->variable->compute_atom(evar,igroup,&efield[0][3],4,0);
+
+    modify->addstep_compute(update->ntimestep + 1);
+
+    // charge interactions
+    // force = qE
+
+    if (qflag) {
+      for (int i = 0; i < nlocal; i++)
+        if (mask[i] & groupbit) {
+          if (xstyle == ATOM) fx = qe2f * q[i]*efield[i][0];
+          else fx = q[i]*ex;
+          f[i][0] += fx;
+          fsum[1] += fx;
+          if (ystyle == ATOM) fy = qe2f * q[i]*efield[i][1];
+          else fy = q[i]*ey;
+          f[i][1] += fy;
+          fsum[2] += fy;
+          if (zstyle == ATOM) fz = qe2f * q[i]*efield[i][2];
+          else fz = q[i]*ez;
+          f[i][2] += fz;
+          fsum[3] += fz;
+          if (estyle == ATOM) fsum[0] += efield[0][3];
+        }
+    }
+
+    // dipole interactions
+    // no force, torque = mu cross E
+
+    if (muflag) {
+      double **mu = atom->mu;
+      double **t = atom->torque;
+      double tx,ty,tz;
+      for (int i = 0; i < nlocal; i++)
+        if (mask[i] & groupbit) {
+          tx = ez*mu[i][1] - ey*mu[i][2];
+          ty = ex*mu[i][2] - ez*mu[i][0];
+          tz = ey*mu[i][0] - ex*mu[i][1];
+          t[i][0] += tx;
+          t[i][1] += ty;
+          t[i][2] += tz;
+        }
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixEfield::post_force_respa(int vflag, int ilevel, int iloop)
+{
+  if (ilevel == nlevels_respa-1) post_force(vflag);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixEfield::min_post_force(int vflag)
+{
+  post_force(vflag);
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based array
+------------------------------------------------------------------------- */
+
+double FixEfield::memory_usage()
+{
+  double bytes = 0.0;
+  if (varflag == ATOM) bytes = atom->nmax*4 * sizeof(double);
+  return bytes;
+} 
+
+/* ----------------------------------------------------------------------
+   return energy added by fix
+------------------------------------------------------------------------- */
+
+double FixEfield::compute_scalar(void)
+{
+  if (force_flag == 0) {
+    MPI_Allreduce(fsum,fsum_all,4,MPI_DOUBLE,MPI_SUM,world);
+    force_flag = 1;
+  }
+  return fsum_all[0];
+}
+
+/* ----------------------------------------------------------------------
+   return total extra force due to fix
+------------------------------------------------------------------------- */
+
+double FixEfield::compute_vector(int n)
+{
+  if (force_flag == 0) {
+    MPI_Allreduce(fsum,fsum_all,4,MPI_DOUBLE,MPI_SUM,world);
+    force_flag = 1;
+  }
+  return fsum_all[n+1];
+}

src/MISC/fix_efield.h

@@ -0,0 +1,101 @@
+/* ----------------------------------------------------------------------
+   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(efield,FixEfield)
+
+#else
+
+#ifndef LMP_FIX_EFIELD_H
+#define LMP_FIX_EFIELD_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixEfield : public Fix {
+ public:
+  FixEfield(class LAMMPS *, int, char **);
+  ~FixEfield();
+  int setmask();
+  void init();
+  void setup(int);
+  void min_setup(int);
+  void post_force(int);
+  void post_force_respa(int, int, int);
+  void min_post_force(int);
+  double memory_usage();
+  double compute_scalar();
+  double compute_vector(int);
+
+ private:
+  double ex,ey,ez;
+  int varflag;
+  char *xstr,*ystr,*zstr,*estr;
+  int xvar,yvar,zvar,evar,xstyle,ystyle,zstyle,estyle;
+  int nlevels_respa;
+  double qe2f;
+  double fdotx;
+  int qflag,muflag;
+
+  int maxatom;
+  double **efield;
+
+  int force_flag;
+  double fsum[4],fsum_all[4];
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Fix efield requires atom attribute q or mu
+
+Self-explanatory.
+
+E: Variable name for fix efield does not exist
+
+Self-explanatory.
+
+E: Variable for fix efield is invalid style
+
+Only equal-style or atom-style variables can be used.
+
+E: Fix efield with dipoles cannot use atom-style variables
+
+This feature is not yet supported.
+
+W: The minimizer does not re-orient dipoles when using fix efield
+
+Self-explanatory.
+
+E: Cannot use variable energy with constant efield in fix efield
+
+Self-explanatory.
+
+E: Must use variable energy with fix efield
+
+One or more variables are defined for fix efield, which require
+variable energy when using the minimizer.
+
+*/

src/MISC/fix_evaporate.cpp

@@ -0,0 +1,386 @@
+/* ----------------------------------------------------------------------
+   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 "math.h"
+#include "stdlib.h"
+#include "string.h"
+#include "fix_evaporate.h"
+#include "atom.h"
+#include "atom_vec.h"
+#include "update.h"
+#include "domain.h"
+#include "region.h"
+#include "comm.h"
+#include "force.h"
+#include "group.h"
+#include "random_park.h"
+#include "random_mars.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+/* ---------------------------------------------------------------------- */
+
+FixEvaporate::FixEvaporate(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg)
+{
+  if (narg < 7) error->all(FLERR,"Illegal fix evaporate command");
+
+  scalar_flag = 1;
+  global_freq = 1;
+  extscalar = 0;
+
+  nevery = force->inumeric(FLERR,arg[3]);
+  nflux = force->inumeric(FLERR,arg[4]);
+  iregion = domain->find_region(arg[5]);
+  int n = strlen(arg[5]) + 1;
+  idregion = new char[n];
+  strcpy(idregion,arg[5]);
+  int seed = force->inumeric(FLERR,arg[6]);
+
+  if (nevery <= 0 || nflux <= 0)
+    error->all(FLERR,"Illegal fix evaporate command");
+  if (iregion == -1)
+    error->all(FLERR,"Region ID for fix evaporate does not exist");
+  if (seed <= 0) error->all(FLERR,"Illegal fix evaporate command");
+
+  // random number generator, same for all procs
+
+  random = new RanPark(lmp,seed);
+
+  // optional args
+
+  molflag = 0;
+
+  int iarg = 7;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"molecule") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix evaporate command");
+      if (strcmp(arg[iarg+1],"no") == 0) molflag = 0;
+      else if (strcmp(arg[iarg+1],"yes") == 0) molflag = 1;
+      else error->all(FLERR,"Illegal fix evaporate command");
+      iarg += 2;
+    } else error->all(FLERR,"Illegal fix evaporate command");
+  }
+
+  // set up reneighboring
+
+  force_reneighbor = 1;
+  next_reneighbor = (update->ntimestep/nevery)*nevery + nevery;
+  ndeleted = 0;
+
+  nmax = 0;
+  list = NULL;
+  mark = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixEvaporate::~FixEvaporate()
+{
+  delete [] idregion;
+  delete random;
+  memory->destroy(list);
+  memory->destroy(mark);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixEvaporate::setmask()
+{
+  int mask = 0;
+  mask |= PRE_EXCHANGE;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixEvaporate::init()
+{
+  // set index and check validity of region
+
+  iregion = domain->find_region(idregion);
+  if (iregion == -1)
+    error->all(FLERR,"Region ID for fix evaporate does not exist");
+
+  // check that no deletable atoms are in atom->firstgroup
+  // deleting such an atom would not leave firstgroup atoms first
+
+  if (atom->firstgroup >= 0) {
+    int *mask = atom->mask;
+    int nlocal = atom->nlocal;
+    int firstgroupbit = group->bitmask[atom->firstgroup];
+
+    int flag = 0;
+    for (int i = 0; i < nlocal; i++)
+      if ((mask[i] & groupbit) && (mask[i] && firstgroupbit)) flag = 1;
+
+    int flagall;
+    MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
+
+    if (flagall)
+      error->all(FLERR,"Cannot evaporate atoms in atom_modify first group");
+  }
+
+  // if molflag not set, warn if any deletable atom has a mol ID
+
+  if (molflag == 0 && atom->molecule_flag) {
+    int *molecule = atom->molecule;
+    int *mask = atom->mask;
+    int nlocal = atom->nlocal;
+    int flag = 0;
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit)
+        if (molecule[i]) flag = 1;
+    int flagall;
+    MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
+    if (flagall && comm->me == 0)
+      error->warning(FLERR,
+                     "Fix evaporate may delete atom with non-zero molecule ID");
+  }
+
+  if (molflag && atom->molecule_flag == 0)
+      error->all(FLERR,
+                 "Fix evaporate molecule requires atom attribute molecule");
+}
+
+/* ----------------------------------------------------------------------
+   perform particle deletion
+   done before exchange, borders, reneighbor
+   so that ghost atoms and neighbor lists will be correct
+------------------------------------------------------------------------- */
+
+void FixEvaporate::pre_exchange()
+{
+  int i,j,m,iwhichglobal,iwhichlocal;
+  int ndel,ndeltopo[4];
+
+  if (update->ntimestep != next_reneighbor) return;
+
+  // grow list and mark arrays if necessary
+
+  if (atom->nlocal > nmax) {
+    memory->destroy(list);
+    memory->destroy(mark);
+    nmax = atom->nmax;
+    memory->create(list,nmax,"evaporate:list");
+    memory->create(mark,nmax,"evaporate:mark");
+  }
+
+  // ncount = # of deletable atoms in region that I own
+  // nall = # on all procs
+  // nbefore = # on procs before me
+  // list[ncount] = list of local indices of atoms I can delete
+
+  double **x = atom->x;
+  int *mask = atom->mask;
+  int *tag = atom->tag;
+  int nlocal = atom->nlocal;
+
+  int ncount = 0;
+  for (i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit)
+      if (domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2]))
+        list[ncount++] = i;
+
+  int nall,nbefore;
+  MPI_Allreduce(&ncount,&nall,1,MPI_INT,MPI_SUM,world);
+  MPI_Scan(&ncount,&nbefore,1,MPI_INT,MPI_SUM,world);
+  nbefore -= ncount;
+
+  // ndel = total # of atom deletions, in or out of region
+  // ndeltopo[1,2,3,4] = ditto for bonds, angles, dihedrals, impropers
+  // mark[] = 1 if deleted
+
+  ndel = 0;
+  for (i = 0; i < nlocal; i++) mark[i] = 0;
+
+  // atomic deletions
+  // choose atoms randomly across all procs and mark them for deletion
+  // shrink eligible list as my atoms get marked
+  // keep ndel,ncount,nall,nbefore current after each atom deletion
+
+  if (molflag == 0) {
+    while (nall && ndel < nflux) {
+      iwhichglobal = static_cast<int> (nall*random->uniform());
+      if (iwhichglobal < nbefore) nbefore--;
+      else if (iwhichglobal < nbefore + ncount) {
+        iwhichlocal = iwhichglobal - nbefore;
+        mark[list[iwhichlocal]] = 1;
+        list[iwhichlocal] = list[ncount-1];
+        ncount--;
+      }
+      ndel++;
+      nall--;
+    }
+
+  // molecule deletions
+  // choose one atom in one molecule randomly across all procs
+  // bcast mol ID and delete all atoms in that molecule on any proc
+  // update deletion count by total # of atoms in molecule
+  // shrink list of eligible candidates as any of my atoms get marked
+  // keep ndel,ndeltopo,ncount,nall,nbefore current after each mol deletion
+
+  } else {
+    int me,proc,iatom,imolecule,ndelone,ndelall;
+    int *molecule = atom->molecule;
+
+    ndeltopo[0] = ndeltopo[1] = ndeltopo[2] = ndeltopo[3] = 0;
+
+    while (nall && ndel < nflux) {
+
+      // pick an iatom,imolecule on proc me to delete
+
+      iwhichglobal = static_cast<int> (nall*random->uniform());
+      if (iwhichglobal >= nbefore && iwhichglobal < nbefore + ncount) {
+        iwhichlocal = iwhichglobal - nbefore;
+        iatom = list[iwhichlocal];
+        imolecule = molecule[iatom];
+        me = comm->me;
+      } else me = -1;
+
+      // bcast mol ID to delete all atoms from
+      // if mol ID > 0, delete any atom in molecule and decrement counters
+      // if mol ID == 0, delete single iatom
+      // be careful to delete correct # of bond,angle,etc for newton on or off
+
+      MPI_Allreduce(&me,&proc,1,MPI_INT,MPI_MAX,world);
+      MPI_Bcast(&imolecule,1,MPI_INT,proc,world);
+      ndelone = 0;
+      for (i = 0; i < nlocal; i++) {
+        if (imolecule && molecule[i] == imolecule) {
+          mark[i] = 1;
+          ndelone++;
+
+          if (atom->avec->bonds_allow) {
+            if (force->newton_bond) ndeltopo[0] += atom->num_bond[i];
+            else {
+              for (j = 0; j < atom->num_bond[i]; j++) {
+                if (tag[i] < atom->bond_atom[i][j]) ndeltopo[0]++;
+              }
+            }
+          }
+          if (atom->avec->angles_allow) {
+            if (force->newton_bond) ndeltopo[1] += atom->num_angle[i];
+            else {
+              for (j = 0; j < atom->num_angle[i]; j++) {
+                m = atom->map(atom->angle_atom2[i][j]);
+                if (m >= 0 && m < nlocal) ndeltopo[1]++;
+              }
+            }
+          }
+          if (atom->avec->dihedrals_allow) {
+            if (force->newton_bond) ndeltopo[2] += atom->num_dihedral[i];
+            else {
+              for (j = 0; j < atom->num_dihedral[i]; j++) {
+                m = atom->map(atom->dihedral_atom2[i][j]);
+                if (m >= 0 && m < nlocal) ndeltopo[2]++;
+              }
+            }
+          }
+          if (atom->avec->impropers_allow) {
+            if (force->newton_bond) ndeltopo[3] += atom->num_improper[i];
+            else {
+              for (j = 0; j < atom->num_improper[i]; j++) {
+                m = atom->map(atom->improper_atom2[i][j]);
+                if (m >= 0 && m < nlocal) ndeltopo[3]++;
+              }
+            }
+          }
+
+        } else if (me == proc && i == iatom) {
+          mark[i] = 1;
+          ndelone++;
+        }
+      }
+
+      // remove any atoms marked for deletion from my eligible list
+
+      i = 0;
+      while (i < ncount) {
+        if (mark[list[i]]) {
+          list[i] = list[ncount-1];
+          ncount--;
+        } else i++;
+      }
+
+      // update ndel,ncount,nall,nbefore
+      // ndelall is total atoms deleted on this iteration
+      // ncount is already correct, so resum to get nall and nbefore
+
+      MPI_Allreduce(&ndelone,&ndelall,1,MPI_INT,MPI_SUM,world);
+      ndel += ndelall;
+      MPI_Allreduce(&ncount,&nall,1,MPI_INT,MPI_SUM,world);
+      MPI_Scan(&ncount,&nbefore,1,MPI_INT,MPI_SUM,world);
+      nbefore -= ncount;
+    }
+  }
+
+  // delete my marked atoms
+  // loop in reverse order to avoid copying marked atoms
+
+  AtomVec *avec = atom->avec;
+
+  for (i = nlocal-1; i >= 0; i--) {
+    if (mark[i]) {
+      avec->copy(atom->nlocal-1,i,1);
+      atom->nlocal--;
+    }
+  }
+
+  // reset global natoms and bonds, angles, etc
+  // if global map exists, reset it now instead of waiting for comm
+  // since deleting atoms messes up ghosts
+
+  atom->natoms -= ndel;
+  if (molflag) {
+    int all[4];
+    MPI_Allreduce(ndeltopo,all,4,MPI_INT,MPI_SUM,world);
+    atom->nbonds -= all[0];
+    atom->nangles -= all[1];
+    atom->ndihedrals -= all[2];
+    atom->nimpropers -= all[3];
+  }
+
+  if (ndel && atom->map_style) {
+    atom->nghost = 0;
+    atom->map_init();
+    atom->map_set();
+  }
+
+  // statistics
+
+  ndeleted += ndel;
+  next_reneighbor = update->ntimestep + nevery;
+}
+
+/* ----------------------------------------------------------------------
+   return number of deleted particles
+------------------------------------------------------------------------- */
+
+double FixEvaporate::compute_scalar()
+{
+  return 1.0*ndeleted;
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based arrays
+------------------------------------------------------------------------- */
+
+double FixEvaporate::memory_usage()
+{
+  double bytes = 2*nmax * sizeof(int);
+  return bytes;
+}

src/MISC/fix_evaporate.h

@@ -11,34 +11,39 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
-#ifdef COMPUTE_CLASS
+#ifdef FIX_CLASS
 
-ComputeStyle(ti,ComputeTI)
+FixStyle(evaporate,FixEvaporate)
 
 #else
 
-#ifndef COMPUTE_TI_H
+#ifndef LMP_FIX_EVAPORATE_H
-#define COMPUTE_TI_H
+#define LMP_FIX_EVAPORATE_H
 
-#include "compute.h"
+#include "fix.h"
 
 namespace LAMMPS_NS {
 
-class ComputeTI : public Compute {
+class FixEvaporate : public Fix {
  public:
-  ComputeTI(class LAMMPS *, int, char **);
+  FixEvaporate(class LAMMPS *, int, char **);
-  ~ComputeTI();
+  ~FixEvaporate();
+  int setmask();
   void init();
+  void pre_exchange();
   double compute_scalar();
+  double memory_usage();
 
  private:
-  int nterms;
+  int nevery,nflux,iregion;
-  int *which;
+  int molflag;
-  int *ivar1,*ivar2;
+  int ndeleted;
-  int *ilo, *ihi; 
+  char *idregion;
-  char **var1,**var2;
+
-  class Pair **pptr;
+  int nmax;
-  char **pstyle;
+  int *list,*mark;
+
+  class RanPark *random;
 };
 
 }
@@ -54,30 +59,22 @@ Self-explanatory.  Check the input script syntax and compare to the
 documentation for the command.  You can use -echo screen as a
 command-line option when running LAMMPS to see the offending line.
 
-E: Variable name for compute ti does not exist
+E: Region ID for fix evaporate does not exist
 
 Self-explanatory.
 
-E: Variable for compute ti is invalid style
+E: Cannot evaporate atoms in atom_modify first group
 
-Self-explanatory.
+This is a restriction due to the way atoms are organized in
+a list to enable the atom_modify first command.
 
-E: Compute ti pair style does not exist
+W: Fix evaporate may delete atom with non-zero molecule ID
 
-Self-explanatory.
+This is probably an error, since you should not delete only one atom
+of a molecule.
 
-E: Compute ti tail when pair style does not compute tail corrections
+E: Fix evaporate molecule requires atom attribute molecule
 
-Self-explanatory.
+The atom style being used does not define a molecule ID.
-
-E: Compute ti kspace style does not exist
-
-Self-explanatory.
-
-E: Energy was not tallied on needed timestep
-
-You are using a thermo keyword that requires potentials to
-have tallied energy, but they didn't on this timestep.  See the
-variable doc page for ideas on how to make this work.
 
 */

src/MISC/fix_orient_fcc.cpp

@@ -0,0 +1,593 @@
+/* ----------------------------------------------------------------------
+   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 authors: Koenraad Janssens and David Olmsted (SNL)
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "string.h"
+#include "stdlib.h"
+#include "mpi.h"
+#include "fix_orient_fcc.h"
+#include "atom.h"
+#include "update.h"
+#include "respa.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "comm.h"
+#include "output.h"
+#include "math_const.h"
+#include "memory.h"
+#include "error.h"
+#include "force.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+using namespace MathConst;
+
+#define BIG 1000000000
+
+/* ---------------------------------------------------------------------- */
+
+FixOrientFCC::FixOrientFCC(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg)
+{
+  MPI_Comm_rank(world,&me);
+
+  if (narg != 11) error->all(FLERR,"Illegal fix orient/fcc command");
+
+  scalar_flag = 1;
+  global_freq = 1;
+  extscalar = 1;
+
+  peratom_flag = 1;
+  size_peratom_cols = 2;
+  peratom_freq = 1;
+
+  nstats = force->inumeric(FLERR,arg[3]);
+  direction_of_motion = force->inumeric(FLERR,arg[4]);
+  a = force->numeric(FLERR,arg[5]);
+  Vxi = force->numeric(FLERR,arg[6]);
+  uxif_low = force->numeric(FLERR,arg[7]);
+  uxif_high = force->numeric(FLERR,arg[8]);
+
+  if (direction_of_motion == 0) {
+    int n = strlen(arg[9]) + 1;
+    chifilename = new char[n];
+    strcpy(chifilename,arg[9]);
+    n = strlen(arg[10]) + 1;
+    xifilename = new char[n];
+    strcpy(xifilename,arg[10]);
+  } else if (direction_of_motion == 1) {
+    int n = strlen(arg[9]) + 1;
+    xifilename = new char[n];
+    strcpy(xifilename,arg[9]);
+    n = strlen(arg[10]) + 1;
+    chifilename = new char[n];
+    strcpy(chifilename,arg[10]);
+  } else error->all(FLERR,"Illegal fix orient/fcc command");
+
+  // initializations
+
+  half_fcc_nn = 6;
+  use_xismooth = false;
+  double xicutoff = 1.57;
+  xicutoffsq = xicutoff * xicutoff;
+  cutsq = 0.5 * a*a*xicutoffsq;
+  nmax = 0;
+
+  // read xi and chi reference orientations from files
+
+  if (me == 0) {
+    char line[IMGMAX];
+    char *result;
+    int count;
+
+    FILE *infile = fopen(xifilename,"r");
+    if (infile == NULL) error->one(FLERR,"Fix orient/fcc file open failed");
+    for (int i = 0; i < 6; i++) {
+      result = fgets(line,IMGMAX,infile);
+      if (!result) error->one(FLERR,"Fix orient/fcc file read failed");
+      count = sscanf(line,"%lg %lg %lg",&Rxi[i][0],&Rxi[i][1],&Rxi[i][2]);
+      if (count != 3) error->one(FLERR,"Fix orient/fcc file read failed");
+    }
+    fclose(infile);
+
+    infile = fopen(chifilename,"r");
+    if (infile == NULL) error->one(FLERR,"Fix orient/fcc file open failed");
+    for (int i = 0; i < 6; i++) {
+      result = fgets(line,IMGMAX,infile);
+      if (!result) error->one(FLERR,"Fix orient/fcc file read failed");
+      count = sscanf(line,"%lg %lg %lg",&Rchi[i][0],&Rchi[i][1],&Rchi[i][2]);
+      if (count != 3) error->one(FLERR,"Fix orient/fcc file read failed");
+    }
+    fclose(infile);
+  }
+
+  MPI_Bcast(&Rxi[0][0],18,MPI_DOUBLE,0,world);
+  MPI_Bcast(&Rchi[0][0],18,MPI_DOUBLE,0,world);
+
+  // make copy of the reference vectors
+
+  for (int i = 0; i < 6; i++)
+    for (int j = 0; j < 3; j++) {
+      half_xi_chi_vec[0][i][j] = Rxi[i][j];
+      half_xi_chi_vec[1][i][j] = Rchi[i][j];
+    }
+
+  // compute xiid,xi0,xi1 for all 12 neighbors
+  // xi is the favored crystal
+  // want order parameter when actual is Rchi
+
+  double xi_sq,dxi[3],rchi[3];
+
+  xiid = 0.0;
+  for (int i = 0; i < 6; i++) {
+    rchi[0] = Rchi[i][0];
+    rchi[1] = Rchi[i][1];
+    rchi[2] = Rchi[i][2];
+    find_best_ref(rchi,0,xi_sq,dxi);
+    xiid += sqrt(xi_sq);
+    for (int j = 0; j < 3; j++) rchi[j] = -rchi[j];
+    find_best_ref(rchi,0,xi_sq,dxi);
+    xiid += sqrt(xi_sq);
+  }
+
+  xiid /= 12.0;
+  xi0 = uxif_low * xiid;
+  xi1 = uxif_high * xiid;
+
+  // set comm size needed by this Fix
+  // NOTE: doesn't seem that use_xismooth is ever true
+
+  if (use_xismooth) comm_forward = 62;
+  else comm_forward = 50;
+
+  added_energy = 0.0;
+
+  nmax = atom->nmax;
+  nbr = (Nbr *) memory->smalloc(nmax*sizeof(Nbr),"orient/fcc:nbr");
+  memory->create(order,nmax,2,"orient/fcc:order");
+  array_atom = order;
+
+  // zero the array since a variable may access it before first run
+
+  int nlocal = atom->nlocal;
+  for (int i = 0; i < nlocal; i++) order[i][0] = order[i][1] = 0.0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixOrientFCC::~FixOrientFCC()
+{
+  delete [] xifilename;
+  delete [] chifilename;
+  memory->sfree(nbr);
+  memory->destroy(order);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixOrientFCC::setmask()
+{
+  int mask = 0;
+  mask |= POST_FORCE;
+  mask |= THERMO_ENERGY;
+  mask |= POST_FORCE_RESPA;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixOrientFCC::init()
+{
+  if (strstr(update->integrate_style,"respa"))
+    nlevels_respa = ((Respa *) update->integrate)->nlevels;
+
+  // need a full neighbor list, built whenever re-neighboring occurs
+
+  int irequest = neighbor->request((void *) this);
+  neighbor->requests[irequest]->pair = 0;
+  neighbor->requests[irequest]->fix = 1;
+  neighbor->requests[irequest]->half = 0;
+  neighbor->requests[irequest]->full = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixOrientFCC::init_list(int id, NeighList *ptr)
+{
+  list = ptr;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixOrientFCC::setup(int vflag)
+{
+  if (strstr(update->integrate_style,"verlet"))
+    post_force(vflag);
+  else {
+    ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
+    post_force_respa(vflag,nlevels_respa-1,0);
+    ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixOrientFCC::post_force(int vflag)
+{
+  int i,j,k,ii,jj,inum,jnum,m,n,nn,nsort,id_self;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+  double edelta,omega;
+  double dx,dy,dz,rsq,xismooth,xi_sq,duxi,duxi_other;
+  double dxi[3];
+  double *dxiptr;
+  bool found_myself;
+
+  // set local ptrs
+
+  double **x = atom->x;
+  double **f = atom->f;
+  int *mask = atom->mask;
+  int *tag = atom->tag;
+  int nlocal = atom->nlocal;
+  int nall = atom->nlocal + atom->nghost;
+
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  // insure nbr and order data structures are adequate size
+
+  if (nall > nmax) {
+    nmax = nall;
+    memory->destroy(nbr);
+    memory->destroy(order);
+    nbr = (Nbr *) memory->smalloc(nmax*sizeof(Nbr),"orient/fcc:nbr");
+    memory->create(order,nmax,2,"orient/fcc:order");
+    array_atom = order;
+  }
+
+  // loop over owned atoms and build Nbr data structure of neighbors
+  // use full neighbor list
+
+  added_energy = 0.0;
+  int count = 0;
+  int mincount = BIG;
+  int maxcount = 0;
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+
+    if (jnum < mincount) mincount = jnum;
+    if (jnum > maxcount) {
+      if (maxcount) delete [] sort;
+      sort = new Sort[jnum];
+      maxcount = jnum;
+    }
+
+    // loop over all neighbors of atom i
+    // for those within cutsq, build sort data structure
+    // store local id, rsq, delta vector, xismooth (if included)
+
+    nsort = 0;
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      j &= NEIGHMASK;
+      count++;
+
+      dx = x[i][0] - x[j][0];
+      dy = x[i][1] - x[j][1];
+      dz = x[i][2] - x[j][2];
+      rsq = dx*dx + dy*dy + dz*dz;
+
+      if (rsq < cutsq) {
+        sort[nsort].id = j;
+        sort[nsort].rsq = rsq;
+        sort[nsort].delta[0] = dx;
+        sort[nsort].delta[1] = dy;
+        sort[nsort].delta[2] = dz;
+        if (use_xismooth) {
+          xismooth = (xicutoffsq - 2.0*rsq/(a*a)) / (xicutoffsq - 1.0);
+          sort[nsort].xismooth = 1.0 - fabs(1.0-xismooth);
+        }
+        nsort++;
+      }
+    }
+
+    // sort neighbors by rsq distance
+    // no need to sort if nsort <= 12
+
+    if (nsort > 12) qsort(sort,nsort,sizeof(Sort),compare);
+
+    // copy up to 12 nearest neighbors into nbr data structure
+    // operate on delta vector via find_best_ref() to compute dxi
+
+    n = MIN(12,nsort);
+    nbr[i].n = n;
+    if (n == 0) continue;
+
+    double xi_total = 0.0;
+    for (j = 0; j < n; j++) {
+      find_best_ref(sort[j].delta,0,xi_sq,dxi);
+      xi_total += sqrt(xi_sq);
+      nbr[i].id[j] = sort[j].id;
+      nbr[i].dxi[j][0] = dxi[0]/n;
+      nbr[i].dxi[j][1] = dxi[1]/n;
+      nbr[i].dxi[j][2] = dxi[2]/n;
+      if (use_xismooth) nbr[i].xismooth[j] = sort[j].xismooth;
+    }
+    xi_total /= n;
+    order[i][0] = xi_total;
+
+    // compute potential derivative to xi
+
+    if (xi_total < xi0) {
+      nbr[i].duxi = 0.0;
+      edelta = 0.0;
+      order[i][1] = 0.0;
+    } else if (xi_total > xi1) {
+      nbr[i].duxi = 0.0;
+      edelta = Vxi;
+      order[i][1] = 1.0;
+    } else {
+      omega = MY_PI2*(xi_total-xi0) / (xi1-xi0);
+      nbr[i].duxi = MY_PI*Vxi*sin(2.0*omega) / (2.0*(xi1-xi0));
+      edelta = Vxi*(1 - cos(2.0*omega)) / 2.0;
+      order[i][1] = omega / MY_PI2;
+    }
+    added_energy += edelta;
+  }
+
+  if (maxcount) delete [] sort;
+
+  // communicate to acquire nbr data for ghost atoms
+
+  comm->forward_comm_fix(this);
+
+  // compute grain boundary force on each owned atom
+  // skip atoms not in group
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    if (!(mask[i] & groupbit)) continue;
+    n = nbr[i].n;
+    duxi = nbr[i].duxi;
+
+    for (j = 0; j < n; j++) {
+      dxiptr = &nbr[i].dxi[j][0];
+      if (use_xismooth) {
+        xismooth = nbr[i].xismooth[j];
+        f[i][0] += duxi * dxiptr[0] * xismooth;
+        f[i][1] += duxi * dxiptr[1] * xismooth;
+        f[i][2] += duxi * dxiptr[2] * xismooth;
+      } else {
+        f[i][0] += duxi * dxiptr[0];
+        f[i][1] += duxi * dxiptr[1];
+        f[i][2] += duxi * dxiptr[2];
+      }
+
+      // m = local index of neighbor
+      // id_self = ID for atom I in atom M's neighbor list
+      // if M is local atom, id_self will be local ID of atom I
+      // if M is ghost atom, id_self will be global ID of atom I
+
+      m = nbr[i].id[j];
+      if (m < nlocal) id_self = i;
+      else id_self = tag[i];
+      found_myself = false;
+      nn = nbr[m].n;
+
+      for (k = 0; k < nn; k++) {
+        if (id_self == nbr[m].id[k]) {
+          if (found_myself) error->one(FLERR,"Fix orient/fcc found self twice");
+          found_myself = true;
+          duxi_other = nbr[m].duxi;
+          dxiptr = &nbr[m].dxi[k][0];
+          if (use_xismooth) {
+            xismooth = nbr[m].xismooth[k];
+            f[i][0] -= duxi_other * dxiptr[0] * xismooth;
+            f[i][1] -= duxi_other * dxiptr[1] * xismooth;
+            f[i][2] -= duxi_other * dxiptr[2] * xismooth;
+          } else {
+            f[i][0] -= duxi_other * dxiptr[0];
+            f[i][1] -= duxi_other * dxiptr[1];
+            f[i][2] -= duxi_other * dxiptr[2];
+          }
+        }
+      }
+    }
+  }
+
+  // print statistics every nstats timesteps
+
+  if (nstats && update->ntimestep % nstats == 0) {
+    int total;
+    MPI_Allreduce(&count,&total,1,MPI_INT,MPI_SUM,world);
+    double ave = total/atom->natoms;
+
+    int min,max;
+    MPI_Allreduce(&mincount,&min,1,MPI_INT,MPI_MIN,world);
+    MPI_Allreduce(&maxcount,&max,1,MPI_INT,MPI_MAX,world);
+
+    if (me == 0) {
+      if (screen) fprintf(screen,
+                          "orient step " BIGINT_FORMAT ": " BIGINT_FORMAT
+                          " atoms have %d neighbors\n",
+                          update->ntimestep,atom->natoms,total);
+      if (logfile) fprintf(logfile,
+                           "orient step " BIGINT_FORMAT ": " BIGINT_FORMAT
+                           " atoms have %d neighbors\n",
+                           update->ntimestep,atom->natoms,total);
+      if (screen)
+        fprintf(screen,"  neighs: min = %d, max = %d, ave = %g\n",
+                min,max,ave);
+      if (logfile)
+        fprintf(logfile,"  neighs: min = %d, max = %d, ave = %g\n",
+                min,max,ave);
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixOrientFCC::post_force_respa(int vflag, int ilevel, int iloop)
+{
+  if (ilevel == nlevels_respa-1) post_force(vflag);
+}
+
+/* ---------------------------------------------------------------------- */
+
+double FixOrientFCC::compute_scalar()
+{
+  double added_energy_total;
+  MPI_Allreduce(&added_energy,&added_energy_total,1,MPI_DOUBLE,MPI_SUM,world);
+  return added_energy_total;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixOrientFCC::pack_comm(int n, int *list, double *buf,
+                            int pbc_flag, int *pbc)
+{
+  int i,j,k,id,num;
+  int *tag = atom->tag;
+  int nlocal = atom->nlocal;
+  int m = 0;
+
+  for (i = 0; i < n; i++) {
+    k = list[i];
+    num = nbr[k].n;
+    buf[m++] = num;
+    buf[m++] = nbr[k].duxi;
+
+    for (j = 0; j < num; j++) {
+      if (use_xismooth) buf[m++] = nbr[m].xismooth[j];
+      buf[m++] = nbr[k].dxi[j][0];
+      buf[m++] = nbr[k].dxi[j][1];
+      buf[m++] = nbr[k].dxi[j][2];
+
+      // id stored in buf needs to be global ID
+      // if k is a local atom, it stores local IDs, so convert to global
+      // if k is a ghost atom (already comm'd), its IDs are already global
+
+      id = nbr[k].id[j];
+      if (k < nlocal) id = tag[id];
+      buf[m++] = id;
+    }
+
+    m += (12-num) * 3;
+    if (use_xismooth) m += 12-num;
+  }
+
+  if (use_xismooth) return 62;
+  return 50;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixOrientFCC::unpack_comm(int n, int first, double *buf)
+{
+  int i,j,num;
+  int last = first + n;
+  int m = 0;
+
+  for (i = first; i < last; i++) {
+    nbr[i].n = num = static_cast<int> (buf[m++]);
+    nbr[i].duxi = buf[m++];
+
+    for (j = 0; j < num; j++) {
+      if (use_xismooth) nbr[i].xismooth[j] = buf[m++];
+      nbr[i].dxi[j][0] = buf[m++];
+      nbr[i].dxi[j][1] = buf[m++];
+      nbr[i].dxi[j][2] = buf[m++];
+      nbr[i].id[j] = static_cast<int> (buf[m++]);
+    }
+
+    m += (12-num) * 3;
+    if (use_xismooth) m += 12-num;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixOrientFCC::find_best_ref(double *displs, int which_crystal,
+                                 double &xi_sq, double *dxi)
+{
+  int i;
+  double dot,tmp;
+
+  double  best_dot  = -1.0;         // best is biggest (smallest angle)
+  int     best_i    = -1;
+  int     best_sign = 0;
+
+  for (i = 0; i < half_fcc_nn; i++) {
+    dot = displs[0] * half_xi_chi_vec[which_crystal][i][0] +
+      displs[1] * half_xi_chi_vec[which_crystal][i][1] +
+      displs[2] * half_xi_chi_vec[which_crystal][i][2];
+    if (fabs(dot) > best_dot) {
+      best_dot = fabs(dot);
+      best_i = i;
+      if (dot < 0.0) best_sign = -1;
+      else best_sign = 1;
+    }
+  }
+
+  xi_sq = 0.0;
+  for (i = 0; i < 3; i++) {
+    tmp = displs[i] - best_sign * half_xi_chi_vec[which_crystal][best_i][i];
+    xi_sq += tmp*tmp;
+  }
+
+  if (xi_sq > 0.0) {
+    double xi = sqrt(xi_sq);
+    for (i = 0; i < 3; i++)
+      dxi[i] = (best_sign * half_xi_chi_vec[which_crystal][best_i][i] -
+                displs[i]) / xi;
+  } else dxi[0] = dxi[1] = dxi[2] = 0.0;
+}
+
+/* ----------------------------------------------------------------------
+   compare two neighbors I and J in sort data structure
+   called via qsort in post_force() method
+   is a static method so can't access sort data structure directly
+   return -1 if I < J, 0 if I = J, 1 if I > J
+   do comparison based on rsq distance
+------------------------------------------------------------------------- */
+
+int FixOrientFCC::compare(const void *pi, const void *pj)
+{
+  FixOrientFCC::Sort *ineigh = (FixOrientFCC::Sort *) pi;
+  FixOrientFCC::Sort *jneigh = (FixOrientFCC::Sort *) pj;
+
+  if (ineigh->rsq < jneigh->rsq) return -1;
+  else if (ineigh->rsq > jneigh->rsq) return 1;
+  return 0;
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based arrays
+------------------------------------------------------------------------- */
+
+double FixOrientFCC::memory_usage()
+{
+  double bytes = nmax * sizeof(Nbr);
+  bytes += 2*nmax * sizeof(double);
+  return bytes;
+}

src/MISC/fix_orient_fcc.h

@@ -0,0 +1,115 @@
+/* ----------------------------------------------------------------------
+   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(orient/fcc,FixOrientFCC)
+
+#else
+
+#ifndef LMP_FIX_ORIENT_FCC_H
+#define LMP_FIX_ORIENT_FCC_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixOrientFCC : public Fix {
+ public:
+  struct Nbr {              // neighbor info for each owned and ghost atom
+    int n;                  // # of closest neighbors (up to 12)
+    int id[12];             // IDs of each neighbor
+                            // if center atom is owned, these are local IDs
+                            // if center atom is ghost, these are global IDs
+    double xismooth[12];    // distance weighting factor for each neighbors
+    double dxi[12][3];      // d order-parameter / dx for each neighbor
+    double duxi;            // d Energy / d order-parameter for atom
+  };
+
+  struct Sort {             // data structure for sorting to find 12 closest
+    int id;                 // ID of neighbor atom
+    double rsq;             // distance between center and neighbor atom
+    double delta[3];        // displacement between center and neighbor atom
+    double xismooth;        // distance weighting factor
+  };
+
+  FixOrientFCC(class LAMMPS *, int, char **);
+  ~FixOrientFCC();
+  int setmask();
+  void init();
+  void init_list(int, class NeighList *);
+  void setup(int);
+  void post_force(int);
+  void post_force_respa(int, int, int);
+  double compute_scalar();
+  int pack_comm(int, int *, double *, int, int *);
+  void unpack_comm(int, int, double *);
+  double memory_usage();
+
+ private:
+  int me;
+  int nlevels_respa;
+
+  int direction_of_motion;         // 1 = center shrinks, 0 = center grows
+  int nstats;                      // stats output every this many steps
+  double a;                        // lattice parameter
+  double Vxi;                      // potential value
+  double uxif_low;                 // cut-off fraction, low order parameter
+  double uxif_high;                // cut-off fraction, high order parameter
+  char *xifilename, *chifilename;  // file names for 2 crystal orientations
+
+  bool use_xismooth;
+  double Rxi[12][3],Rchi[12][3],half_xi_chi_vec[2][6][3];
+  double xiid,xi0,xi1,xicutoffsq,cutsq,added_energy;
+  int half_fcc_nn;
+
+  int nmax;                        // expose 2 per-atom quantities
+  double **order;                  // order param and normalized order param
+
+  Nbr *nbr;
+  Sort *sort;
+  class NeighList *list;
+
+  void find_best_ref(double *, int, double &, double *);
+  static int compare(const void *, const void *);
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Fix orient/fcc file open failed
+
+The fix orient/fcc command could not open a specified file.
+
+E: Fix orient/fcc file read failed
+
+The fix orient/fcc command could not read the needed parameters from a
+specified file.
+
+E: Fix orient/fcc found self twice
+
+The neighbor lists used by fix orient/fcc are messed up.  If this
+error occurs, it is likely a bug, so send an email to the
+"developers"_http://lammps.sandia.gov/authors.html.
+
+*/

src/MISC/fix_thermal_conductivity.cpp

@@ -0,0 +1,330 @@
+/* ----------------------------------------------------------------------
+   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: Craig Tenney (UND) added support
+                        for swapping atoms of different masses
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "mpi.h"
+#include "string.h"
+#include "stdlib.h"
+#include "fix_thermal_conductivity.h"
+#include "atom.h"
+#include "force.h"
+#include "domain.h"
+#include "modify.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+#define BIG 1.0e10
+
+/* ---------------------------------------------------------------------- */
+
+FixThermalConductivity::FixThermalConductivity(LAMMPS *lmp,
+                                               int narg, char **arg) :
+  Fix(lmp, narg, arg)
+{
+  if (narg < 6) error->all(FLERR,"Illegal fix thermal/conductivity command");
+
+  MPI_Comm_rank(world,&me);
+
+  nevery = force->inumeric(FLERR,arg[3]);
+  if (nevery <= 0) error->all(FLERR,"Illegal fix thermal/conductivity command");
+
+  scalar_flag = 1;
+  global_freq = nevery;
+  extscalar = 0;
+
+  if (strcmp(arg[4],"x") == 0) edim = 0;
+  else if (strcmp(arg[4],"y") == 0) edim = 1;
+  else if (strcmp(arg[4],"z") == 0) edim = 2;
+  else error->all(FLERR,"Illegal fix thermal/conductivity command");
+
+  nbin = force->inumeric(FLERR,arg[5]);
+  if (nbin % 2 || nbin <= 2)
+    error->all(FLERR,"Illegal fix thermal/conductivity command");
+
+  // optional keywords
+
+  nswap = 1;
+
+  int iarg = 6;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"swap") == 0) {
+      if (iarg+2 > narg)
+        error->all(FLERR,"Illegal fix thermal/conductivity command");
+      nswap = force->inumeric(FLERR,arg[iarg+1]);
+      if (nswap <= 0)
+        error->all(FLERR,
+                   "Fix thermal/conductivity swap value must be positive");
+      iarg += 2;
+    } else error->all(FLERR,"Illegal fix thermal/conductivity command");
+  }
+
+  // initialize array sizes to nswap+1 so have space to shift values down
+
+  index_lo = new int[nswap+1];
+  index_hi = new int[nswap+1];
+  ke_lo = new double[nswap+1];
+  ke_hi = new double[nswap+1];
+
+  e_exchange = 0.0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixThermalConductivity::~FixThermalConductivity()
+{
+  delete [] index_lo;
+  delete [] index_hi;
+  delete [] ke_lo;
+  delete [] ke_hi;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixThermalConductivity::setmask()
+{
+  int mask = 0;
+  mask |= END_OF_STEP;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixThermalConductivity::init()
+{
+  // warn if any fix ave/spatial comes after this fix
+  // can cause glitch in averaging since ave will happen after swap
+
+  int foundme = 0;
+  for (int i = 0; i < modify->nfix; i++) {
+    if (modify->fix[i] == this) foundme = 1;
+    if (foundme && strcmp(modify->fix[i]->style,"ave/spatial") == 0 && me == 0)
+      error->warning(FLERR,
+                     "Fix thermal/conductivity comes before fix ave/spatial");
+  }
+
+  // set bounds of 2 slabs in edim
+  // only necessary for static box, else re-computed in end_of_step()
+  // lo bin is always bottom bin
+  // hi bin is just above half height
+
+  if (domain->box_change == 0) {
+    prd = domain->prd[edim];
+    boxlo = domain->boxlo[edim];
+    boxhi = domain->boxhi[edim];
+    double binsize = (boxhi-boxlo) / nbin;
+    slablo_lo = boxlo;
+    slablo_hi = boxlo + binsize;
+    slabhi_lo = boxlo + (nbin/2)*binsize;
+    slabhi_hi = boxlo + (nbin/2+1)*binsize;
+  }
+
+  periodicity = domain->periodicity[edim];
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixThermalConductivity::end_of_step()
+{
+  int i,j,m,insert;
+  double coord,ke;
+  MPI_Status status;
+  struct {
+    double value;
+    int proc;
+  } mine[2],all[2];
+
+  // if box changes, recompute bounds of 2 slabs in edim
+
+  if (domain->box_change) {
+    prd = domain->prd[edim];
+    boxlo = domain->boxlo[edim];
+    boxhi = domain->boxhi[edim];
+    double binsize = (boxhi-boxlo) / nbin;
+    slablo_lo = boxlo;
+    slablo_hi = boxlo + binsize;
+    slabhi_lo = boxlo + (nbin/2)*binsize;
+    slabhi_hi = boxlo + (nbin/2+1)*binsize;
+  }
+
+  // make 2 lists of up to nswap atoms
+  // hottest atoms in lo slab, coldest atoms in hi slab (really mid slab)
+  // lo slab list is sorted by hottest, hi slab is sorted by coldest
+  // map atoms back into periodic box if necessary
+  // insert = location in list to insert new atom
+
+  double **x = atom->x;
+  double **v = atom->v;
+  double *mass = atom->mass;
+  double *rmass = atom->rmass;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  nhi = nlo = 0;
+
+  for (i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      coord = x[i][edim];
+      if (coord < boxlo && periodicity) coord += prd;
+      else if (coord >= boxhi && periodicity) coord -= prd;
+
+      if (coord >= slablo_lo && coord < slablo_hi) {
+        ke = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
+        if (rmass) ke *= 0.5*rmass[i];
+        else ke *= 0.5*mass[type[i]];
+        if (nlo < nswap || ke > ke_lo[nswap-1]) {
+          for (insert = nlo-1; insert >= 0; insert--)
+            if (ke < ke_lo[insert]) break;
+          insert++;
+          for (m = nlo-1; m >= insert; m--) {
+            ke_lo[m+1] = ke_lo[m];
+            index_lo[m+1] = index_lo[m];
+          }
+          ke_lo[insert] = ke;
+          index_lo[insert] = i;
+          if (nlo < nswap) nlo++;
+        }
+      }
+
+      if (coord >= slabhi_lo && coord < slabhi_hi) {
+        ke = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
+        if (rmass) ke *= 0.5*rmass[i];
+        else ke *= 0.5*mass[type[i]];
+        if (nhi < nswap || ke < ke_hi[nswap-1]) {
+          for (insert = nhi-1; insert >= 0; insert--)
+            if (ke > ke_hi[insert]) break;
+          insert++;
+          for (m = nhi-1; m >= insert; m--) {
+            ke_hi[m+1] = ke_hi[m];
+            index_hi[m+1] = index_hi[m];
+          }
+          ke_hi[insert] = ke;
+          index_hi[insert] = i;
+          if (nhi < nswap) nhi++;
+        }
+      }
+    }
+
+  // loop over nswap pairs
+  // pair 2 global atoms at beginning of sorted lo/hi slab lists via Allreduce
+  // BIG values are for procs with no atom to contribute
+  // use negative of hottest KE since is doing a MINLOC
+  // MINLOC also communicates which procs own them
+  // exchange kinetic energy between the 2 particles
+  // if I own both particles just swap, else point2point comm of velocities
+
+  double sbuf[4],rbuf[4],vcm[3];
+  double eswap = 0.0;
+
+  mine[0].proc = mine[1].proc = me;
+  int ilo = 0;
+  int ihi = 0;
+
+  for (m = 0; m < nswap; m++) {
+    if (ilo < nlo) mine[0].value = -ke_lo[ilo];
+    else mine[0].value = BIG;
+    if (ihi < nhi) mine[1].value = ke_hi[ihi];
+    else mine[1].value = BIG;
+
+    MPI_Allreduce(mine,all,2,MPI_DOUBLE_INT,MPI_MINLOC,world);
+    if (all[0].value == BIG || all[1].value == BIG) continue;
+
+    if (me == all[0].proc && me == all[1].proc) {
+      i = index_lo[ilo++];
+      j = index_hi[ihi++];
+      sbuf[0] = v[j][0];
+      sbuf[1] = v[j][1];
+      sbuf[2] = v[j][2];
+      if (rmass) sbuf[3] = rmass[j];
+      else sbuf[3] = mass[type[j]];
+      rbuf[0] = v[i][0];
+      rbuf[1] = v[i][1];
+      rbuf[2] = v[i][2];
+      if (rmass) rbuf[3] = rmass[i];
+      else rbuf[3] = mass[type[i]];
+      vcm[0] = (sbuf[3]*sbuf[0] + rbuf[3]*rbuf[0]) / (sbuf[3] + rbuf[3]);
+      vcm[1] = (sbuf[3]*sbuf[1] + rbuf[3]*rbuf[1]) / (sbuf[3] + rbuf[3]);
+      vcm[2] = (sbuf[3]*sbuf[2] + rbuf[3]*rbuf[2]) / (sbuf[3] + rbuf[3]);
+      v[j][0] = 2.0 * vcm[0] - sbuf[0];
+      v[j][1] = 2.0 * vcm[1] - sbuf[1];
+      v[j][2] = 2.0 * vcm[2] - sbuf[2];
+      eswap += sbuf[3] * (vcm[0] * (vcm[0] - sbuf[0]) +
+                          vcm[1] * (vcm[1] - sbuf[1]) +
+                          vcm[2] * (vcm[2] - sbuf[2]));
+      v[i][0] = 2.0 * vcm[0] - rbuf[0];
+      v[i][1] = 2.0 * vcm[1] - rbuf[1];
+      v[i][2] = 2.0 * vcm[2] - rbuf[2];
+      eswap -= rbuf[3] * (vcm[0] * (vcm[0] - rbuf[0]) +
+                          vcm[1] * (vcm[1] - rbuf[1]) +
+                          vcm[2] * (vcm[2] - rbuf[2]));
+
+    } else if (me == all[0].proc) {
+      j = index_lo[ilo++];
+      sbuf[0] = v[j][0];
+      sbuf[1] = v[j][1];
+      sbuf[2] = v[j][2];
+      if (rmass) sbuf[3] = rmass[j];
+      else sbuf[3] = mass[type[j]];
+      MPI_Sendrecv(sbuf,4,MPI_DOUBLE,all[1].proc,0,
+                   rbuf,4,MPI_DOUBLE,all[1].proc,0,world,&status);
+      vcm[0] = (sbuf[3]*sbuf[0] + rbuf[3]*rbuf[0]) / (sbuf[3] + rbuf[3]);
+      vcm[1] = (sbuf[3]*sbuf[1] + rbuf[3]*rbuf[1]) / (sbuf[3] + rbuf[3]);
+      vcm[2] = (sbuf[3]*sbuf[2] + rbuf[3]*rbuf[2]) / (sbuf[3] + rbuf[3]);
+      v[j][0] = 2.0 * vcm[0] - sbuf[0];
+      v[j][1] = 2.0 * vcm[1] - sbuf[1];
+      v[j][2] = 2.0 * vcm[2] - sbuf[2];
+      eswap -= sbuf[3] * (vcm[0] * (vcm[0] - sbuf[0]) +
+                          vcm[1] * (vcm[1] - sbuf[1]) +
+                          vcm[2] * (vcm[2] - sbuf[2]));
+
+    } else if (me == all[1].proc) {
+      j = index_hi[ihi++];
+      sbuf[0] = v[j][0];
+      sbuf[1] = v[j][1];
+      sbuf[2] = v[j][2];
+      if (rmass) sbuf[3] = rmass[j];
+      else sbuf[3] = mass[type[j]];
+      MPI_Sendrecv(sbuf,4,MPI_DOUBLE,all[0].proc,0,
+                   rbuf,4,MPI_DOUBLE,all[0].proc,0,world,&status);
+      vcm[0] = (sbuf[3]*sbuf[0] + rbuf[3]*rbuf[0]) / (sbuf[3] + rbuf[3]);
+      vcm[1] = (sbuf[3]*sbuf[1] + rbuf[3]*rbuf[1]) / (sbuf[3] + rbuf[3]);
+      vcm[2] = (sbuf[3]*sbuf[2] + rbuf[3]*rbuf[2]) / (sbuf[3] + rbuf[3]);
+      v[j][0] = 2.0 * vcm[0] - sbuf[0];
+      v[j][1] = 2.0 * vcm[1] - sbuf[1];
+      v[j][2] = 2.0 * vcm[2] - sbuf[2];
+      eswap += sbuf[3] * (vcm[0] * (vcm[0] - sbuf[0]) +
+                          vcm[1] * (vcm[1] - sbuf[1]) +
+                          vcm[2] * (vcm[2] - sbuf[2]));
+    }
+  }
+
+  // tally energy exchange from all swaps
+
+  double eswap_all;
+  MPI_Allreduce(&eswap,&eswap_all,1,MPI_DOUBLE,MPI_SUM,world);
+  e_exchange += force->mvv2e * eswap_all;
+}
+
+/* ---------------------------------------------------------------------- */
+
+double FixThermalConductivity::compute_scalar()
+{
+  return e_exchange;
+}

src/MISC/fix_thermal_conductivity.h

@@ -0,0 +1,75 @@
+/* ----------------------------------------------------------------------
+   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(thermal/conductivity,FixThermalConductivity)
+
+#else
+
+#ifndef LMP_FIX_THERMAL_CONDUCTIVITY_H
+#define LMP_FIX_THERMAL_CONDUCTIVITY_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixThermalConductivity : public Fix {
+ public:
+  FixThermalConductivity(class LAMMPS *, int, char **);
+  ~FixThermalConductivity();
+  int setmask();
+  void init();
+  void end_of_step();
+  double compute_scalar();
+
+ private:
+  int me;
+  int edim,nbin,periodicity;
+  int nswap;
+  double prd,boxlo,boxhi;
+  double slablo_lo,slablo_hi,slabhi_lo,slabhi_hi;
+  double e_exchange;
+
+  int nlo,nhi;
+  int *index_lo,*index_hi;
+  double *ke_lo,*ke_hi;
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Fix thermal/conductivity swap value must be positive
+
+Self-explanatory.
+
+W: Fix thermal/conductivity comes before fix ave/spatial
+
+The order of these 2 fixes in your input script is such that fix
+thermal/conductivity comes first.  If you are using fix ave/spatial to
+measure the temperature profile induced by fix viscosity, then this
+may cause a glitch in the profile since you are averaging immediately
+after swaps have occurred.  Flipping the order of the 2 fixes
+typically helps.
+
+*/

src/MISC/fix_ttm.cpp

@@ -0,0 +1,685 @@
+/* ----------------------------------------------------------------------
+   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 authors: Paul Crozier (SNL)
+                         Carolyn Phillips (University of Michigan)
+------------------------------------------------------------------------- */
+
+#include "lmptype.h"
+#include "mpi.h"
+#include "math.h"
+#include "string.h"
+#include "stdlib.h"
+#include "fix_ttm.h"
+#include "atom.h"
+#include "force.h"
+#include "update.h"
+#include "domain.h"
+#include "region.h"
+#include "respa.h"
+#include "comm.h"
+#include "random_mars.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+#define MAXLINE 1024
+
+/* ---------------------------------------------------------------------- */
+
+FixTTM::FixTTM(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg)
+{
+  if (narg < 15) error->all(FLERR,"Illegal fix ttm command");
+
+  vector_flag = 1;
+  size_vector = 2;
+  global_freq = 1;
+  extvector = 1;
+  nevery = 1;
+  restart_peratom = 1;
+  restart_global = 1;
+
+  seed = force->inumeric(FLERR,arg[3]);
+  electronic_specific_heat = force->numeric(FLERR,arg[4]);
+  electronic_density = force->numeric(FLERR,arg[5]);
+  electronic_thermal_conductivity = force->numeric(FLERR,arg[6]);
+  gamma_p = force->numeric(FLERR,arg[7]);
+  gamma_s = force->numeric(FLERR,arg[8]);
+  v_0 = force->numeric(FLERR,arg[9]);
+  nxnodes = force->inumeric(FLERR,arg[10]);
+  nynodes = force->inumeric(FLERR,arg[11]);
+  nznodes = force->inumeric(FLERR,arg[12]);
+
+  fpr = fopen(arg[13],"r");
+  if (fpr == NULL) {
+    char str[128];
+    sprintf(str,"Cannot open file %s",arg[13]);
+    error->one(FLERR,str);
+  }
+
+  nfileevery = force->inumeric(FLERR,arg[14]);
+
+  if (nfileevery) {
+    if (narg != 16) error->all(FLERR,"Illegal fix ttm command");
+    MPI_Comm_rank(world,&me);
+    if (me == 0) {
+      fp = fopen(arg[15],"w");
+      if (fp == NULL) {
+        char str[128];
+        sprintf(str,"Cannot open fix ttm file %s",arg[15]);
+        error->one(FLERR,str);
+      }
+    }
+  }
+
+  // error check
+
+  if (seed <= 0) error->all(FLERR,"Invalid random number seed in fix ttm command");
+  if (electronic_specific_heat <= 0.0)
+    error->all(FLERR,"Fix ttm electronic_specific_heat must be > 0.0");
+  if (electronic_density <= 0.0)
+    error->all(FLERR,"Fix ttm electronic_density must be > 0.0");
+  if (electronic_thermal_conductivity < 0.0)
+    error->all(FLERR,"Fix ttm electronic_thermal_conductivity must be >= 0.0");
+  if (gamma_p <= 0.0) error->all(FLERR,"Fix ttm gamma_p must be > 0.0");
+  if (gamma_s < 0.0) error->all(FLERR,"Fix ttm gamma_s must be >= 0.0");
+  if (v_0 < 0.0) error->all(FLERR,"Fix ttm v_0 must be >= 0.0");
+  if (nxnodes <= 0 || nynodes <= 0 || nznodes <= 0)
+    error->all(FLERR,"Fix ttm number of nodes must be > 0");
+
+  v_0_sq = v_0*v_0;
+
+  // initialize Marsaglia RNG with processor-unique seed
+
+  random = new RanMars(lmp,seed + comm->me);
+
+  // allocate per-type arrays for force prefactors
+
+  gfactor1 = new double[atom->ntypes+1];
+  gfactor2 = new double[atom->ntypes+1];
+
+  // allocate 3d grid variables
+
+  total_nnodes = nxnodes*nynodes*nznodes;
+
+  memory->create(nsum,nxnodes,nynodes,nznodes,"ttm:nsum");
+  memory->create(nsum_all,nxnodes,nynodes,nznodes,"ttm:nsum_all");
+  memory->create(T_initial_set,nxnodes,nynodes,nznodes,"ttm:T_initial_set");
+  memory->create(sum_vsq,nxnodes,nynodes,nznodes,"ttm:sum_vsq");
+  memory->create(sum_mass_vsq,nxnodes,nynodes,nznodes,"ttm:sum_mass_vsq");
+  memory->create(sum_vsq_all,nxnodes,nynodes,nznodes,"ttm:sum_vsq_all");
+  memory->create(sum_mass_vsq_all,nxnodes,nynodes,nznodes,
+                 "ttm:sum_mass_vsq_all");
+  memory->create(T_electron_old,nxnodes,nynodes,nznodes,"ttm:T_electron_old");
+  memory->create(T_electron,nxnodes,nynodes,nznodes,"ttm:T_electron");
+  memory->create(net_energy_transfer,nxnodes,nynodes,nznodes,
+                 "TTM:net_energy_transfer");
+  memory->create(net_energy_transfer_all,nxnodes,nynodes,nznodes,
+                 "TTM:net_energy_transfer_all");
+
+  flangevin = NULL;
+  grow_arrays(atom->nmax);
+
+  // zero out the flangevin array
+
+  for (int i = 0; i < atom->nmax; i++) {
+    flangevin[i][0] = 0;
+    flangevin[i][1] = 0;
+    flangevin[i][2] = 0;
+  }
+
+  atom->add_callback(0);
+  atom->add_callback(1);
+
+  // set initial electron temperatures from user input file
+
+  if (me == 0) read_initial_electron_temperatures();
+  MPI_Bcast(&T_electron[0][0][0],total_nnodes,MPI_DOUBLE,0,world);
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixTTM::~FixTTM()
+{
+  if (nfileevery && me == 0) fclose(fp);
+
+  delete random;
+
+  delete [] gfactor1;
+  delete [] gfactor2;
+
+  memory->destroy(nsum);
+  memory->destroy(nsum_all);
+  memory->destroy(T_initial_set);
+  memory->destroy(sum_vsq);
+  memory->destroy(sum_mass_vsq);
+  memory->destroy(sum_vsq_all);
+  memory->destroy(sum_mass_vsq_all);
+  memory->destroy(T_electron_old);
+  memory->destroy(T_electron);
+  memory->destroy(flangevin);
+  memory->destroy(net_energy_transfer);
+  memory->destroy(net_energy_transfer_all);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixTTM::setmask()
+{
+  int mask = 0;
+  mask |= POST_FORCE;
+  mask |= POST_FORCE_RESPA;
+  mask |= END_OF_STEP;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixTTM::init()
+{
+  if (domain->dimension == 2)
+    error->all(FLERR,"Cannot use fix ttm with 2d simulation");
+  if (domain->nonperiodic != 0)
+    error->all(FLERR,"Cannot use nonperiodic boundares with fix ttm");
+  if (domain->triclinic)
+    error->all(FLERR,"Cannot use fix ttm with triclinic box");
+
+  // set force prefactors
+
+  for (int i = 1; i <= atom->ntypes; i++) {
+    gfactor1[i] = - gamma_p / force->ftm2v;
+    gfactor2[i] =
+      sqrt(24.0*force->boltz*gamma_p/update->dt/force->mvv2e) / force->ftm2v;
+  }
+
+  for (int ixnode = 0; ixnode < nxnodes; ixnode++)
+    for (int iynode = 0; iynode < nynodes; iynode++)
+      for (int iznode = 0; iznode < nznodes; iznode++)
+        net_energy_transfer_all[ixnode][iynode][iznode] = 0;
+
+  if (strstr(update->integrate_style,"respa"))
+    nlevels_respa = ((Respa *) update->integrate)->nlevels;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixTTM::setup(int vflag)
+{
+  if (strstr(update->integrate_style,"verlet"))
+    post_force_setup(vflag);
+  else {
+    ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
+    post_force_respa_setup(vflag,nlevels_respa-1,0);
+    ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixTTM::post_force(int vflag)
+{
+  double **x = atom->x;
+  double **v = atom->v;
+  double **f = atom->f;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  double gamma1,gamma2;
+
+  // apply damping and thermostat to all atoms in fix group
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+
+      double xscale = (x[i][0] - domain->boxlo[0])/domain->xprd;
+      double yscale = (x[i][1] - domain->boxlo[1])/domain->yprd;
+      double zscale = (x[i][2] - domain->boxlo[2])/domain->zprd;
+      int ixnode = static_cast<int>(xscale*nxnodes);
+      int iynode = static_cast<int>(yscale*nynodes);
+      int iznode = static_cast<int>(zscale*nznodes);
+      while (ixnode > nxnodes-1) ixnode -= nxnodes;
+      while (iynode > nynodes-1) iynode -= nynodes;
+      while (iznode > nznodes-1) iznode -= nznodes;
+      while (ixnode < 0) ixnode += nxnodes;
+      while (iynode < 0) iynode += nynodes;
+      while (iznode < 0) iznode += nznodes;
+
+      if (T_electron[ixnode][iynode][iznode] < 0)
+        error->all(FLERR,"Electronic temperature dropped below zero");
+
+      double tsqrt = sqrt(T_electron[ixnode][iynode][iznode]);
+
+      gamma1 = gfactor1[type[i]];
+      double vsq = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
+      if (vsq > v_0_sq) gamma1 *= (gamma_p + gamma_s)/gamma_p;
+      gamma2 = gfactor2[type[i]] * tsqrt;
+
+      flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
+      flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
+      flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
+
+      f[i][0] += flangevin[i][0];
+      f[i][1] += flangevin[i][1];
+      f[i][2] += flangevin[i][2];
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixTTM::post_force_setup(int vflag)
+{
+  double **f = atom->f;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  // apply langevin forces that have been stored from previous run
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      f[i][0] += flangevin[i][0];
+      f[i][1] += flangevin[i][1];
+      f[i][2] += flangevin[i][2];
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixTTM::post_force_respa(int vflag, int ilevel, int iloop)
+{
+  if (ilevel == nlevels_respa-1) post_force(vflag);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixTTM::post_force_respa_setup(int vflag, int ilevel, int iloop)
+{
+  if (ilevel == nlevels_respa-1) post_force_setup(vflag);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixTTM::reset_dt()
+{
+  for (int i = 1; i <= atom->ntypes; i++)
+    gfactor2[i] =
+      sqrt(24.0*force->boltz*gamma_p/update->dt/force->mvv2e) / force->ftm2v;
+}
+
+/* ----------------------------------------------------------------------
+   read in initial electron temperatures from a user-specified file
+   only called by proc 0
+------------------------------------------------------------------------- */
+
+void FixTTM::read_initial_electron_temperatures()
+{
+  char line[MAXLINE];
+
+  for (int ixnode = 0; ixnode < nxnodes; ixnode++)
+    for (int iynode = 0; iynode < nynodes; iynode++)
+      for (int iznode = 0; iznode < nznodes; iznode++)
+        T_initial_set[ixnode][iynode][iznode] = 0;
+
+  // read initial electron temperature values from file
+
+  int ixnode,iynode,iznode;
+  double T_tmp;
+  while (1) {
+    if (fgets(line,MAXLINE,fpr) == NULL) break;
+    sscanf(line,"%d %d %d %lg",&ixnode,&iynode,&iznode,&T_tmp);
+    if (T_tmp < 0.0) error->one(FLERR,"Fix ttm electron temperatures must be > 0.0");
+    T_electron[ixnode][iynode][iznode] = T_tmp;
+    T_initial_set[ixnode][iynode][iznode] = 1;
+  }
+
+  for (int ixnode = 0; ixnode < nxnodes; ixnode++)
+    for (int iynode = 0; iynode < nynodes; iynode++)
+      for (int iznode = 0; iznode < nznodes; iznode++)
+        if (T_initial_set[ixnode][iynode][iznode] == 0)
+          error->one(FLERR,"Initial temperatures not all set in fix ttm");
+
+  // close file
+
+  fclose(fpr);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixTTM::end_of_step()
+{
+  double **x = atom->x;
+  double **v = atom->v;
+  double *mass = atom->mass;
+  double *rmass = atom->rmass;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int ixnode = 0; ixnode < nxnodes; ixnode++)
+    for (int iynode = 0; iynode < nynodes; iynode++)
+      for (int iznode = 0; iznode < nznodes; iznode++)
+        net_energy_transfer[ixnode][iynode][iznode] = 0;
+
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      double xscale = (x[i][0] - domain->boxlo[0])/domain->xprd;
+      double yscale = (x[i][1] - domain->boxlo[1])/domain->yprd;
+      double zscale = (x[i][2] - domain->boxlo[2])/domain->zprd;
+      int ixnode = static_cast<int>(xscale*nxnodes);
+      int iynode = static_cast<int>(yscale*nynodes);
+      int iznode = static_cast<int>(zscale*nznodes);
+      while (ixnode > nxnodes-1) ixnode -= nxnodes;
+      while (iynode > nynodes-1) iynode -= nynodes;
+      while (iznode > nznodes-1) iznode -= nznodes;
+      while (ixnode < 0) ixnode += nxnodes;
+      while (iynode < 0) iynode += nynodes;
+      while (iznode < 0) iznode += nznodes;
+      net_energy_transfer[ixnode][iynode][iznode] +=
+        (flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] +
+         flangevin[i][2]*v[i][2]);
+    }
+
+  MPI_Allreduce(&net_energy_transfer[0][0][0],
+                &net_energy_transfer_all[0][0][0],
+                total_nnodes,MPI_DOUBLE,MPI_SUM,world);
+
+  double dx = domain->xprd/nxnodes;
+  double dy = domain->yprd/nynodes;
+  double dz = domain->zprd/nznodes;
+  double del_vol = dx*dy*dz;
+
+  // num_inner_timesteps = # of inner steps (thermal solves)
+  // required this MD step to maintain a stable explicit solve
+
+  int num_inner_timesteps = 1;
+  double inner_dt = update->dt;
+  double stability_criterion = 1.0 -
+    2.0*inner_dt/(electronic_specific_heat*electronic_density) *
+    (electronic_thermal_conductivity*(1.0/dx/dx + 1.0/dy/dy + 1.0/dz/dz));
+  if (stability_criterion < 0.0) {
+    inner_dt = 0.5*(electronic_specific_heat*electronic_density) /
+      (electronic_thermal_conductivity*(1.0/dx/dx + 1.0/dy/dy + 1.0/dz/dz));
+    num_inner_timesteps = static_cast<int>(update->dt/inner_dt) + 1;
+    inner_dt = update->dt/double(num_inner_timesteps);
+    if (num_inner_timesteps > 1000000)
+      error->warning(FLERR,"Too many inner timesteps in fix ttm",0);
+  }
+
+  for (int ith_inner_timestep = 0; ith_inner_timestep < num_inner_timesteps;
+       ith_inner_timestep++) {
+
+    for (int ixnode = 0; ixnode < nxnodes; ixnode++)
+      for (int iynode = 0; iynode < nynodes; iynode++)
+        for (int iznode = 0; iznode < nznodes; iznode++)
+          T_electron_old[ixnode][iynode][iznode] =
+            T_electron[ixnode][iynode][iznode];
+
+    // compute new electron T profile
+
+    for (int ixnode = 0; ixnode < nxnodes; ixnode++)
+      for (int iynode = 0; iynode < nynodes; iynode++)
+        for (int iznode = 0; iznode < nznodes; iznode++) {
+          int right_xnode = ixnode + 1;
+          int right_ynode = iynode + 1;
+          int right_znode = iznode + 1;
+          if (right_xnode == nxnodes) right_xnode = 0;
+          if (right_ynode == nynodes) right_ynode = 0;
+          if (right_znode == nznodes) right_znode = 0;
+          int left_xnode = ixnode - 1;
+          int left_ynode = iynode - 1;
+          int left_znode = iznode - 1;
+          if (left_xnode == -1) left_xnode = nxnodes - 1;
+          if (left_ynode == -1) left_ynode = nynodes - 1;
+          if (left_znode == -1) left_znode = nznodes - 1;
+          T_electron[ixnode][iynode][iznode] =
+            T_electron_old[ixnode][iynode][iznode] +
+            inner_dt/(electronic_specific_heat*electronic_density) *
+            (electronic_thermal_conductivity *
+             ((T_electron_old[right_xnode][iynode][iznode] +
+               T_electron_old[left_xnode][iynode][iznode] -
+               2*T_electron_old[ixnode][iynode][iznode])/dx/dx +
+              (T_electron_old[ixnode][right_ynode][iznode] +
+               T_electron_old[ixnode][left_ynode][iznode] -
+               2*T_electron_old[ixnode][iynode][iznode])/dy/dy +
+              (T_electron_old[ixnode][iynode][right_znode] +
+               T_electron_old[ixnode][iynode][left_znode] -
+               2*T_electron_old[ixnode][iynode][iznode])/dz/dz) -
+              (net_energy_transfer_all[ixnode][iynode][iznode])/del_vol);
+        }
+  }
+
+  // output nodal temperatures for current timestep
+
+  if ((nfileevery) && !(update->ntimestep % nfileevery)) {
+
+    // compute atomic Ta for each grid point
+
+    for (int ixnode = 0; ixnode < nxnodes; ixnode++)
+      for (int iynode = 0; iynode < nynodes; iynode++)
+        for (int iznode = 0; iznode < nznodes; iznode++) {
+          nsum[ixnode][iynode][iznode] = 0;
+          nsum_all[ixnode][iynode][iznode] = 0;
+          sum_vsq[ixnode][iynode][iznode] = 0.0;
+          sum_mass_vsq[ixnode][iynode][iznode] = 0.0;
+          sum_vsq_all[ixnode][iynode][iznode] = 0.0;
+          sum_mass_vsq_all[ixnode][iynode][iznode] = 0.0;
+        }
+
+    double massone;
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) {
+        if (rmass) massone = rmass[i];
+        else massone = mass[type[i]];
+        double xscale = (x[i][0] - domain->boxlo[0])/domain->xprd;
+        double yscale = (x[i][1] - domain->boxlo[1])/domain->yprd;
+        double zscale = (x[i][2] - domain->boxlo[2])/domain->zprd;
+        int ixnode = static_cast<int>(xscale*nxnodes);
+        int iynode = static_cast<int>(yscale*nynodes);
+        int iznode = static_cast<int>(zscale*nznodes);
+        while (ixnode > nxnodes-1) ixnode -= nxnodes;
+        while (iynode > nynodes-1) iynode -= nynodes;
+        while (iznode > nznodes-1) iznode -= nznodes;
+        while (ixnode < 0) ixnode += nxnodes;
+        while (iynode < 0) iynode += nynodes;
+        while (iznode < 0) iznode += nznodes;
+        double vsq = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
+        nsum[ixnode][iynode][iznode] += 1;
+        sum_vsq[ixnode][iynode][iznode] += vsq;
+        sum_mass_vsq[ixnode][iynode][iznode] += massone*vsq;
+      }
+
+    MPI_Allreduce(&nsum[0][0][0],&nsum_all[0][0][0],total_nnodes,
+                  MPI_INT,MPI_SUM,world);
+    MPI_Allreduce(&sum_vsq[0][0][0],&sum_vsq_all[0][0][0],total_nnodes,
+                  MPI_DOUBLE,MPI_SUM,world);
+    MPI_Allreduce(&sum_mass_vsq[0][0][0],&sum_mass_vsq_all[0][0][0],
+                  total_nnodes,MPI_DOUBLE,MPI_SUM,world);
+
+    if (me == 0) {
+      fprintf(fp,BIGINT_FORMAT,update->ntimestep);
+
+      double T_a;
+      for (int ixnode = 0; ixnode < nxnodes; ixnode++)
+        for (int iynode = 0; iynode < nynodes; iynode++)
+          for (int iznode = 0; iznode < nznodes; iznode++) {
+            T_a = 0;
+            if (nsum_all[ixnode][iynode][iznode] > 0)
+              T_a = sum_mass_vsq_all[ixnode][iynode][iznode]/
+                (3.0*force->boltz*nsum_all[ixnode][iynode][iznode]/force->mvv2e);
+            fprintf(fp," %f",T_a);
+          }
+
+      fprintf(fp,"\t");
+      for (int ixnode = 0; ixnode < nxnodes; ixnode++)
+        for (int iynode = 0; iynode < nynodes; iynode++)
+          for (int iznode = 0; iznode < nznodes; iznode++)
+            fprintf(fp,"%f ",T_electron[ixnode][iynode][iznode]);
+      fprintf(fp,"\n");
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of 3d grid
+------------------------------------------------------------------------- */
+
+double FixTTM::memory_usage()
+{
+  double bytes = 0.0;
+  bytes += 5*total_nnodes * sizeof(int);
+  bytes += 14*total_nnodes * sizeof(double);
+  return bytes;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixTTM::grow_arrays(int ngrow)
+{
+
+ memory->grow(flangevin,ngrow,3,"TTM:flangevin");
+
+}
+
+/* ----------------------------------------------------------------------
+  return the energy of the electronic subsystem or the net_energy transfer
+   between the subsystems
+------------------------------------------------------------------------- */
+
+double FixTTM::compute_vector(int n)
+{
+  double e_energy = 0.0;
+  double transfer_energy = 0.0;
+
+  double dx = domain->xprd/nxnodes;
+  double dy = domain->yprd/nynodes;
+  double dz = domain->zprd/nznodes;
+  double del_vol = dx*dy*dz;
+
+  for (int ixnode = 0; ixnode < nxnodes; ixnode++)
+    for (int iynode = 0; iynode < nynodes; iynode++)
+      for (int iznode = 0; iznode < nznodes; iznode++) {
+        e_energy +=
+          T_electron[ixnode][iynode][iznode]*electronic_specific_heat*
+          electronic_density*del_vol;
+        transfer_energy +=
+          net_energy_transfer_all[ixnode][iynode][iznode]*update->dt;
+  }
+
+  if (n == 0) return e_energy;
+  if (n == 1) return transfer_energy;
+  return 0.0;
+}
+
+/* ----------------------------------------------------------------------
+   pack entire state of Fix into one write
+------------------------------------------------------------------------- */
+
+void FixTTM::write_restart(FILE *fp)
+{
+  double *rlist;
+  memory->create(rlist,nxnodes*nynodes*nznodes+1,"TTM:rlist");
+
+  int n = 0;
+  rlist[n++] = seed;
+
+  for (int ixnode = 0; ixnode < nxnodes; ixnode++)
+    for (int iynode = 0; iynode < nynodes; iynode++)
+      for (int iznode = 0; iznode < nznodes; iznode++)
+        rlist[n++] =  T_electron[ixnode][iynode][iznode];
+
+  if (comm->me == 0) {
+    int size = n * sizeof(double);
+    fwrite(&size,sizeof(int),1,fp);
+    fwrite(rlist,sizeof(double),n,fp);
+  }
+
+  memory->destroy(rlist);
+}
+
+/* ----------------------------------------------------------------------
+   use state info from restart file to restart the Fix
+------------------------------------------------------------------------- */
+
+void FixTTM::restart(char *buf)
+{
+  int n = 0;
+  double *rlist = (double *) buf;
+
+  // the seed must be changed from the initial seed
+
+  seed = static_cast<int> (0.5*rlist[n++]);
+
+  for (int ixnode = 0; ixnode < nxnodes; ixnode++)
+    for (int iynode = 0; iynode < nynodes; iynode++)
+      for (int iznode = 0; iznode < nznodes; iznode++)
+        T_electron[ixnode][iynode][iznode] = rlist[n++];
+
+  delete random;
+  random = new RanMars(lmp,seed+comm->me);
+}
+
+/* ----------------------------------------------------------------------
+   pack values in local atom-based arrays for restart file
+------------------------------------------------------------------------- */
+
+int FixTTM::pack_restart(int i, double *buf)
+{
+  buf[0] = 4;
+  buf[1] = flangevin[i][0];
+  buf[2] = flangevin[i][1];
+  buf[3] = flangevin[i][2];
+  return 4;
+}
+
+/* ----------------------------------------------------------------------
+   unpack values from atom->extra array to restart the fix
+------------------------------------------------------------------------- */
+
+void FixTTM::unpack_restart(int nlocal, int nth)
+{
+  double **extra = atom->extra;
+
+  // skip to Nth set of extra values
+
+  int m = 0;
+  for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);
+  m++;
+
+  flangevin[nlocal][0] = extra[nlocal][m++];
+  flangevin[nlocal][1] = extra[nlocal][m++];
+  flangevin[nlocal][2] = extra[nlocal][m++];
+}
+
+/* ----------------------------------------------------------------------
+   maxsize of any atom's restart data
+------------------------------------------------------------------------- */
+
+int FixTTM::maxsize_restart()
+{
+  return 4;
+}
+
+/* ----------------------------------------------------------------------
+   size of atom nlocal's restart data
+------------------------------------------------------------------------- */
+
+int FixTTM::size_restart(int nlocal)
+{
+  return 4;
+}

src/MISC/fix_ttm.h

@@ -0,0 +1,156 @@
+/* ----------------------------------------------------------------------
+   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(ttm,FixTTM)
+
+#else
+
+#ifndef LMP_FIX_TTM_H
+#define LMP_FIX_TTM_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixTTM : public Fix {
+ public:
+  FixTTM(class LAMMPS *, int, char **);
+  ~FixTTM();
+  int setmask();
+  void init();
+  void setup(int);
+  void post_force(int);
+  void post_force_respa(int, int, int);
+  void post_force_setup(int);
+  void post_force_respa_setup(int, int, int);
+  void end_of_step();
+  void reset_dt();
+  void write_restart(FILE *);
+  void restart(char *);
+  int pack_restart(int, double *);
+  void unpack_restart(int, int);
+  int size_restart(int);
+  int maxsize_restart();
+  double memory_usage();
+  void grow_arrays(int);
+  double compute_vector(int);
+
+ private:
+  int me;
+  int nfileevery;
+  int nlevels_respa;
+  int seed;
+  class RanMars *random;
+  FILE *fp,*fpr;
+  int nxnodes,nynodes,nznodes,total_nnodes;
+  int ***nsum;
+  int ***nsum_all,***T_initial_set;
+  double *gfactor1,*gfactor2,*ratio;
+  double **flangevin;
+  double ***T_electron,***T_electron_old;
+  double ***sum_vsq,***sum_mass_vsq;
+  double ***sum_vsq_all,***sum_mass_vsq_all;
+  double ***net_energy_transfer,***net_energy_transfer_all;
+  double electronic_specific_heat,electronic_density;
+  double electronic_thermal_conductivity;
+  double gamma_p,gamma_s,v_0,v_0_sq;
+
+  void read_initial_electron_temperatures();
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Cannot open file %s
+
+The specified file cannot be opened.  Check that the path and name are
+correct.
+
+E: Cannot open fix ttm file %s
+
+The output file for the fix ttm command cannot be opened.  Check that
+the path and name are correct.
+
+E: Invalid random number seed in fix ttm command
+
+Random number seed must be > 0.
+
+E: Fix ttm electronic_specific_heat must be > 0.0
+
+Self-explanatory.
+
+E: Fix ttm electronic_density must be > 0.0
+
+Self-explanatory.
+
+E: Fix ttm electronic_thermal_conductivity must be >= 0.0
+
+Self-explanatory.
+
+E: Fix ttm gamma_p must be > 0.0
+
+Self-explanatory.
+
+E: Fix ttm gamma_s must be >= 0.0
+
+Self-explanatory.
+
+E: Fix ttm v_0 must be >= 0.0
+
+Self-explanatory.
+
+E: Fix ttm number of nodes must be > 0
+
+Self-explanatory.
+
+E: Cannot use fix ttm with 2d simulation
+
+This is a current restriction of this fix due to the grid it creates.
+
+E: Cannot use nonperiodic boundares with fix ttm
+
+This fix requires a fully periodic simulation box.
+
+E: Cannot use fix ttm with triclinic box
+
+This is a current restriction of this fix due to the grid it creates.
+
+E: Electronic temperature dropped below zero
+
+Something has gone wrong with the fix ttm electron temperature model.
+
+E: Fix ttm electron temperatures must be > 0.0
+
+Self-explanatory.
+
+E: Initial temperatures not all set in fix ttm
+
+Self-explantory.
+
+W: Too many inner timesteps in fix ttm
+
+Self-explanatory.
+
+*/

src/MISC/fix_viscosity.cpp

@@ -0,0 +1,309 @@
+/* ----------------------------------------------------------------------
+   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: Craig Tenney (UND) added support
+                        for swapping atoms of different masses
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "mpi.h"
+#include "string.h"
+#include "stdlib.h"
+#include "fix_viscosity.h"
+#include "atom.h"
+#include "domain.h"
+#include "modify.h"
+#include "error.h"
+#include "force.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+// needs to be big, but not so big that lose precision when subtract velocity
+
+#define BIG 1.0e10
+
+/* ---------------------------------------------------------------------- */
+
+FixViscosity::FixViscosity(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg)
+{
+  if (narg < 7) error->all(FLERR,"Illegal fix viscosity command");
+
+  MPI_Comm_rank(world,&me);
+
+  nevery = force->inumeric(FLERR,arg[3]);
+  if (nevery <= 0) error->all(FLERR,"Illegal fix viscosity command");
+
+  scalar_flag = 1;
+  global_freq = nevery;
+  extscalar = 0;
+
+  if (strcmp(arg[4],"x") == 0) vdim = 0;
+  else if (strcmp(arg[4],"y") == 0) vdim = 1;
+  else if (strcmp(arg[4],"z") == 0) vdim = 2;
+  else error->all(FLERR,"Illegal fix viscosity command");
+
+  if (strcmp(arg[5],"x") == 0) pdim = 0;
+  else if (strcmp(arg[5],"y") == 0) pdim = 1;
+  else if (strcmp(arg[5],"z") == 0) pdim = 2;
+  else error->all(FLERR,"Illegal fix viscosity command");
+
+  nbin = force->inumeric(FLERR,arg[6]);
+  if (nbin % 2 || nbin <= 2) error->all(FLERR,"Illegal fix viscosity command");
+
+  // optional keywords
+
+  nswap = 1;
+  vtarget = BIG;
+
+  int iarg = 7;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"swap") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix viscosity command");
+      nswap = force->inumeric(FLERR,arg[iarg+1]);
+      if (nswap <= 0) 
+        error->all(FLERR,"Fix viscosity swap value must be positive");
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"vtarget") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix viscosity command");
+      if (strcmp(arg[iarg+1],"INF") == 0) vtarget = BIG;
+      else vtarget = force->numeric(FLERR,arg[iarg+1]);
+      if (vtarget <= 0.0)
+        error->all(FLERR,"Fix viscosity vtarget value must be positive");
+      iarg += 2;
+    } else error->all(FLERR,"Illegal fix viscosity command");
+  }
+
+  // initialize array sizes to nswap+1 so have space to shift values down
+
+  pos_index = new int[nswap+1];
+  neg_index = new int[nswap+1];
+  pos_delta = new double[nswap+1];
+  neg_delta = new double[nswap+1];
+
+  p_exchange = 0.0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixViscosity::~FixViscosity()
+{
+  delete [] pos_index;
+  delete [] neg_index;
+  delete [] pos_delta;
+  delete [] neg_delta;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixViscosity::setmask()
+{
+  int mask = 0;
+  mask |= END_OF_STEP;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixViscosity::init()
+{
+  // warn if any fix ave/spatial comes after this fix
+  // can cause glitch in averaging since ave will happen after swap
+
+  int foundme = 0;
+  for (int i = 0; i < modify->nfix; i++) {
+    if (modify->fix[i] == this) foundme = 1;
+    if (foundme && strcmp(modify->fix[i]->style,"ave/spatial") == 0 && me == 0)
+      error->warning(FLERR,"Fix viscosity comes before fix ave/spatial");
+  }
+
+  // set bounds of 2 slabs in pdim
+  // only necessary for static box, else re-computed in end_of_step()
+  // lo bin is always bottom bin
+  // hi bin is just above half height
+
+  if (domain->box_change == 0) {
+    prd = domain->prd[pdim];
+    boxlo = domain->boxlo[pdim];
+    boxhi = domain->boxhi[pdim];
+    double binsize = (boxhi-boxlo) / nbin;
+    slablo_lo = boxlo;
+    slablo_hi = boxlo + binsize;
+    slabhi_lo = boxlo + (nbin/2)*binsize;
+    slabhi_hi = boxlo + (nbin/2+1)*binsize;
+  }
+
+  periodicity = domain->periodicity[pdim];
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixViscosity::end_of_step()
+{
+  int i,m,insert;
+  double coord,delta;
+  MPI_Status status;
+  struct {
+    double value;
+    int proc;
+  } mine[2],all[2];
+
+  // if box changes, recompute bounds of 2 slabs in pdim
+
+  if (domain->box_change) {
+    prd = domain->prd[pdim];
+    boxlo = domain->boxlo[pdim];
+    boxhi = domain->boxhi[pdim];
+    double binsize = (boxhi-boxlo) / nbin;
+    slablo_lo = boxlo;
+    slablo_hi = boxlo + binsize;
+    slabhi_lo = boxlo + (nbin/2)*binsize;
+    slabhi_hi = boxlo + (nbin/2+1)*binsize;
+  }
+
+  // make 2 lists of up to nswap atoms with velocity closest to +/- vtarget
+  // lists are sorted by closeness to vtarget
+  // only consider atoms in the bottom/middle slabs
+  // map atoms back into periodic box if necessary
+  // insert = location in list to insert new atom
+
+  double **x = atom->x;
+  double **v = atom->v;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  npositive = nnegative = 0;
+
+  for (i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      coord = x[i][pdim];
+      if (coord < boxlo && periodicity) coord += prd;
+      else if (coord >= boxhi && periodicity) coord -= prd;
+
+      if (coord >= slablo_lo && coord < slablo_hi) {
+        if (v[i][vdim] < 0.0) continue;
+        delta = fabs(v[i][vdim] - vtarget);
+        if (npositive < nswap || delta < pos_delta[nswap-1]) {
+          for (insert = npositive-1; insert >= 0; insert--)
+            if (delta > pos_delta[insert]) break;
+          insert++;
+          for (m = npositive-1; m >= insert; m--) {
+            pos_delta[m+1] = pos_delta[m];
+            pos_index[m+1] = pos_index[m];
+          }
+          pos_delta[insert] = delta;
+          pos_index[insert] = i;
+          if (npositive < nswap) npositive++;
+        }
+      }
+
+      if (coord >= slabhi_lo && coord < slabhi_hi) {
+        if (v[i][vdim] > 0.0) continue;
+        delta = fabs(v[i][vdim] + vtarget);
+        if (nnegative < nswap || delta < neg_delta[nswap-1]) {
+          for (insert = nnegative-1; insert >= 0; insert--)
+            if (delta > neg_delta[insert]) break;
+          insert++;
+          for (m = nnegative-1; m >= insert; m--) {
+            neg_delta[m+1] = neg_delta[m];
+            neg_index[m+1] = neg_index[m];
+          }
+          neg_delta[insert] = delta;
+          neg_index[insert] = i;
+          if (nnegative < nswap) nnegative++;
+        }
+      }
+    }
+
+  // loop over nswap pairs
+  // find 2 global atoms with smallest delta in bottom/top slabs
+  // BIG values are for procs with no atom to contribute
+  // MINLOC also communicates which procs own them
+  // exchange momenta between the 2 particles
+  // if I own both particles just swap, else point2point comm of vel,mass
+
+  double *mass = atom->mass;
+  double *rmass = atom->rmass;
+
+  int ipos,ineg;
+  double sbuf[2],rbuf[2],vcm;
+
+  double pswap = 0.0;
+  mine[0].proc = mine[1].proc = me;
+  int ipositive = 0;
+  int inegative = 0;
+
+  for (m = 0; m < nswap; m++) {
+    if (ipositive < npositive) mine[0].value = pos_delta[ipositive];
+    else mine[0].value = BIG;
+    if (inegative < nnegative) mine[1].value = neg_delta[inegative];
+    else mine[1].value = BIG;
+
+    MPI_Allreduce(mine,all,2,MPI_DOUBLE_INT,MPI_MINLOC,world);
+
+    if (all[0].value == BIG || all[1].value == BIG) continue;
+
+    if (me == all[0].proc && me == all[1].proc) {
+      ipos = pos_index[ipositive++];
+      ineg = neg_index[inegative++];
+      rbuf[0] = v[ipos][vdim];
+      if (rmass) rbuf[1] = rmass[ipos];
+      else rbuf[1] = mass[type[ipos]];
+      sbuf[0] = v[ineg][vdim];
+      if (rmass) sbuf[1] = rmass[ineg];
+      else sbuf[1] = mass[type[ineg]];
+      vcm = (sbuf[1]*sbuf[0] + rbuf[1]*rbuf[0]) / (sbuf[1] + rbuf[1]);
+      v[ineg][vdim] = 2.0 * vcm - sbuf[0];
+      v[ipos][vdim] = 2.0 * vcm - rbuf[0];
+      pswap += rbuf[1] * (vcm - rbuf[0]) - sbuf[1] * (vcm - sbuf[0]);
+
+    } else if (me == all[0].proc) {
+      ipos = pos_index[ipositive++];
+      sbuf[0] = v[ipos][vdim];
+      if (rmass) sbuf[1] = rmass[ipos];
+      else sbuf[1] = mass[type[ipos]];
+      MPI_Sendrecv(sbuf,2,MPI_DOUBLE,all[1].proc,0,
+                   rbuf,2,MPI_DOUBLE,all[1].proc,0,world,&status);
+      vcm = (sbuf[1]*sbuf[0] + rbuf[1]*rbuf[0]) / (sbuf[1] + rbuf[1]);
+      v[ipos][vdim] = 2.0 * vcm - sbuf[0];
+      pswap += sbuf[1] * (vcm - sbuf[0]);
+
+    } else if (me == all[1].proc) {
+      ineg = neg_index[inegative++];
+      sbuf[0] = v[ineg][vdim];
+      if (rmass) sbuf[1] = rmass[ineg];
+      else sbuf[1] = mass[type[ineg]];
+      MPI_Sendrecv(sbuf,2,MPI_DOUBLE,all[0].proc,0,
+                   rbuf,2,MPI_DOUBLE,all[0].proc,0,world,&status);
+      vcm = (sbuf[1]*sbuf[0] + rbuf[1]*rbuf[0]) / (sbuf[1] + rbuf[1]);
+      v[ineg][vdim] = 2.0 * vcm - sbuf[0];
+      pswap -= sbuf[1] * (vcm - sbuf[0]);
+    }
+  }
+
+  // tally momentum exchange from all swaps
+
+  double pswap_all;
+  MPI_Allreduce(&pswap,&pswap_all,1,MPI_DOUBLE,MPI_SUM,world);
+  p_exchange += pswap_all;
+}
+
+/* ---------------------------------------------------------------------- */
+
+double FixViscosity::compute_scalar()
+{
+  return p_exchange;
+}

src/MISC/fix_viscosity.h

@@ -0,0 +1,80 @@
+/* ----------------------------------------------------------------------
+   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(viscosity,FixViscosity)
+
+#else
+
+#ifndef LMP_FIX_VISCOSITY_H
+#define LMP_FIX_VISCOSITY_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixViscosity : public Fix {
+ public:
+  FixViscosity(class LAMMPS *, int, char **);
+  ~FixViscosity();
+  int setmask();
+  void init();
+  void end_of_step();
+  double compute_scalar();
+
+ private:
+  int me;
+  int vdim,pdim,nbin,periodicity;
+  int nswap;
+  double vtarget;
+  double prd,boxlo,boxhi;
+  double slablo_lo,slablo_hi,slabhi_lo,slabhi_hi;
+  double p_exchange;
+
+  int npositive,nnegative;
+  int *pos_index,*neg_index;
+  double *pos_delta,*neg_delta;
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Fix viscosity swap value must be positive
+
+Self-explanatory.
+
+E: Fix viscosity vtarget value must be positive
+
+Self-explanatory.
+
+W: Fix viscosity comes before fix ave/spatial
+
+The order of these 2 fixes in your input script is such that
+fix viscosity comes first.  If you are using fix ave/spatial
+to measure the velocity profile induced by fix viscosity, then
+this may cause a glitch in the profile since you are averaging
+immediately after swaps have occurred.  Flipping the order
+of the 2 fixes typically helps.
+
+*/

src/compute_msd_nongauss.cpp

@@ -1,106 +0,0 @@
-/* ----------------------------------------------------------------------
-   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
-   http://lammps.sandia.gov, Sandia National Laboratories
-   Steve Plimpton, sjplimp@sandia.gov
-
-   Copyright (2003) Sandia Corporation.  Under the terms of Contract
-   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
-   certain rights in this software.  This software is distributed under
-   the GNU General Public License.
-
-   See the README file in the top-level LAMMPS directory.
-------------------------------------------------------------------------- */
-
-/* ----------------------------------------------------------------------
-   Contributing authors: Rob Hoy
-------------------------------------------------------------------------- */
-
-#include "string.h"
-#include "compute_msd_nongauss.h"
-#include "atom.h"
-#include "update.h"
-#include "group.h"
-#include "domain.h"
-#include "fix_store.h"
-#include "error.h"
-
-using namespace LAMMPS_NS;
-
-/* ---------------------------------------------------------------------- */
-
-ComputeMSDNonGauss::ComputeMSDNonGauss(LAMMPS *lmp, int narg, char **arg) :
-  ComputeMSD(lmp, narg, arg)
-{
-  size_vector = 3;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void ComputeMSDNonGauss::compute_vector()
-{
-  invoked_vector = update->ntimestep;
-
-  // cm = current center of mass
-
-  double cm[3];
-  if (comflag) group->xcm(igroup,masstotal,cm);
-  else cm[0] = cm[1] = cm[2] = 0.0;
-
-  // dx,dy,dz = displacement of atom from original position
-  // original unwrapped position is stored by fix
-  // relative to center of mass if comflag is set
-  // for triclinic, need to unwrap current atom coord via h matrix
-
-  double **xoriginal = fix->astore;
-
-  double **x = atom->x;
-  int *mask = atom->mask;
-  tagint *image = atom->image;
-  int nlocal = atom->nlocal;
-
-  double *h = domain->h;
-  double xprd = domain->xprd;
-  double yprd = domain->yprd;
-  double zprd = domain->zprd;
-
-  double dx,dy,dz;
-  int xbox,ybox,zbox;
-
-  double msd[2];
-  msd[0] = msd[1] = 0.0;
-
-  if (domain->triclinic == 0) {
-    for (int i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) {
-        xbox = (image[i] & IMGMASK) - IMGMAX;
-        ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
-        zbox = (image[i] >> IMG2BITS) - IMGMAX;
-        dx = x[i][0] + xbox*xprd - cm[0] - xoriginal[i][0];
-        dy = x[i][1] + ybox*yprd - cm[1] - xoriginal[i][1];
-        dz = x[i][2] + zbox*zprd - cm[2] - xoriginal[i][2];
-        msd[0] += dx*dx + dy*dy + dz*dz;
-        msd[1] += (dx*dx + dy*dy + dz*dz)*(dx*dx + dy*dy + dz*dz);
-      }
-
-  } else {
-    for (int i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) {
-        xbox = (image[i] & IMGMASK) - IMGMAX;
-        ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
-        zbox = (image[i] >> IMG2BITS) - IMGMAX;
-        dx = x[i][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox -
-          cm[0] - xoriginal[i][0];
-        dy = x[i][1] + h[1]*ybox + h[3]*zbox - cm[1] - xoriginal[i][1];
-        dz = x[i][2] + h[2]*zbox - cm[2] - xoriginal[i][2];
-        msd[0] += dx*dx + dy*dy + dz*dz;
-        msd[1] += (dx*dx + dy*dy + dz*dz)*(dx*dx + dy*dy + dz*dz);
-      }
-  }
-
-  MPI_Allreduce(msd,vector,2,MPI_DOUBLE,MPI_SUM,world);
-  if (nmsd) {
-    vector[0] /= nmsd;
-    vector[1] /= nmsd;
-    vector[2] = (3*vector[1])/(5*vector[0]*vector[0]) - 1;
-  }
-}

src/compute_msd_nongauss.h

@@ -1,51 +0,0 @@
-/* ----------------------------------------------------------------------
-   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
-   http://lammps.sandia.gov, Sandia National Laboratories
-   Steve Plimpton, sjplimp@sandia.gov
-
-   Copyright (2003) Sandia Corporation.  Under the terms of Contract
-   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
-   certain rights in this software.  This software is distributed under
-   the GNU General Public License.
-
-   See the README file in the top-level LAMMPS directory.
-------------------------------------------------------------------------- */
-
-#ifdef COMPUTE_CLASS
-
-ComputeStyle(msd/nongauss,ComputeMSDNonGauss)
-
-#else
-
-#ifndef LMP_COMPUTE_MSD_NONGAUSS_H
-#define LMP_COMPUTE_MSD_NONGAUSS_H
-
-#include "compute_msd.h"
-
-namespace LAMMPS_NS {
-
-class ComputeMSDNonGauss : public ComputeMSD {
- public:
-  ComputeMSDNonGauss(class LAMMPS *, int, char **);
-  ~ComputeMSDNonGauss() {}
-  void compute_vector();
-};
-
-}
-
-#endif
-#endif
-
-/* ERROR/WARNING messages:
-
-E: Illegal ... command
-
-Self-explanatory.  Check the input script syntax and compare to the
-documentation for the command.  You can use -echo screen as a
-command-line option when running LAMMPS to see the offending line.
-
-E: Could not find compute msd fix ID
-
-Self-explanatory.
-
-*/

src/compute_ti.cpp

@@ -1,241 +0,0 @@
-/* ----------------------------------------------------------------------
-   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
-   http://lammps.sandia.gov, Sandia National Laboratories
-   Steve Plimpton, sjplimp@sandia.gov
-
-   Copyright (2003) Sandia Corporation.  Under the terms of Contract
-   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
-   certain rights in this software.  This software is distributed under
-   the GNU General Public License.
-
-   See the README file in the top-level LAMMPS directory.
-------------------------------------------------------------------------- */
-
-/* ----------------------------------------------------------------------
-   Contributing author: Sai Jayaraman (University of Notre Dame)
-------------------------------------------------------------------------- */
-
-#include "mpi.h"
-#include "atom.h"
-#include "string.h"
-#include "compute_ti.h"
-#include "update.h"
-#include "modify.h"
-#include "domain.h"
-#include "force.h"
-#include "pair.h"
-#include "kspace.h"
-#include "input.h"
-#include "variable.h"
-#include "error.h"
-
-using namespace LAMMPS_NS;
-
-enum{PAIR,TAIL,KSPACE};
-
-/* ---------------------------------------------------------------------- */
-
-ComputeTI::ComputeTI(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg)
-{
-  if (narg < 4) error->all(FLERR,"Illegal compute ti command");
-
-  peflag = 1;
-  peratom_flag = 1;
-  peatomflag = 1; 
-  scalar_flag = 1;
-  extscalar = 1;
-  timeflag = 1;
-
-  // terms come in triplets
-  // changed to quadruplets to include atom type
-
-  nterms = (narg-3) / 4;
-  if (narg != 4*nterms + 3) error->all(FLERR,"Illegal compute ti command");
-
-  which = new int[nterms];
-  ivar1 = new int[nterms];
-  ivar2 = new int[nterms];
-  ilo = new int[nterms];
-  ihi = new int[nterms];
-  var1 = new char*[nterms];
-  var2 = new char*[nterms];
-  pptr = new Pair*[nterms];
-  pstyle = new char*[nterms];
-
-  for (int m = 0; m < nterms; m++) pstyle[m] = NULL;
-
-  // parse keywords
-
-  nterms = 0;
-
-  int iarg = 3;
-  while (iarg < narg) {
-    if (iarg+4 > narg) error->all(FLERR,"Illegal compute ti command");
-    if (strcmp(arg[iarg],"kspace") == 0) which[nterms] = KSPACE;
-    else if (strcmp(arg[iarg],"tail") == 0) which[nterms] = TAIL;
-    else which[nterms] = PAIR;
-
-    int n = strlen(arg[iarg]) + 1;
-    pstyle[nterms] = new char[n];
-    strcpy(pstyle[nterms],arg[iarg]);
-    force->bounds(arg[iarg+1],atom->ntypes,ilo[nterms],ihi[nterms]);
-    iarg += 1;
-
-    if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
-      int n = strlen(&arg[iarg+1][2]) + 1;
-      var1[nterms] = new char[n];
-      strcpy(var1[nterms],&arg[iarg+1][2]);
-    } else error->all(FLERR,"Illegal compute ti command");
-    if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) {
-      int n = strlen(&arg[iarg+2][2]) + 1;
-      var2[nterms] = new char[n];
-      strcpy(var2[nterms],&arg[iarg+2][2]);
-    } else error->all(FLERR,"Illegal compute ti command");
-
-    nterms++;
-    iarg += 3;
-  }
-}
-
-/* --------------------------------------------------------------------- */
-
-ComputeTI::~ComputeTI()
-{
-  for (int m = 0; m < nterms; m++) {
-    delete [] var1[m];
-    delete [] var2[m];
-    delete [] pstyle[m];
-  }
-  delete [] which;
-  delete [] ivar1;
-  delete [] ivar2;
-  delete [] var1;
-  delete [] var2;
-  delete [] ilo;
-  delete [] ihi;
-  delete [] pptr;
-  delete [] pstyle;
-}
-
-/* --------------------------------------------------------------------- */
-
-void ComputeTI::init()
-{
-  // setup and error checks
-
-  for (int m = 0; m < nterms; m++) {
-    ivar1[m] = input->variable->find(var1[m]);
-    ivar2[m] = input->variable->find(var2[m]);
-    if (ivar1[m] < 0 || ivar2 < 0)
-      error->all(FLERR,"Variable name for compute ti does not exist");
-    if (!input->variable->equalstyle(ivar1[m]) ||
-        !input->variable->equalstyle(ivar2[m]))
-      error->all(FLERR,"Variable for compute ti is invalid style");
-
-    if (which[m] == PAIR) {
-      pptr[m] = force->pair_match(pstyle[m],1);
-      if (pptr[m] == NULL)
-        error->all(FLERR,"Compute ti pair style does not exist");
-
-    } else if (which[m] == TAIL) {
-      if (force->pair == NULL || force->pair->tail_flag == 0)
-        error->all(FLERR,"Compute ti tail when pair style does not "
-                   "compute tail corrections");
-
-    } else if (which[m] == KSPACE) {
-      if (force->kspace == NULL)
-        error->all(FLERR,"Compute ti kspace style does not exist");
-    }
-  }
-}
-
-/* --------------------------------------------------------------------- */
-
-double ComputeTI::compute_scalar()
-{
-  double eng,engall,value1,value2;
-
-  invoked_scalar = update->ntimestep;
-  if (update->eflag_global != invoked_scalar)
-    error->all(FLERR,"Energy was not tallied on needed timestep");
-
-  double dUdl = 0.0;
-
-  for (int m = 0; m < nterms; m++) {
-    int total_flag = 0;
-    if ((ihi[m]-ilo[m])==atom->ntypes) total_flag = 1; 
-    eng = 0.0;
-    value1 = input->variable->compute_equal(ivar1[m]);
-    value2 = input->variable->compute_equal(ivar2[m]);
-    if (value1 == 0.0) continue;
-
-    if (which[m] == PAIR) {
-      int ntypes = atom->ntypes;
-      int *mask = atom->mask; 
-      if (total_flag) {
-        eng = pptr[m]->eng_vdwl + pptr[m]->eng_coul;
-        MPI_Allreduce(&eng,&engall,1,MPI_DOUBLE,MPI_SUM,world);
-      } 
-      else { 
-        int nlocal = atom->nlocal;
-        int npair = nlocal;
-        int *mask = atom->mask;
-        int *type = atom->type;
-        
-        double *eatom = pptr[m]->eatom;
-        
-        if (force->newton) npair += atom->nghost;
-        for (int i = 0; i < npair; i++)    
-          if ((ilo[m]<=type[i])&(ihi[m]>=type[i])) eng += eatom[i];  
-        MPI_Allreduce(&eng,&engall,1,MPI_DOUBLE,MPI_SUM,world);
-      }
-      dUdl += engall/value1 * value2;
-      
-    } else if (which[m] == TAIL) {
-      double vol = domain->xprd*domain->yprd*domain->zprd;
-      if (total_flag) 
-        eng = force->pair->etail / vol;
-      else {
-        eng = 0; 
-        for (int it = 1; it <= atom->ntypes; it++) {
-          int jt; 
-          if ((it >= ilo[m])&&(it <=ihi[m])) jt = it;
-          else jt = ilo[m];
-          for (; jt <=ihi[m];jt++) {
-            if ((force->pair->tail_flag)&&(force->pair->setflag[it][jt])) {
-              double cut = force->pair->init_one(it,jt);
-              eng += force->pair->etail_ij;
-            }
-            if (it !=jt) eng += force->pair->etail_ij; 
-          }
-        }
-        eng /= vol; 
-      }
-      dUdl += eng/value1 * value2;
-      
-    } else if (which[m] == KSPACE) {
-      int ntypes = atom->ntypes;
-      int *mask = atom->mask; 
-      if (total_flag) 
-        eng = force->kspace->energy;
-      else { 
-        int nlocal = atom->nlocal;
-        int npair = nlocal;
-        int *mask = atom->mask;
-        int *type = atom->type;
-        double *eatom = force->kspace->eatom;
-        eng = 0;
-        for(int i = 0; i < nlocal; i++)
-          if ((ilo[m]<=type[i])&(ihi[m]>=type[i])) 
-            eng += eatom[i];  
-        MPI_Allreduce(&eng,&engall,1,MPI_DOUBLE,MPI_SUM,world);
-        eng = engall;
-      }
-      dUdl += eng/value1 * value2;
-    }
-  }
-  
-  scalar = dUdl;
-  return scalar;
-}

src/fix_adapt.cpp

@@ -1,389 +0,0 @@
-/* ----------------------------------------------------------------------
-   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
-   http://lammps.sandia.gov, Sandia National Laboratories
-   Steve Plimpton, sjplimp@sandia.gov
-
-   Copyright (2003) Sandia Corporation.  Under the terms of Contract
-   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
-   certain rights in this software.  This software is distributed under
-   the GNU General Public License.
-
-   See the README file in the top-level LAMMPS directory.
-------------------------------------------------------------------------- */
-
-#include "math.h"
-#include "string.h"
-#include "stdlib.h"
-#include "fix_adapt.h"
-#include "atom.h"
-#include "update.h"
-#include "modify.h"
-#include "force.h"
-#include "pair.h"
-#include "pair_hybrid.h"
-#include "kspace.h"
-#include "input.h"
-#include "variable.h"
-#include "math_const.h"
-#include "memory.h"
-#include "error.h"
-
-using namespace LAMMPS_NS;
-using namespace FixConst;
-using namespace MathConst;
-
-enum{PAIR,KSPACE,ATOM};
-enum{DIAMETER,CHARGE};
-
-/* ---------------------------------------------------------------------- */
-
-FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
-{
-  if (narg < 5) error->all(FLERR,"Illegal fix adapt command");
-  nevery = force->inumeric(FLERR,arg[3]);
-  if (nevery < 0) error->all(FLERR,"Illegal fix adapt command");
-
-  // count # of adaptations
-
-  nadapt = 0;
-
-  int iarg = 4;
-  while (iarg < narg) {
-    if (strcmp(arg[iarg],"pair") == 0) {
-      if (iarg+6 > narg) error->all(FLERR,"Illegal fix adapt command");
-      nadapt++;
-      iarg += 6;
-    } else if (strcmp(arg[iarg],"kspace") == 0) {
-      if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt command");
-      nadapt++;
-      iarg += 2;
-    } else if (strcmp(arg[iarg],"atom") == 0) {
-      if (iarg+3 > narg) error->all(FLERR,"Illegal fix adapt command");
-      nadapt++;
-      iarg += 3;
-    } else break;
-  }
-
-  if (nadapt == 0) error->all(FLERR,"Illegal fix adapt command");
-  adapt = new Adapt[nadapt];
-
-  // parse keywords
-
-  nadapt = 0;
-  diamflag = 0;
-  chgflag = 0;
-
-  iarg = 4;
-  while (iarg < narg) {
-    if (strcmp(arg[iarg],"pair") == 0) {
-      if (iarg+6 > narg) error->all(FLERR,"Illegal fix adapt command");
-      adapt[nadapt].which = PAIR;
-      int n = strlen(arg[iarg+1]) + 1;
-      adapt[nadapt].pstyle = new char[n];
-      strcpy(adapt[nadapt].pstyle,arg[iarg+1]);
-      n = strlen(arg[iarg+2]) + 1;
-      adapt[nadapt].pparam = new char[n];
-      strcpy(adapt[nadapt].pparam,arg[iarg+2]);
-      force->bounds(arg[iarg+3],atom->ntypes,
-                    adapt[nadapt].ilo,adapt[nadapt].ihi);
-      force->bounds(arg[iarg+4],atom->ntypes,
-                    adapt[nadapt].jlo,adapt[nadapt].jhi);
-      if (strstr(arg[iarg+5],"v_") == arg[iarg+5]) {
-        n = strlen(&arg[iarg+5][2]) + 1;
-        adapt[nadapt].var = new char[n];
-        strcpy(adapt[nadapt].var,&arg[iarg+5][2]);
-      } else error->all(FLERR,"Illegal fix adapt command");
-      nadapt++;
-      iarg += 6;
-    } else if (strcmp(arg[iarg],"kspace") == 0) {
-      if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt command");
-      adapt[nadapt].which = KSPACE;
-      if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
-        int n = strlen(&arg[iarg+1][2]) + 1;
-        adapt[nadapt].var = new char[n];
-        strcpy(adapt[nadapt].var,&arg[iarg+1][2]);
-      } else error->all(FLERR,"Illegal fix adapt command");
-      nadapt++;
-      iarg += 2;
-    } else if (strcmp(arg[iarg],"atom") == 0) {
-      if (iarg+3 > narg) error->all(FLERR,"Illegal fix adapt command");
-      adapt[nadapt].which = ATOM;
-      if (strcmp(arg[iarg+1],"diameter") == 0) {
-        adapt[nadapt].aparam = DIAMETER;
-        diamflag = 1;
-      } else if (strcmp(arg[iarg+1],"charge") == 0) {
-        adapt[nadapt].aparam = CHARGE; 
-        chgflag = 1; 
-      } else error->all(FLERR,"Illegal fix adapt command");
-      if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) {
-        int n = strlen(&arg[iarg+2][2]) + 1;
-        adapt[nadapt].var = new char[n];
-        strcpy(adapt[nadapt].var,&arg[iarg+2][2]);
-      } else error->all(FLERR,"Illegal fix adapt command");
-      nadapt++;
-      iarg += 3;
-    } else break;
-  }
-
-  // optional keywords
-
-  resetflag = 0;
-  scaleflag = 0;
-
-  while (iarg < narg) {
-    if (strcmp(arg[iarg],"reset") == 0) {
-      if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt command");
-      if (strcmp(arg[iarg+1],"no") == 0) resetflag = 0;
-      else if (strcmp(arg[iarg+1],"yes") == 0) resetflag = 1;
-      else error->all(FLERR,"Illegal fix adapt command");
-      iarg += 2;
-    } else if (strcmp(arg[iarg],"scale") == 0) {
-      if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt command");
-      if (strcmp(arg[iarg+1],"no") == 0) scaleflag = 0;
-      else if (strcmp(arg[iarg+1],"yes") == 0) scaleflag = 1;
-      else error->all(FLERR,"Illegal fix adapt command");
-      iarg += 2;
-    } else error->all(FLERR,"Illegal fix adapt command");
-  }
-
-  // allocate pair style arrays
-
-  int n = atom->ntypes;
-  for (int m = 0; m < nadapt; m++)
-    if (adapt[m].which == PAIR)
-      memory->create(adapt[m].array_orig,n+1,n+1,"adapt:array_orig");
-}
-
-/* ---------------------------------------------------------------------- */
-
-FixAdapt::~FixAdapt()
-{
-  for (int m = 0; m < nadapt; m++) {
-    delete [] adapt[m].var;
-    if (adapt[m].which == PAIR) {
-      delete [] adapt[m].pstyle;
-      delete [] adapt[m].pparam;
-      memory->destroy(adapt[m].array_orig);
-    }
-  }
-  delete [] adapt;
-}
-
-/* ---------------------------------------------------------------------- */
-
-int FixAdapt::setmask()
-{
-  int mask = 0;
-  mask |= PRE_FORCE;
-  mask |= POST_RUN;
-  return mask;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixAdapt::init()
-{
-  int i,j;
-
-  // setup and error checks
-
-  anypair = 0;
-
-  for (int m = 0; m < nadapt; m++) {
-    Adapt *ad = &adapt[m];
-
-    ad->ivar = input->variable->find(ad->var);
-    if (ad->ivar < 0)
-      error->all(FLERR,"Variable name for fix adapt does not exist");
-    if (!input->variable->equalstyle(ad->ivar))
-      error->all(FLERR,"Variable for fix adapt is invalid style");
-
-    if (ad->which == PAIR) {
-      anypair = 1;
-
-      Pair *pair = force->pair_match(ad->pstyle,1);
-      if (pair == NULL) error->all(FLERR,"Fix adapt pair style does not exist");
-      void *ptr = pair->extract(ad->pparam,ad->pdim);
-      if (ptr == NULL) 
-        error->all(FLERR,"Fix adapt pair style param not supported");
-
-      ad->pdim = 2;
-      if (ad->pdim == 0) ad->scalar = (double *) ptr;
-      if (ad->pdim == 2) ad->array = (double **) ptr;
-
-      // if pair hybrid, test that ilo,ihi,jlo,jhi are valid for sub-style
-
-      if (ad->pdim == 2 && (strcmp(force->pair_style,"hybrid") == 0 ||
-                            strcmp(force->pair_style,"hybrid/overlay") == 0)) {
-        PairHybrid *pair = (PairHybrid *) force->pair;
-        for (i = ad->ilo; i <= ad->ihi; i++)
-          for (j = MAX(ad->jlo,i); j <= ad->jhi; j++)
-            if (!pair->check_ijtype(i,j,ad->pstyle))
-              error->all(FLERR,"Fix adapt type pair range is not valid for "
-                         "pair hybrid sub-style");
-      }
-
-    } else if (ad->which == KSPACE) {
-      if (force->kspace == NULL)
-        error->all(FLERR,"Fix adapt kspace style does not exist");
-      kspace_scale = (double *) force->kspace->extract("scale");
-
-    } else if (ad->which == ATOM) {
-      if (ad->aparam == DIAMETER) {
-        if (!atom->radius_flag)
-          error->all(FLERR,"Fix adapt requires atom attribute diameter");
-      }
-      if (ad->aparam == CHARGE) {
-	if (!atom->q_flag)
-	  error->all(FLERR,"Fix adapt requires atom attribute charge");
-      }
-    }
-  }
-
-  // make copy of original pair array values
-
-  for (int m = 0; m < nadapt; m++) {
-    Adapt *ad = &adapt[m];
-    if (ad->which == PAIR && ad->pdim == 2) {
-      for (i = ad->ilo; i <= ad->ihi; i++)
-        for (j = MAX(ad->jlo,i); j <= ad->jhi; j++)
-          ad->array_orig[i][j] = ad->array[i][j];
-    }
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixAdapt::setup_pre_force(int vflag)
-{
-  change_settings();
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixAdapt::pre_force(int vflag)
-{
-  if (nevery == 0) return;
-  if (update->ntimestep % nevery) return;
-  change_settings();
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixAdapt::post_run()
-{
-  if (resetflag) restore_settings();
-}
-
-/* ----------------------------------------------------------------------
-   change pair,kspace,atom parameters based on variable evaluation
-------------------------------------------------------------------------- */
-
-void FixAdapt::change_settings()
-{
-  int i,j;
-
-  // variable evaluation may invoke computes so wrap with clear/add
-
-  modify->clearstep_compute();
-
-  for (int m = 0; m < nadapt; m++) {
-    Adapt *ad = &adapt[m];
-    double value = input->variable->compute_equal(ad->ivar);
-
-    // set global scalar or type pair array values
-
-    if (ad->which == PAIR) {
-      if (ad->pdim == 0) {
-        if (scaleflag) *ad->scalar = value * ad->scalar_orig;
-        else *ad->scalar = value;
-      } else if (ad->pdim == 2) {
-        if (scaleflag)
-          for (i = ad->ilo; i <= ad->ihi; i++)
-            for (j = MAX(ad->jlo,i); j <= ad->jhi; j++)
-              ad->array[i][j] = value*ad->array_orig[i][j];
-        else
-          for (i = ad->ilo; i <= ad->ihi; i++)
-            for (j = MAX(ad->jlo,i); j <= ad->jhi; j++)
-              ad->array[i][j] = value;
-      }
-
-    // set kspace scale factor
-
-    } else if (ad->which == KSPACE) {
-      *kspace_scale = value;
-
-    } else if (ad->which == ATOM) {
-
-      // set radius from diameter
-      // also scale rmass to new value
-
-      if (ad->aparam == DIAMETER) {
-        int mflag = 0;
-        if (atom->rmass_flag) mflag = 1;
-        double density;
-
-        double *radius = atom->radius;
-        double *rmass = atom->rmass;
-        int *mask = atom->mask;
-        int nlocal = atom->nlocal;
-
-        if (mflag == 0) {
-          for (i = 0; i < nlocal; i++)
-            if (mask[i] & groupbit)
-              radius[i] = 0.5*value;
-        } else {
-          for (i = 0; i < nlocal; i++)
-            if (mask[i] & groupbit) {
-              density = rmass[i] / (4.0*MY_PI/3.0 *
-                                    radius[i]*radius[i]*radius[i]);
-              radius[i] = 0.5*value;
-              rmass[i] = 4.0*MY_PI/3.0 *
-                radius[i]*radius[i]*radius[i] * density;
-            }
-        }
-      } else if (ad->aparam == CHARGE) {
-        double *q = atom->q; 
-	int *mask = atom->mask;
-	int nlocal = atom->nlocal;
-	for (i = 0; i < nlocal; i++)
-	  if (mask[i] & groupbit) q[i] = value; 
-      }
-    }
-  }
-
-  modify->addstep_compute(update->ntimestep + nevery);
-
-  // re-initialize pair styles if any PAIR settings were changed
-  // this resets other coeffs that may depend on changed values,
-  // and also offset and tail corrections
-
-  if (anypair) force->pair->reinit();
-}
-
-/* ----------------------------------------------------------------------
-   restore pair,kspace.atom parameters to original values
-------------------------------------------------------------------------- */
-
-void FixAdapt::restore_settings()
-{
-  for (int m = 0; m < nadapt; m++) {
-    Adapt *ad = &adapt[m];
-    if (ad->which == PAIR) {
-      if (ad->pdim == 0) *ad->scalar = ad->scalar_orig;
-      else if (ad->pdim == 2) {
-        for (int i = ad->ilo; i <= ad->ihi; i++)
-          for (int j = MAX(ad->jlo,i); j <= ad->jhi; j++)
-            ad->array[i][j] = ad->array_orig[i][j];
-      }
-
-    } else if (ad->which == KSPACE) {
-      *kspace_scale = 1.0;
-
-    } else if (ad->which == ATOM) {
-
-    }
-  }
-
-  if (anypair) force->pair->reinit();
-}

src/fix_adapt.h

@@ -1,107 +0,0 @@
-/* ----------------------------------------------------------------------
-   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
-   http://lammps.sandia.gov, Sandia National Laboratories
-   Steve Plimpton, sjplimp@sandia.gov
-
-   Copyright (2003) Sandia Corporation.  Under the terms of Contract
-   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
-   certain rights in this software.  This software is distributed under
-   the GNU General Public License.
-
-   See the README file in the top-level LAMMPS directory.
-------------------------------------------------------------------------- */
-
-#ifdef FIX_CLASS
-
-FixStyle(adapt,FixAdapt)
-
-#else
-
-#ifndef LMP_FIX_ADAPT_H
-#define LMP_FIX_ADAPT_H
-
-#include "fix.h"
-
-namespace LAMMPS_NS {
-
-class FixAdapt : public Fix {
- public:
-  int diamflag;        // 1 if atom diameters will vary, for AtomVecGranular
-  int chgflag; 
-
-  FixAdapt(class LAMMPS *, int, char **);
-  ~FixAdapt();
-  int setmask();
-  void init();
-  void setup_pre_force(int);
-  void pre_force(int);
-  void post_run();
-
- private:
-  int nadapt,resetflag,scaleflag;
-  int anypair;
-
-  struct Adapt {
-    int which,ivar;
-    char *var;
-    char *pstyle,*pparam;
-    int ilo,ihi,jlo,jhi;
-    int pdim;
-    double *scalar,scalar_orig;
-    double **array,**array_orig;
-    int aparam;
-  };
-
-  Adapt *adapt;
-  double *kspace_scale;
-
-  void change_settings();
-  void restore_settings();
-};
-
-}
-
-#endif
-#endif
-
-/* ERROR/WARNING messages:
-
-E: Illegal ... command
-
-Self-explanatory.  Check the input script syntax and compare to the
-documentation for the command.  You can use -echo screen as a
-command-line option when running LAMMPS to see the offending line.
-
-E: Variable name for fix adapt does not exist
-
-Self-explanatory.
-
-E: Variable for fix adapt is invalid style
-
-Only equal-style variables can be used.
-
-E: Fix adapt pair style does not exist
-
-Self-explanatory
-
-E: Fix adapt pair style param not supported
-
-The pair style does not know about the parameter you specified.
-
-E: Fix adapt type pair range is not valid for pair hybrid sub-style
-
-Self-explanatory.
-
-E: Fix adapt kspace style does not exist
-
-Self-explanatory.
-
-E: Fix adapt requires atom attribute diameter
-
-The atom style being used does not specify an atom diameter.
-
-E: Fix adapt requires atom attribute charge
-
-The atom style being used does not specify an atom charge.
-
-*/