diff --git a/src/SHOCK/Install.sh b/src/SHOCK/Install.sh
index e2c94514e..642f0839c 100644
--- a/src/SHOCK/Install.sh
+++ b/src/SHOCK/Install.sh
@@ -1,15 +1,15 @@
# Install/unInstall package files in LAMMPS
if (test $1 = 1) then
- cp fix_msst.cpp ..
+ cp fix_msst.cpp fix_nphug.cpp ..
- cp fix_msst.h ..
+ cp fix_msst.h fix_nphug.h ..
elif (test $1 = 0) then
- rm -f ../fix_msst.cpp
+ rm -f ../fix_msst.cpp ../fix_nphug.cpp
- rm -f ../fix_msst.h
+ rm -f ../fix_msst.h ../fix_nphug.h
fi
diff --git a/src/SHOCK/fix_nphug.cpp b/src/SHOCK/fix_nphug.cpp
new file mode 100644
index 000000000..f22f67382
--- /dev/null
+++ b/src/SHOCK/fix_nphug.cpp
@@ -0,0 +1,224 @@
+/* ----------------------------------------------------------------------
+ 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 "string.h"
+#include "fix_nphug.h"
+#include "modify.h"
+#include "error.h"
+#include "update.h"
+#include "compute.h"
+#include "force.h"
+#include "domain.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+FixNPHug::FixNPHug(LAMMPS *lmp, int narg, char **arg) :
+ FixNH(lmp, narg, arg)
+{
+ // Hard-code to use initial state of system
+
+ v0_set = p0_set = e0_set = 0;
+
+ // Hard-code to use z direction
+
+ direction = 2;
+ if (p_flag[0] == 1 || p_flag[1] == 1 ||
+ p_flag[3] == 1 || p_flag[4] == 1 || p_flag[5] == 1)
+ error->all("Only pressure control in z direction to be used with fix nphug");
+ if (p_flag[2] == 0)
+ error->all("Pressure control in z direction must be used with fix nphug");
+
+ if (!tstat_flag)
+ error->all("Temperature control must be used with fix nphug");
+ if (!pstat_flag)
+ error->all("Pressure control must be used with fix nphug");
+
+ // 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;
+
+ // create a new compute potential energy compute
+
+ n = strlen(id) + 3;
+ 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;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixNPHug::~FixNPHug()
+{
+
+ // temp and press computes handled by base class
+ // delete pe compute
+
+ if (peflag) modify->delete_compute(id_pe);
+ delete [] id_pe;
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNPHug::init()
+{
+ // Call base class init()
+
+ FixNH::init();
+
+ // set pe ptr
+
+ int icompute = modify->find_compute(id_pe);
+ if (icompute < 0)
+ error->all("Potential energy ID for fix nvt/nph/npt does not exist");
+ pe = modify->compute[icompute];
+}
+
+
+/* ----------------------------------------------------------------------
+ compute initial state before integrator starts
+------------------------------------------------------------------------- */
+
+void FixNPHug::setup(int vflag)
+{
+ FixNH::setup(vflag);
+
+ if ( v0_set == 0 ) {
+ v0 = compute_vol();
+ v0_set = 1;
+ }
+
+ if ( p0_set == 0 ) {
+ p0 = p_current[direction];
+ p0_set = 1;
+ }
+
+ if ( e0_set == 0 ) {
+ e0 = compute_etotal();
+ e0_set = 1;
+ }
+
+}
+
+/* ----------------------------------------------------------------------
+ compute target temperature and kinetic energy
+-----------------------------------------------------------------------*/
+
+void FixNPHug::compute_temp_target()
+{
+ t_target = t_current + compute_hugoniot();
+ ke_target = tdof * boltz * t_target;
+ if (ke_target < 0.0) ke_target = 0.0;
+
+ // If t_target is very small, need to choose
+ // more reasonable value for use by barostat and
+ // thermostat masses. ke_target is left as is.
+
+ if (t_target <= 1.0e-6) {
+ if (strcmp(update->unit_style,"lj") == 0) t0 = 1.0;
+ else t0 = 300.0;
+ }
+ pressure->addstep(update->ntimestep+1);
+ pe->addstep(update->ntimestep+1);
+}
+
+
+/* ---------------------------------------------------------------------- */
+
+double FixNPHug::compute_etotal()
+{
+ double epot,ekin,etot;
+ epot = pe->compute_scalar();
+ if (thermo_energy) epot -= compute_scalar();
+ ekin = temperature->compute_scalar();
+ ekin *= 0.5 * temperature->dof * force->boltz;
+ etot = epot+ekin;
+ return etot;
+}
+
+/* ---------------------------------------------------------------------- */
+
+double FixNPHug::compute_vol()
+{
+ if (domain->dimension == 3)
+ return domain->xprd * domain->yprd * domain->zprd;
+ else
+ return domain->xprd * domain->yprd;
+}
+
+/* ----------------------------------------------------------------------
+ Computes the deviation of the current point
+ from the Hugoniot in energy units.
+------------------------------------------------------------------------- */
+
+double FixNPHug::compute_hugoniot()
+{
+ double v, e, p;
+ double dhugo;
+
+ e = compute_etotal();
+
+ temperature->compute_vector();
+ pressure->compute_vector();
+ p = pressure->vector[direction];
+
+ v = compute_vol();
+
+ dhugo = (0.5 * (p + p0 ) * ( v0 - v)) /
+ force->nktv2p + e0 - e;
+
+ return dhugo;
+}
diff --git a/src/SHOCK/fix_nphug.h b/src/SHOCK/fix_nphug.h
new file mode 100644
index 000000000..37ca8451e
--- /dev/null
+++ b/src/SHOCK/fix_nphug.h
@@ -0,0 +1,52 @@
+/* ----------------------------------------------------------------------
+ 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(nphug,FixNPHug)
+
+#else
+
+#ifndef LMP_FIX_NPHUG_H
+#define LMP_FIX_NPHUG_H
+
+#include "fix_nh.h"
+
+namespace LAMMPS_NS {
+
+class FixNPHug : public FixNH {
+ public:
+ FixNPHug(class LAMMPS *, int, char **);
+ ~FixNPHug();
+ void init();
+ void setup(int);
+
+ private:
+ class Compute *pe; // PE compute pointer
+
+ void compute_temp_target();
+ double compute_etotal();
+ double compute_vol();
+ double compute_hugoniot();
+
+ char *id_pe;
+ int peflag;
+ int v0_set,p0_set,e0_set;
+ double v0,p0,e0;
+ int direction;
+};
+
+}
+
+#endif
+#endif