diff --git a/src/fix_box_relax.cpp b/src/fix_box_relax.cpp
new file mode 100644
index 0000000000..3042b86f02
--- /dev/null
+++ b/src/fix_box_relax.cpp
@@ -0,0 +1,444 @@
+/* ----------------------------------------------------------------------
+   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_box_relax.h"
+#include "atom.h"
+#include "domain.h"
+#include "update.h"
+#include "force.h"
+#include "modify.h"
+#include "compute.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+enum{XYZ,XY,YZ,XZ,ANISO};
+
+/* ---------------------------------------------------------------------- */
+
+FixBoxRelax::FixBoxRelax(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg)
+{
+  if (narg < 4) error->all("Illegal fix box/relax command");
+
+  box_change = 1;
+
+  if (strcmp(arg[3],"xyz") == 0) {
+    if (narg < 5) error->all("Illegal fix box/relax command");
+    press_couple = XYZ;
+    p_target[0] = p_target[1] = p_target[2] = atof(arg[4]);
+    p_flag[0] = p_flag[1] = p_flag[2] = 1;
+    if (domain->dimension == 2) p_flag[2] = 0;
+
+  } else {
+    if (strcmp(arg[3],"xy") == 0) press_couple = XY;
+    else if (strcmp(arg[3],"yz") == 0) press_couple = YZ;
+    else if (strcmp(arg[3],"xz") == 0) press_couple = XZ;
+    else if (strcmp(arg[3],"aniso") == 0) press_couple = ANISO;
+    else error->all("Illegal fix box/relax command");
+
+    if (narg < 7) error->all("Illegal fix box/relax command");
+
+    if (domain->dimension == 2 && 
+	(press_couple == XY || press_couple == YZ || press_couple == XZ))
+      error->all("Invalid fix box/relax command for a 2d simulation");
+
+    if (strcmp(arg[4],"NULL") == 0) p_flag[0] = 0;
+    else {
+      p_target[0] = atof(arg[4]);
+      p_flag[0] = 1;
+    }
+    if (strcmp(arg[5],"NULL") == 0) p_flag[1] = 0;
+    else {
+      p_target[1] = atof(arg[5]);
+      p_flag[1] = 1;
+    }
+    if (strcmp(arg[6],"NULL") == 0) p_flag[2] = 0;
+    else {
+      if (domain->dimension == 2)
+	error->all("Invalid fix box/relax command for a 2d simulation");
+      p_target[2] = atof(arg[6]);
+      p_flag[2] = 1;
+    }
+  }
+
+  // process extra keywords
+
+  allremap = 0;
+
+  int iarg;
+  if (press_couple == XYZ) iarg = 5;
+  else iarg = 7;
+
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"dilate") == 0) {
+      if (iarg+2 > narg) error->all("Illegal fix box/relax command");
+      if (strcmp(arg[iarg+1],"all") == 0) allremap = 0;
+      else if (strcmp(arg[iarg+1],"partial") == 0) allremap = 1;
+      else error->all("Illegal fix box/relax command");
+      iarg += 2;
+    } else error->all("Illegal fix box/relax command");
+  }
+
+  // error checks
+
+  if (press_couple == XY && (p_flag[0] == 0 || p_flag[1] == 0))
+    error->all("Invalid fix npt command");
+  if (press_couple == YZ && (p_flag[1] == 0 || p_flag[2] == 0))
+    error->all("Invalid fix npt command");
+  if (press_couple == XZ && (p_flag[0] == 0 || p_flag[2] == 0))
+    error->all("Invalid fix npt command");
+
+  if (press_couple == XY && p_target[0] != p_target[1])
+    error->all("Invalid fix npt command");
+  if (press_couple == YZ && p_target[1] != p_target[2])
+    error->all("Invalid fix npt command");
+  if (press_couple == XZ && p_target[0] != p_target[2])
+    error->all("Invalid fix npt command");
+
+  if (p_flag[0] && domain->xperiodic == 0)
+    error->all("Cannot use fix box/relax on a non-periodic dimension");
+  if (p_flag[1] && domain->yperiodic == 0)
+    error->all("Cannot use fix box/relax on a non-periodic dimension");
+  if (p_flag[2] && domain->zperiodic == 0)
+    error->all("Cannot use fix box/relax on a non-periodic dimension");
+
+  // create a new compute temp style
+  // id = fix-ID + temp
+  // compute group = all since pressure is always global (group all)
+  //   and thus its KE/temperature contribution should use group all
+
+  int n = strlen(id) + 6;
+  id_temp = new char[n];
+  strcpy(id_temp,id);
+  strcat(id_temp,"_temp");
+
+  char **newarg = new char*[3];
+  newarg[0] = id_temp;
+  newarg[1] = (char *) "all";
+  newarg[2] = (char *) "temp";
+  modify->add_compute(3,newarg);
+  delete [] newarg;
+  tflag = 1;
+
+  // create a new compute pressure style
+  // id = fix-ID + press, compute group = all
+  // pass id_temp as 4th arg to pressure constructor
+
+  n = strlen(id) + 7;
+  id_press = new char[n];
+  strcpy(id_press,id);
+  strcat(id_press,"_press");
+
+  newarg = new char*[4];
+  newarg[0] = id_press;
+  newarg[1] = (char *) "all";
+  newarg[2] = (char *) "pressure";
+  newarg[3] = id_temp;
+  modify->add_compute(4,newarg);
+  delete [] newarg;
+  pflag = 1;
+
+  nrigid = 0;
+  rfix = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixBoxRelax::~FixBoxRelax()
+{
+  delete [] rfix;
+
+  // delete temperature and pressure if fix created them
+
+  if (tflag) modify->delete_compute(id_temp);
+  if (pflag) modify->delete_compute(id_press);
+  delete [] id_temp;
+  delete [] id_press;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixBoxRelax::setmask()
+{
+  int mask = 0;
+  mask |= MIN_ENERGY;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixBoxRelax::init()
+{
+  // set temperature and pressure ptrs
+
+  int icompute = modify->find_compute(id_temp);
+  if (icompute < 0) error->all("Temp ID for fix npt does not exist");
+  temperature = modify->compute[icompute];
+
+  icompute = modify->find_compute(id_press);
+  if (icompute < 0) error->all("Press ID for fix npt does not exist");
+  pressure = modify->compute[icompute];
+
+  // initial box dimensions
+
+  xprdinit = domain->xprd;
+  yprdinit = domain->yprd;
+  zprdinit = domain->zprd;
+  if (domain->dimension == 3) volinit = domain->xprd*domain->yprd*domain->zprd;
+  else volinit = domain->xprd*domain->yprd;
+  pv2e = 1.0 / force->nktv2p;
+
+  // detect if any rigid fixes exist so rigid bodies move when box is remapped
+  // rfix[] = indices to each fix rigid
+
+  delete [] rfix;
+  nrigid = 0;
+  rfix = NULL;
+
+  for (int i = 0; i < modify->nfix; i++)
+    if (modify->fix[i]->rigid_flag) nrigid++;
+  if (nrigid) {
+    rfix = new int[nrigid];
+    nrigid = 0;
+    for (int i = 0; i < modify->nfix; i++)
+      if (modify->fix[i]->rigid_flag) rfix[nrigid++] = i;
+  }
+}
+
+/* ----------------------------------------------------------------------
+   compute energy and force due to extra degrees of freedom
+------------------------------------------------------------------------- */
+
+double FixBoxRelax::min_energy(double *fextra)
+{
+  double scale, scalex, scaley, scalez;
+
+  double t_current = temperature->compute_scalar();
+  if (press_couple == XYZ) {
+    double tmp = pressure->compute_scalar();
+  } else {
+    temperature->compute_vector();
+    pressure->compute_vector();
+  }
+  couple();
+
+  // trigger virial computation on every iteration of minimizer
+
+  pressure->addstep(update->ntimestep+1);
+
+  // compute new energy, forces for each extra degree of freedom
+  // PV must be in units of energy
+  
+  double eng;
+  if (press_couple == XYZ) {
+    scale = domain->xprd/xprdinit;
+    if (domain->dimension == 3) {
+      eng = pv2e * p_target[0] * (scale*scale*scale-1.0)*volinit;
+      fextra[0] = pv2e * (p_current[0] - p_target[0])*3.0*volinit*scale*scale;
+    } else {
+      eng = pv2e * p_target[0] * (scale*scale-1.0)*volinit;
+      fextra[0] = pv2e * (p_current[0] - p_target[0])*2.0*volinit*scale;
+    }
+  }
+  if (press_couple == ANISO) {
+    scalex = domain->xprd/xprdinit;
+    scaley = domain->yprd/yprdinit;
+    scalez = domain->zprd/zprdinit;
+    eng = pv2e * p_target[0] * (scalex*scaley*scalez-1.0)*volinit;
+    fextra[0] = fextra[1] = fextra[2] = 0.0;
+    if (p_flag[0]) {
+      fextra[0] = pv2e * (p_current[0] - p_target[0])*scaley*scalez*volinit;
+    }
+    if (p_flag[1]) {
+      fextra[1] = pv2e * (p_current[1] - p_target[1])*scalex*scalez*volinit;
+    }
+    if (p_flag[2]) {
+      fextra[2] = pv2e * (p_current[2] - p_target[2])*scalex*scaley*volinit;
+    }
+  }
+  if (press_couple == XY) {
+    scale = domain->xprd/xprdinit;
+    eng = pv2e * p_target[0] * (scale*scale-1.0)*volinit;
+    fextra[0] = fextra[1] = pv2e * (p_current[0] - p_target[0])*scale*volinit;
+    fextra[2] = 0.0;
+  }
+  if (press_couple == YZ) {
+    scale = domain->yprd/yprdinit;
+    eng = pv2e * p_target[1] * (scale*scale-1.0)*volinit;
+    fextra[1] = fextra[2] = pv2e * (p_current[1] - p_target[1])*scale*volinit;
+    fextra[0] = 0.0;
+  }
+  if (press_couple == XZ) {
+    scale = domain->xprd/xprdinit;
+    eng = pv2e * p_target[0] * (scale*scale-1.0)*volinit;
+    fextra[0] = fextra[2] = pv2e * (p_current[0] - p_target[0])*scale*volinit;
+    fextra[1] = 0.0;
+  }
+
+  return eng;
+}
+
+/* ----------------------------------------------------------------------
+   store extra dof values for linesearch starting point
+------------------------------------------------------------------------- */
+
+void FixBoxRelax::min_store()
+{
+  for (int i = 0; i < 3; i++) {
+    boxlo0[i] = domain->boxlo[i];
+    boxhi0[i] = domain->boxhi[i];
+  }
+  s0[0] = (boxhi0[0]-boxlo0[0])/xprdinit;
+  s0[1] = (boxhi0[1]-boxlo0[1])/yprdinit;
+  s0[2] = (boxhi0[2]-boxlo0[2])/zprdinit;
+}
+
+/* ----------------------------------------------------------------------
+   change the box dimensions by alpha from 
+   x0extra along linesearch direction fextra
+------------------------------------------------------------------------- */
+
+void FixBoxRelax::min_step(double alpha, double *fextra)
+{
+  double htmp;
+  if (press_couple == XYZ) {
+    ds[0] = alpha*fextra[0];
+    ds[1] = alpha*fextra[0];
+    ds[2] = alpha*fextra[0];
+  }
+  if (press_couple == ANISO) {
+    ds[0] = ds[1] = ds[2] = 0.0;
+    if (p_flag[0]) ds[0] = alpha*fextra[0];
+    if (p_flag[1]) ds[1] = alpha*fextra[1];
+    if (p_flag[2]) ds[2] = alpha*fextra[2];
+  }
+  if (press_couple == XY) {
+    ds[0] = ds[1] = alpha*fextra[0];
+    ds[2] = 0.0;
+  }
+  if (press_couple == XZ) {
+    ds[0] = ds[2] = alpha*fextra[0];
+    ds[1] = 0.0;
+  }
+  if (press_couple == YZ) {
+    ds[1] = ds[2] = alpha*fextra[1];
+    ds[0] = 0.0;
+  }
+
+  remap();
+}
+
+/* ----------------------------------------------------------------------
+   return number of degrees of freedom added by this fix
+------------------------------------------------------------------------- */
+
+int FixBoxRelax::min_dof()
+{
+  if (press_couple == XYZ) return 1;
+  else return 3;
+}
+
+/* ----------------------------------------------------------------------
+   dilate the box and owned/ghost atoms around center of box
+------------------------------------------------------------------------- */
+
+void FixBoxRelax::remap()
+{
+  int i,n;
+  double ctr;
+  
+  // ctr = geometric center of box in a dimension
+  // rescale simulation box from linesearch starting point
+  // scale atom coords for all atoms or only for fix group atoms
+  // don't do anything if non-periodic or style is constant volume
+
+  // we assume initially the atoms positions and box dimensions 
+  // are at linesearch starting point. 
+  // If not, atom positions will not be scaled correctly.
+
+  double **x = atom->x;
+  int *mask = atom->mask;
+  n = atom->nlocal + atom->nghost;
+
+  // convert pertinent atoms and rigid bodies to lamda coords
+
+  if (allremap) domain->x2lamda(n);
+  else {
+    for (i = 0; i < n; i++)
+      if (mask[i] & groupbit)
+	domain->x2lamda(x[i],x[i]);
+  }
+
+  if (nrigid)
+    for (i = 0; i < nrigid; i++)
+      modify->fix[rfix[i]]->deform(0);
+
+  // reset global and local box to new size/shape
+
+  for (i = 0; i < 3; i++) {
+    if (p_flag[i]) {
+      ctr = 0.5 * ( boxlo0[i] + boxhi0[i]);
+      domain->boxlo[i] = boxlo0[i] + (boxlo0[i]-ctr)*ds[i]/s0[i];
+      domain->boxhi[i] = boxhi0[i] + (boxhi0[i]-ctr)*ds[i]/s0[i];
+    }
+  }
+
+  domain->set_global_box();
+  domain->set_local_box();
+
+  // convert pertinent atoms and rigid bodies back to box coords
+
+  if (allremap) domain->lamda2x(n);
+  else {
+    for (i = 0; i < n; i++)
+      if (mask[i] & groupbit)
+	domain->lamda2x(x[i],x[i]);
+  }
+
+  if (nrigid)
+    for (i = 0; i < nrigid; i++)
+      modify->fix[rfix[i]]->deform(1);
+
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixBoxRelax::couple()
+{
+  double *tensor = pressure->vector;
+
+  if (press_couple == XYZ)
+    p_current[0] = p_current[1] = p_current[2] = pressure->scalar;
+  else if (press_couple == XY) {
+    double ave = 0.5 * (tensor[0] + tensor[1]);
+    p_current[0] = p_current[1] = ave;
+    p_current[2] = tensor[2];
+  } else if (press_couple == YZ) {
+    double ave = 0.5 * (tensor[1] + tensor[2]);
+    p_current[1] = p_current[2] = ave;
+    p_current[0] = tensor[0];
+  } else if (press_couple == XZ) {
+    double ave = 0.5 * (tensor[0] + tensor[2]);
+    p_current[0] = p_current[2] = ave;
+    p_current[1] = tensor[1];
+  } else if (press_couple == ANISO) {
+    p_current[0] = tensor[0];
+    p_current[1] = tensor[1];
+    p_current[2] = tensor[2];
+  }
+}
diff --git a/src/fix_box_relax.h b/src/fix_box_relax.h
new file mode 100644
index 0000000000..ef2074b5aa
--- /dev/null
+++ b/src/fix_box_relax.h
@@ -0,0 +1,58 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifndef FIX_BOX_RELAX_H
+#define FIX_BOX_RELAX_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixBoxRelax : public Fix {
+ public:
+  FixBoxRelax(class LAMMPS *, int, char **);
+  ~FixBoxRelax();
+  int setmask();
+  void init();
+  double min_energy(double *);
+  void min_store();
+  void min_step(double, double *);
+  int min_dof();
+
+ private:
+  int p_flag[3];
+  int press_couple,allremap;
+  double p_target[3],p_current[3];
+  double dilation[3];
+  double volinit,xprdinit,yprdinit,zprdinit;
+  double pv2e;
+
+  // These variables are for linesearch
+  double boxlo0[3],boxhi0[3];              // x0extra box bounds
+  double s0[6];	          	           // x0extra scale matrix in Voigt notation
+  double ds[6];                            // increment in scale matrix in Voigt notation
+
+  char *id_temp,*id_press;
+  class Compute *temperature,*pressure;
+  int tflag,pflag;
+
+  int nrigid;
+  int *rfix;
+
+  void remap();
+  void couple();
+};
+
+}
+
+#endif