diff --git a/src/.gitignore b/src/.gitignore
index 02b05b364b32e565abefb94aa5fdbaedded82396..0db6a21a754b20e4f9c8a7fe61cc1a41993c7ed2 100644
--- a/src/.gitignore
+++ b/src/.gitignore
@@ -205,6 +205,8 @@
 /compute_pe_tally.h
 /compute_plasticity_atom.cpp
 /compute_plasticity_atom.h
+/compute_pressure_grem.cpp
+/compute_pressure_grem.h
 /compute_rigid_local.cpp
 /compute_rigid_local.h
 /compute_spec_atom.cpp
@@ -343,6 +345,8 @@
 /fix_gle.h
 /fix_gpu.cpp
 /fix_gpu.h
+/fix_grem.cpp
+/fix_grem.h
 /fix_imd.cpp
 /fix_imd.h
 /fix_ipi.cpp
diff --git a/src/USER-MISC/compute_pressure_grem.cpp b/src/USER-MISC/compute_pressure_grem.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..387bde6427d6a9c2169c45fddcbaaf08dd1a311b
--- /dev/null
+++ b/src/USER-MISC/compute_pressure_grem.cpp
@@ -0,0 +1,149 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include <mpi.h>
+#include <string.h>
+#include <stdlib.h>
+#include "compute_pressure_grem.h"
+#include "atom.h"
+#include "update.h"
+#include "domain.h"
+#include "modify.h"
+#include "fix.h"
+#include "force.h"
+#include "pair.h"
+#include "kspace.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ----------------------------------------------------------------------
+   Last argument is the id of the gREM fix
+------------------------------------------------------------------------- */
+
+ComputePressureGrem::ComputePressureGrem(LAMMPS *lmp, int narg, char **arg) :
+  ComputePressure(lmp, narg-1, arg)
+{
+  fix_grem = arg[narg-1];
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePressureGrem::init()
+{
+  ComputePressure::init();
+
+  // Initialize hook to gREM fix
+  int ifix = modify->find_fix(fix_grem);
+  if (ifix < 0)
+    error->all(FLERR,"Fix grem ID for compute pressure/grem does not exist");
+
+  int dim;
+  scale_grem = (double *)modify->fix[ifix]->extract("scale_grem",dim);
+
+  if (scale_grem == NULL || dim != 0)
+    error->all(FLERR,"Cannot extract gREM scale factor from fix grem");
+}
+
+/* ----------------------------------------------------------------------
+   compute total pressure, averaged over Pxx, Pyy, Pzz
+------------------------------------------------------------------------- */
+
+double ComputePressureGrem::compute_scalar()
+{
+  invoked_scalar = update->ntimestep;
+  if (update->vflag_global != invoked_scalar)
+    error->all(FLERR,"Virial was not tallied on needed timestep");
+
+  // invoke temperature if it hasn't been already
+
+  double t;
+  if (keflag) {
+    if (temperature->invoked_scalar != update->ntimestep)
+      t = temperature->compute_scalar() * (*scale_grem);
+    else t = temperature->scalar * (*scale_grem);
+  }
+
+  if (dimension == 3) {
+    inv_volume = 1.0 / (domain->xprd * domain->yprd * domain->zprd);
+    virial_compute(3,3);
+    if (keflag)
+      scalar = (temperature->dof * boltz * t +
+                virial[0] + virial[1] + virial[2]) / 3.0 * inv_volume * nktv2p;
+    else
+      scalar = (virial[0] + virial[1] + virial[2]) / 3.0 * inv_volume * nktv2p;
+  } else {
+    inv_volume = 1.0 / (domain->xprd * domain->yprd);
+    virial_compute(2,2);
+    if (keflag)
+      scalar = (temperature->dof * boltz * t +
+                virial[0] + virial[1]) / 2.0 * inv_volume * nktv2p;
+    else
+      scalar = (virial[0] + virial[1]) / 2.0 * inv_volume * nktv2p;
+  }
+
+  return scalar;
+}
+
+/* ----------------------------------------------------------------------
+   compute pressure tensor
+   assume KE tensor has already been computed
+------------------------------------------------------------------------- */
+
+void ComputePressureGrem::compute_vector()
+{
+  invoked_vector = update->ntimestep;
+  if (update->vflag_global != invoked_vector)
+    error->all(FLERR,"Virial was not tallied on needed timestep");
+
+  if (force->kspace && kspace_virial && force->kspace->scalar_pressure_flag)
+    error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' for "
+	       "tensor components with kspace_style msm");
+
+  // invoke temperature if it hasn't been already
+
+  double *ke_tensor;
+  if (keflag) {
+    if (temperature->invoked_vector != update->ntimestep)
+      temperature->compute_vector();
+    ke_tensor = temperature->vector;
+  }
+  for (int i = 0; i < 6; i++)
+    ke_tensor[i] *= *scale_grem;
+
+  if (dimension == 3) {
+    inv_volume = 1.0 / (domain->xprd * domain->yprd * domain->zprd);
+    virial_compute(6,3);
+    if (keflag) {
+      for (int i = 0; i < 6; i++)
+        vector[i] = (ke_tensor[i] + virial[i]) * inv_volume * nktv2p;
+    } else
+      for (int i = 0; i < 6; i++)
+        vector[i] = virial[i] * inv_volume * nktv2p;
+  } else {
+    inv_volume = 1.0 / (domain->xprd * domain->yprd);
+    virial_compute(4,2);
+    if (keflag) {
+      vector[0] = (ke_tensor[0] + virial[0]) * inv_volume * nktv2p;
+      vector[1] = (ke_tensor[1] + virial[1]) * inv_volume * nktv2p;
+      vector[3] = (ke_tensor[3] + virial[3]) * inv_volume * nktv2p;
+      vector[2] = vector[4] = vector[5] = 0.0;
+    } else {
+      vector[0] = virial[0] * inv_volume * nktv2p;
+      vector[1] = virial[1] * inv_volume * nktv2p;
+      vector[3] = virial[3] * inv_volume * nktv2p;
+      vector[2] = vector[4] = vector[5] = 0.0;
+    }
+  }
+}
+
diff --git a/src/USER-MISC/compute_pressure_grem.h b/src/USER-MISC/compute_pressure_grem.h
new file mode 100644
index 0000000000000000000000000000000000000000..8d30af93222b9e72a9deb0bb6c39e5b7c27ac80a
--- /dev/null
+++ b/src/USER-MISC/compute_pressure_grem.h
@@ -0,0 +1,91 @@
+/* -*- 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 COMPUTE_CLASS
+
+ComputeStyle(pressure/grem,ComputePressureGrem)
+
+#else
+
+#ifndef LMP_COMPUTE_PRESSURE_GREM_H
+#define LMP_COMPUTE_PRESSURE_GREM_H
+
+#include "compute_pressure.h"
+
+namespace LAMMPS_NS {
+
+class ComputePressureGrem : public ComputePressure {
+ public:
+  ComputePressureGrem(class LAMMPS *, int, char **);
+  virtual ~ComputePressureGrem() {};
+  virtual void init();
+  virtual double compute_scalar();
+  virtual void compute_vector();
+
+ protected:
+  // Access to gREM fix scale factor
+  char   *fix_grem;
+  double *scale_grem;
+};
+
+}
+
+#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: Compute pressure must use group all
+
+Virial contributions computed by potentials (pair, bond, etc) are
+computed on all atoms.
+
+E: Could not find compute pressure temperature ID
+
+The compute ID for calculating temperature does not exist.
+
+E: Compute pressure temperature ID does not compute temperature
+
+The compute ID assigned to a pressure computation must compute
+temperature.
+
+E: Compute pressure requires temperature ID to include kinetic energy
+
+The keflag cannot be used unless a temperature compute is provided.
+
+E: Virial was not tallied on needed timestep
+
+You are using a thermo keyword that requires potentials to
+have tallied the virial, but they didn't on this timestep.  See the
+variable doc page for ideas on how to make this work.
+
+E: Must use 'kspace_modify pressure/scalar no' for tensor components with kspace_style msm
+
+Otherwise MSM will compute only a scalar pressure.  See the kspace_modify
+command for details on this setting.
+
+E: Fix grem ID for compute pressure/grem does not exist
+
+Compute pressure/grem was passed an invalid fix id
+
+E: Cannot extract gREM scale factor from fix grem
+
+The fix id passed to compute pressure/grem refers to an incompatible fix
+
+*/
diff --git a/src/USER-MISC/fix_grem.cpp b/src/USER-MISC/fix_grem.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7778052d23fbf39ad245eef7f1092914423694b5
--- /dev/null
+++ b/src/USER-MISC/fix_grem.cpp
@@ -0,0 +1,239 @@
+/* ----------------------------------------------------------------------
+   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.
+
+   Force scaling fix for gREM.
+   Cite: http://dx.doi.org/10.1063/1.3432176
+   Cite: http://dx.doi.org/10.1021/acs.jpcb.5b07614
+   
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Edyta Malolepsza (Broad Institute)
+   Contributing author: David Stelter (Boston University)
+   Contributing author: Tom Keyes (Boston University)
+------------------------------------------------------------------------- */
+
+#include <string.h>
+#include <stdlib.h>
+#include <math.h>
+#include "comm.h"
+#include "fix_grem.h"
+#include "atom.h"
+#include "force.h"
+#include "update.h"
+#include "modify.h"
+#include "domain.h"
+#include "input.h"
+#include "compute.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+enum{NONE,CONSTANT,EQUAL,ATOM};
+
+/* ---------------------------------------------------------------------- */
+
+FixGrem::FixGrem(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg)
+{
+  if (narg < 8) error->all(FLERR,"Illegal fix grem command");
+
+  scalar_flag = 1;
+  vector_flag = 1;
+  size_vector = 3;
+  global_freq = 1;
+  extscalar = 1;
+  extvector = 1;
+
+  // tbath - temp of bath, the same as defined in thermostat
+
+  lambda = force->numeric(FLERR,arg[3]);
+  eta = force->numeric(FLERR,arg[4]);
+  h0 = force->numeric(FLERR,arg[5]);
+  tbath = force->numeric(FLERR,arg[6]);
+  pressref = force->numeric(FLERR,arg[7]);
+
+  // 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*[5];
+  newarg[0] = id_press;
+  newarg[1] = (char *) "all";
+  newarg[2] = (char *) "pressure/grem";
+  newarg[3] = id_temp;
+  newarg[4] = id;
+  modify->add_compute(5,newarg);
+  delete [] newarg;
+  
+  // create a new compute ke style
+  // id = fix-ID + ke
+
+  n = strlen(id) + 8;
+  id_ke = new char[n];
+  strcpy(id_ke,id);
+  strcat(id_ke,"_ke");
+
+  newarg = new char*[3];
+  newarg[0] = id_ke;
+  newarg[1] = (char *) "all";
+  newarg[2] = (char *) "ke";
+  modify->add_compute(3,newarg);
+  delete [] newarg;
+  keflag = 1;
+
+  // create a new compute pe style
+  // id = fix-ID + pe
+
+  n = strlen(id) + 9;
+  id_pe = new char[n];
+  strcpy(id_pe,id);
+  strcat(id_pe,"_pe");
+
+  newarg = new char*[3];
+  newarg[0] = id_pe;
+  newarg[1] = (char *) "all";
+  newarg[2] = (char *) "pe";
+  modify->add_compute(3,newarg);
+  delete [] newarg;
+  peflag = 1;
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixGrem::~FixGrem()
+{
+  // delete temperature, pressure and energies if fix created them
+
+  if (tflag) modify->delete_compute(id_temp);
+  if (pflag) modify->delete_compute(id_press);
+  if (keflag) modify->delete_compute(id_ke);
+  if (peflag) modify->delete_compute(id_pe);
+  delete [] id_temp;
+  delete [] id_press;
+  delete [] id_ke;
+  delete [] id_pe;
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixGrem::setmask()
+{
+  int mask = 0;
+  mask |= POST_FORCE;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixGrem::init()
+{
+
+  // add check of sign of eta
+
+  // set temperature and pressure ptrs
+
+  int icompute = modify->find_compute(id_temp);
+  if (icompute < 0)
+    error->all(FLERR,"Temperature ID for fix scaledforce does not exist");
+  temperature = modify->compute[icompute];
+
+  icompute = modify->find_compute(id_ke);
+  if (icompute < 0)
+    error->all(FLERR,"Ke ID for fix scaledforce does not exist");
+  ke = modify->compute[icompute];
+
+  icompute = modify->find_compute(id_pe);
+  if (icompute < 0)
+    error->all(FLERR,"Pe ID for fix scaledforce does not exist");
+  pe = modify->compute[icompute];  
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixGrem::setup(int vflag)
+{
+  if (strstr(update->integrate_style,"verlet"))
+    post_force(vflag);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixGrem::min_setup(int vflag)
+{
+  post_force(vflag);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixGrem::post_force(int vflag)
+{
+  double **f = atom->f;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  
+  double tmpvolume = domain->xprd * domain->yprd * domain->zprd;
+  double tmppe = pe->compute_scalar();
+  // potential energy
+  double tmpenthalpy = tmppe+pressref*tmpvolume/(force->nktv2p);
+
+  double teffective = lambda+eta*(tmpenthalpy-h0);
+  scale_grem = tbath/teffective;
+  
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      f[i][0] *= scale_grem;
+      f[i][1] *= scale_grem;
+      f[i][2] *= scale_grem;
+    }
+  pe->addstep(update->ntimestep+1);
+}
+
+/* ----------------------------------------------------------------------
+   extract scale factor
+------------------------------------------------------------------------- */
+
+void *FixGrem::extract(const char *str, int &dim)
+{
+  dim=0;
+  if (strcmp(str,"scale_grem") == 0) {
+    return &scale_grem;
+  }
+  return NULL;
+}
diff --git a/src/USER-MISC/fix_grem.h b/src/USER-MISC/fix_grem.h
new file mode 100644
index 0000000000000000000000000000000000000000..e5087436b592be4ad7bfcc27e71b9847a355fd4b
--- /dev/null
+++ b/src/USER-MISC/fix_grem.h
@@ -0,0 +1,84 @@
+/* ----------------------------------------------------------------------
+   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(grem,FixGrem)
+
+#else
+
+#ifndef LMP_FIX_GREM_H
+#define LMP_FIX_GREM_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixGrem : public Fix {
+ public:
+  FixGrem(class LAMMPS *, int, char **);
+  ~FixGrem();
+  int setmask();
+  void init();
+  void setup(int);
+  void min_setup(int);
+  void post_force(int);
+  void *extract(const char *, int &);
+  double scale_grem;
+
+ private:
+  double lambda,eta,h0,tbath,pressref;
+
+ protected:
+  char *id_temp,*id_press,*id_ke,*id_pe;
+  class Compute *temperature,*pressure,*ke,*pe;
+  int pflag,tflag,keflag,peflag;
+
+};
+
+}
+
+#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: Region ID for fix grem does not exist
+
+Self-explanatory.
+
+E: Variable name for fix grem does not exist
+
+Self-explanatory.
+
+E: Variable for fix grem is invalid style
+
+Self-explanatory.
+
+E: Cannot use variable energy with constant force in fix grem
+
+This is because for constant force, LAMMPS can compute the change
+in energy directly.
+
+E: Must use variable energy with fix grem
+
+Must define an energy vartiable when applyting a dynamic
+force during minimization.
+
+*/
diff --git a/src/compute_pressure.h b/src/compute_pressure.h
index 4ff11683750941d7ddcdcdd808753df849e0f444..a59a64e634fbd5ed5c7194e5ca6703949e7dcb05 100644
--- a/src/compute_pressure.h
+++ b/src/compute_pressure.h
@@ -28,9 +28,9 @@ class ComputePressure : public Compute {
  public:
   ComputePressure(class LAMMPS *, int, char **);
   virtual ~ComputePressure();
-  void init();
-  double compute_scalar();
-  void compute_vector();
+  virtual void init();
+  virtual double compute_scalar();
+  virtual void compute_vector();
   void reset_extra_compute_fix(const char *);
 
  protected: