intial refactoring of THERMO_ENERGY mask

src/RIGID/fix_shake.cpp

@@ -33,8 +33,6 @@
 #include "memory.h"
 #include "error.h"
 
-
-
 using namespace LAMMPS_NS;
 using namespace FixConst;
 using namespace MathConst;

src/compute_pe.cpp

@@ -99,7 +99,8 @@ double ComputePE::compute_scalar()
     scalar += force->pair->etail / volume;
   }
 
-  if (fixflag && modify->n_thermo_energy) scalar += modify->thermo_energy();
+  if (fixflag && modify->n_energy_global)
+    scalar += modify->energy_global();
 
   return scalar;
 }

src/compute_pe_atom.cpp

@@ -152,8 +152,8 @@ void ComputePEAtom::compute_peratom()
   // add in per-atom contributions from relevant fixes
   // always only for owned atoms, not ghost
 
-  if (fixflag && modify->n_thermo_energy_atom)
+  if (fixflag && modify->n_energy_atom)
-    modify->thermo_energy_atom(nlocal,energy);
+    modify->energy_atom(nlocal,energy);
 
   // communicate ghost energy between neighbor procs
 

src/fix.cpp

@@ -63,9 +63,11 @@ Fix::Fix(LAMMPS *lmp, int /*narg*/, char **arg) :
   box_change = NO_BOX_CHANGE;
   thermo_energy = 0;
   thermo_virial = 0;
+  energy_global_flag = energy_atom_flag = 0;
+  virial_global_flag = virial_atom_flag = 0;
+  ecouple_flag = 0;
   rigid_flag = 0;
   peatom_flag = 0;
-  virial_flag = 0;
   no_change_box = 0;
   time_integrate = 0;
   time_depend = 0;
@@ -150,8 +152,7 @@ void Fix::modify_params(int narg, char **arg)
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command");
       if (strcmp(arg[iarg+1],"no") == 0) thermo_energy = 0;
       else if (strcmp(arg[iarg+1],"yes") == 0) {
-        if (!(THERMO_ENERGY & setmask()))
+        if (energy_flag == 0) error->all(FLERR,"Illegal fix_modify command");
-          error->all(FLERR,"Illegal fix_modify command");
         thermo_energy = 1;
       } else error->all(FLERR,"Illegal fix_modify command");
       iarg += 2;
@@ -159,8 +160,7 @@ void Fix::modify_params(int narg, char **arg)
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command");
       if (strcmp(arg[iarg+1],"no") == 0) thermo_virial = 0;
       else if (strcmp(arg[iarg+1],"yes") == 0) {
-        if (virial_flag == 0)
+        if (virial_flag == 0) error->all(FLERR,"Illegal fix_modify command");
-          error->all(FLERR,"Illegal fix_modify command");
         thermo_virial = 1;
       } else error->all(FLERR,"Illegal fix_modify command");
       iarg += 2;
@@ -180,9 +180,12 @@ void Fix::modify_params(int narg, char **arg)
 }
 
 /* ----------------------------------------------------------------------
-   setup for energy, virial computation
+   setup for peratom energy and global/peratom virial computation
    see integrate::ev_set() for values of eflag (0-3) and vflag (0-6)
-   fixes call this if they use ev_tally()
+   fixes call Fix::ev_init() if tally energy and virial values
+   if thermo_energy is not set, energy tallying is disabled
+   if thermo_virial is not set, virial tallying is disabled
+   global energy is tallied separately, output by compute_scalar() method
 ------------------------------------------------------------------------- */
 
 void Fix::ev_setup(int eflag, int vflag)
@@ -191,13 +194,19 @@ void Fix::ev_setup(int eflag, int vflag)
 
   evflag = 1;
 
-  eflag_either = eflag;
+  if (!thermo_energy) eflag_either = eflag_global = eflag_atom = 0;
-  eflag_global = eflag & ENERGY_GLOBAL;
+  else {
-  eflag_atom = eflag & ENERGY_ATOM;
+    eflag_either = eflag;
+    eflag_global = eflag & ENERGY_GLOBAL;
+    eflag_atom = eflag & ENERGY_ATOM;
+  }
 
-  vflag_either = vflag;
+  if (!thermo_virial) vflag_either = vflag_global = vflag_atom = 0;
-  vflag_global = vflag & (VIRIAL_PAIR | VIRIAL_FDOTR);
+  else {
-  vflag_atom = vflag & (VIRIAL_ATOM | VIRIAL_CENTROID);
+    vflag_either = vflag;
+    vflag_global = vflag & (VIRIAL_PAIR | VIRIAL_FDOTR);
+    vflag_atom = vflag & (VIRIAL_ATOM | VIRIAL_CENTROID);
+  }
 
   // reallocate per-atom arrays if necessary
 
@@ -235,24 +244,17 @@ void Fix::ev_setup(int eflag, int vflag)
 }
 
 /* ----------------------------------------------------------------------
-   if thermo_virial is on:
+   setup for global/peratom virial computation
-     setup for virial computation
+   see integrate::ev_set() for values of vflag (0-6)
-     see integrate::ev_set() for values of vflag
+   fixes call Fix::v_init() if tally virial values but not energy
-     fixes call this if use v_tally()
+   if thermo_virial is not set, virial tallying is disabled
-   else: set evflag=0
 ------------------------------------------------------------------------- */
 
 void Fix::v_setup(int vflag)
 {
   int i,n;
 
-  if (!thermo_virial) {
-    evflag = 0;
-    return;
-  }
-
   evflag = 1;
-
   vflag_global = vflag & (VIRIAL_PAIR | VIRIAL_FDOTR);
   vflag_atom = vflag & (VIRIAL_ATOM | VIRIAL_CENTROID);
 
@@ -304,7 +306,6 @@ void Fix::ev_tally(int n, int *list, double total, double eng, double *v)
   v_tally(n,list,total,v);
 }
 
-
 /* ----------------------------------------------------------------------
    tally virial into global and per-atom accumulators
    n = # of local owned atoms involved, with local indices in list

src/fix.h

@@ -40,14 +40,17 @@ class Fix : protected Pointers {
   };
 
   bigint next_reneighbor;        // next timestep to force a reneighboring
-  int thermo_energy;             // 1 if fix_modify enabled ThEng, 0 if not
-  int thermo_virial;             // 1 if fix_modify enabled ThVir, 0 if not
   int nevery;                    // how often to call an end_of_step fix
-  int rigid_flag;                // 1 if Fix integrates rigid bodies, 0 if not
+  int thermo_energy;             // 1 if fix_modify energy enabled, 0 if not
-  int peatom_flag;               // 1 if Fix contributes per-atom eng, 0 if not
+  int thermo_virial;             // 1 if fix_modify virial enabled, 0 if not
-  int virial_flag;               // 1 if Fix contributes to virial, 0 if not
+  int energy_global_flag;        // 1 if contributes to global eng
+  int energy_atom_flag;          // 1 if contributes to peratom eng
+  int virial_global_flag;        // 1 if contributes to global virial
+  int virial_peratom_flag;       // 1 if contributes to peratom virial
+  int ecouple_flag;              // 1 if thermostat which accumulates eng to ecouple
+  int time_integrate;            // 1 if performs time integration, 0 if no
+  int rigid_flag;                // 1 if integrates rigid bodies, 0 if not
   int no_change_box;             // 1 if cannot swap ortho <-> triclinic
-  int time_integrate;            // 1 if fix performs time integration, 0 if no
   int time_depend;               // 1 if requires continuous timestepping
   int create_attribute;          // 1 if fix stores attributes that need
                                  //      setting when a new atom is created
@@ -99,8 +102,9 @@ class Fix : protected Pointers {
   int comm_reverse;              // size of reverse communication (0 if none)
   int comm_border;               // size of border communication (0 if none)
 
-  double virial[6];              // accumulated virial
+  double ecouple;            // cumulative energy added to reservoir by thermostatting
-  double *eatom,**vatom;         // accumulated per-atom energy/virial
+  double virial[6];          // virial for this timestep
+  double *eatom,**vatom;     // per-atom energy/virial for this timestep
 
   int centroidstressflag;        // centroid stress compared to two-body stress
                                  // CENTROID_SAME = same as two-body stress
@@ -239,11 +243,17 @@ class Fix : protected Pointers {
   int dynamic;    // recount atoms for temperature computes
 
   void ev_init(int eflag, int vflag) {
-    if (eflag||vflag) ev_setup(eflag, vflag);
+    if ((eflag && thermo_energy) || (vflag && thermo_virial)) ev_setup(eflag, vflag);
-    else evflag = eflag_either = eflag_global = eflag_atom = vflag_either = vflag_global = vflag_atom = 0;
+    else evflag = eflag_either = eflag_global = eflag_atom =
+	 vflag_either = vflag_global = vflag_atom = 0;
   }
   void ev_setup(int, int);
   void ev_tally(int, int *, double, double, double *);
+
+  void v_init(int vflag) {
+    if (vflag && thermo_virial) v_setup(vflag);
+    else evflag = vflag_either = vflag_global = vflag_atom = 0;
+  }
   void v_setup(int);
   void v_tally(int, int *, double, double *);
   void v_tally(int, double *);
@@ -263,7 +273,7 @@ namespace FixConst {
     FINAL_INTEGRATE =         1<<8,
     END_OF_STEP =             1<<9,
     POST_RUN =                1<<10,
-    THERMO_ENERGY =           1<<11,
+    //THERMO_ENERGY =           1<<11,   // remove when done with refactoring
     INITIAL_INTEGRATE_RESPA = 1<<12,
     POST_INTEGRATE_RESPA =    1<<13,
     PRE_FORCE_RESPA =         1<<14,

src/fix_addforce.cpp

@@ -352,7 +352,7 @@ void FixAddForce::post_force(int vflag)
           v[3] = xstyle ? xvalue*unwrap[1] : 0.0;
           v[4] = xstyle ? xvalue*unwrap[2] : 0.0;
           v[5] = ystyle ? yvalue*unwrap[2] : 0.0;
-          v_tally(i, v);
+          v_tally(i,v);
         }
       }
     }

src/fix_temp_rescale.cpp

@@ -12,9 +12,10 @@
 ------------------------------------------------------------------------- */
 
 #include "fix_temp_rescale.h"
-#include <cstring>
 
+#include <cstring>
 #include <cmath>
+
 #include "atom.h"
 #include "force.h"
 #include "group.h"
@@ -26,7 +27,6 @@
 #include "compute.h"
 #include "error.h"
 
-
 using namespace LAMMPS_NS;
 using namespace FixConst;
 

src/modify.cpp

@@ -46,8 +46,8 @@ Modify::Modify(LAMMPS *lmp) : Pointers(lmp)
   n_initial_integrate = n_post_integrate = 0;
   n_pre_exchange = n_pre_neighbor = n_post_neighbor = 0;
   n_pre_force = n_pre_reverse = n_post_force = 0;
-  n_final_integrate = n_end_of_step = n_thermo_energy = 0;
+  n_final_integrate = n_end_of_step = 0;
-  n_thermo_energy_atom = 0;
+  n_energy_couple = n_energy_global = n_energy_atom = 0;
   n_initial_integrate_respa = n_post_integrate_respa = 0;
   n_pre_force_respa = n_post_force_respa = n_final_integrate_respa = 0;
   n_min_pre_exchange = n_min_pre_force = n_min_pre_reverse = 0;
@@ -60,7 +60,7 @@ Modify::Modify(LAMMPS *lmp) : Pointers(lmp)
   list_pre_exchange = list_pre_neighbor = list_post_neighbor = nullptr;
   list_pre_force = list_pre_reverse = list_post_force = nullptr;
   list_final_integrate = list_end_of_step = nullptr;
-  list_thermo_energy = list_thermo_energy_atom = nullptr;
+  list_energy_couple = list_energy_global = list_energy_atom = nullptr;
   list_initial_integrate_respa = list_post_integrate_respa = nullptr;
   list_pre_force_respa = list_post_force_respa = nullptr;
   list_final_integrate_respa = nullptr;
@@ -137,8 +137,9 @@ Modify::~Modify()
   delete [] list_post_force;
   delete [] list_final_integrate;
   delete [] list_end_of_step;
-  delete [] list_thermo_energy;
+  delete [] list_energy_couple;
-  delete [] list_thermo_energy_atom;
+  delete [] list_energy_global;
+  delete [] list_energy_atom;
   delete [] list_initial_integrate_respa;
   delete [] list_post_integrate_respa;
   delete [] list_pre_force_respa;
@@ -218,8 +219,9 @@ void Modify::init()
   list_init(POST_FORCE,n_post_force,list_post_force);
   list_init(FINAL_INTEGRATE,n_final_integrate,list_final_integrate);
   list_init_end_of_step(END_OF_STEP,n_end_of_step,list_end_of_step);
-  list_init_thermo_energy(THERMO_ENERGY,n_thermo_energy,list_thermo_energy);
+  list_init_energy_couple(n_energy_couple,list_energy_couple);
-  list_init_thermo_energy_atom(n_thermo_energy_atom,list_thermo_energy_atom);
+  list_init_energy_global(n_energy_global,list_energy_global);
+  list_init_energy_atom(n_energy_atom,list_energy_atom);
 
   list_init(INITIAL_INTEGRATE_RESPA,
             n_initial_integrate_respa,list_initial_integrate_respa);
@@ -486,31 +488,45 @@ void Modify::end_of_step()
 }
 
 /* ----------------------------------------------------------------------
-   thermo energy call, only for relevant fixes
+   coupling energy call, only for relevant fixes
-   called by Thermo class
+   stored by each fix in ecouple variable
-   compute_scalar() is fix call to return energy
+   ecouple = cumulative energy added to reservoir by thermostatting
 ------------------------------------------------------------------------- */
 
-double Modify::thermo_energy()
+double Modify::energy_couple()
 {
   double energy = 0.0;
-  for (int i = 0; i < n_thermo_energy; i++)
+  for (int i = 0; i < n_energy_couple; i++)
-    energy += fix[list_thermo_energy[i]]->compute_scalar();
+    energy += fix[list_energy_couple[i]]->ecouple;
   return energy;
 }
 
 /* ----------------------------------------------------------------------
-   per-atom thermo energy call, only for relevant fixes
+   global energy call, only for relevant fixes
+   they return energy via compute_scalar()
+   called by compute pe
+------------------------------------------------------------------------- */
+
+double Modify::energy_global()
+{
+  double energy = 0.0;
+  for (i = 0; i < n_energy_global; i++)
+    energy += fix[list_energy_global[i]]->compute_scalar();
+  return energy;
+}
+
+/* ----------------------------------------------------------------------
+   peratom energy call, only for relevant fixes
    called by compute pe/atom
 ------------------------------------------------------------------------- */
 
-void Modify::thermo_energy_atom(int nlocal, double *energy)
+void Modify::energy_atom(int nlocal, double *energy)
 {
   int i,j;
   double *eatom;
 
-  for (i = 0; i < n_thermo_energy_atom; i++) {
+  for (i = 0; i < n_energy_atom; i++) {
-    eatom = fix[list_thermo_energy_atom[i]]->eatom;
+    eatom = fix[list_energy_atom[i]]->eatom;
     if (!eatom) continue;
     for (j = 0; j < nlocal; j++) energy[j] += eatom[j];
   }
@@ -1632,43 +1648,79 @@ void Modify::list_init_end_of_step(int mask, int &n, int *&list)
 }
 
 /* ----------------------------------------------------------------------
-   create list of fix indices for thermo energy fixes
+   create list of fix indices for fixes that compute reservoir coupling energy
-   only added to list if fix has THERMO_ENERGY mask set,
+   only added to list if fix has ecouple_flag set
-   and its thermo_energy flag was set via fix_modify
 ------------------------------------------------------------------------- */
 
-void Modify::list_init_thermo_energy(int mask, int &n, int *&list)
+void Modify::list_init_energy_couple(int &n, int *&list)
 {
   delete [] list;
 
   n = 0;
   for (int i = 0; i < nfix; i++)
-    if (fmask[i] & mask && fix[i]->thermo_energy) n++;
+    if (fix[i]->ecouple_flag) n++;
   list = new int[n];
 
   n = 0;
   for (int i = 0; i < nfix; i++)
-    if (fmask[i] & mask && fix[i]->thermo_energy) list[n++] = i;
+    if (fix[i]->ecouple_flag) list[n++] = i;
 }
 
 /* ----------------------------------------------------------------------
-   create list of fix indices for peratom thermo energy fixes
+   create list of fix indices for fixes that compute peratom energy
-   only added to list if fix has its peatom_flag set,
+   only added to list if fix has energy_atom_flag and thermo_energy set
-   and its thermo_energy flag was set via fix_modify
 ------------------------------------------------------------------------- */
 
-void Modify::list_init_thermo_energy_atom(int &n, int *&list)
+void Modify::list_init_energy_atom(int &n, int *&list)
 {
   delete [] list;
 
   n = 0;
   for (int i = 0; i < nfix; i++)
-    if (fix[i]->peatom_flag && fix[i]->thermo_energy) n++;
+    if (fix[i]->energy_atom_flag && fix[i]->thermo_energy) n++;
   list = new int[n];
 
   n = 0;
   for (int i = 0; i < nfix; i++)
-    if (fix[i]->peatom_flag && fix[i]->thermo_energy) list[n++] = i;
+    if (fix[i]->energy_atom_flag && fix[i]->thermo_energy) list[n++] = i;
+}
+
+/* ----------------------------------------------------------------------
+   create list of fix indices for fixes that compute global energy
+   only added to list if fix has energy_global_flag and thermo_energy set
+------------------------------------------------------------------------- */
+
+void Modify::list_init_energy_global(int &n, int *&list)
+{
+  delete [] list;
+
+  n = 0;
+  for (int i = 0; i < nfix; i++)
+    if (fix[i]->energy_global_flag && fix[i]->thermo_energy) n++;
+  list = new int[n];
+
+  n = 0;
+  for (int i = 0; i < nfix; i++)
+    if (fix[i]->energy_global_flag && fix[i]->thermo_energy) list[n++] = i;
+}
+
+/* ----------------------------------------------------------------------
+   create list of fix indices for fixes that compute peratom energy
+   only added to list if fix has energy_atom_flag and thermo_energy set
+------------------------------------------------------------------------- */
+
+void Modify::list_init_energy_atom(int &n, int *&list)
+{
+  delete [] list;
+
+  n = 0;
+  for (int i = 0; i < nfix; i++)
+    if (fix[i]->energy_atom_flag && fix[i]->thermo_energy) n++;
+  list = new int[n];
+
+  n = 0;
+  for (int i = 0; i < nfix; i++)
+    if (fix[i]->energy_atom_flag && fix[i]->thermo_energy) list[n++] = i;
 }
 
 /* ----------------------------------------------------------------------

src/modify.h

@@ -33,7 +33,8 @@ class Modify : protected Pointers {
   int n_initial_integrate,n_post_integrate,n_pre_exchange;
   int n_pre_neighbor,n_post_neighbor;
   int n_pre_force,n_pre_reverse,n_post_force;
-  int n_final_integrate,n_end_of_step,n_thermo_energy,n_thermo_energy_atom;
+  int n_final_integrate,n_end_of_step;
+  int n_energy_couple,n_energy_global,n_energy_atom;
   int n_initial_integrate_respa,n_post_integrate_respa;
   int n_pre_force_respa,n_post_force_respa,n_final_integrate_respa;
   int n_min_pre_exchange,n_min_pre_neighbor,n_min_post_neighbor;
@@ -69,8 +70,9 @@ class Modify : protected Pointers {
   virtual void post_force(int);
   virtual void final_integrate();
   virtual void end_of_step();
-  virtual double thermo_energy();
+  virtual double energy_couple();
-  virtual void thermo_energy_atom(int, double *);
+  virtual double energy_global();
+  virtual void energy_atom(int, double *);
   virtual void post_run();
   virtual void create_attribute(int);
 
@@ -135,8 +137,8 @@ class Modify : protected Pointers {
   int *list_initial_integrate,*list_post_integrate;
   int *list_pre_exchange,*list_pre_neighbor,*list_post_neighbor;
   int *list_pre_force,*list_pre_reverse,*list_post_force;
-  int *list_final_integrate,*list_end_of_step,*list_thermo_energy;
+  int *list_final_integrate,*list_end_of_step;
-  int *list_thermo_energy_atom;
+  int *list_energy_couple,*list_energy_global,*list_energy_atom;
   int *list_initial_integrate_respa,*list_post_integrate_respa;
   int *list_pre_force_respa,*list_post_force_respa;
   int *list_final_integrate_respa;
@@ -163,8 +165,9 @@ class Modify : protected Pointers {
 
   void list_init(int, int &, int *&);
   void list_init_end_of_step(int, int &, int *&);
-  void list_init_thermo_energy(int, int &, int *&);
+  void list_init_energy_couple(int &, int *&);
-  void list_init_thermo_energy_atom(int &, int *&);
+  void list_init_energy_global(int &, int *&);
+  void list_init_energy_atom(int &, int *&);
   void list_init_dofflag(int &, int *&);
   void list_init_compute();
 

src/thermo.cpp

@@ -54,13 +54,14 @@ using namespace MathConst;
 
 // step, elapsed, elaplong, dt, time, cpu, tpcpu, spcpu, cpuremain,
 // part, timeremain
-// atoms, temp, press, pe, ke, etotal, enthalpy
+// atoms, temp, press, pe, ke, etotal
 // evdwl, ecoul, epair, ebond, eangle, edihed, eimp, emol, elong, etail
-// vol, density, lx, ly, lz, xlo, xhi, ylo, yhi, zlo, zhi, xy, xz, yz,
+// enthalpy, ecouple, econserve
+// vol, density, lx, ly, lz, xlo, xhi, ylo, yhi, zlo, zhi, xy, xz, yz
 // xlat, ylat, zlat
-// bonds, angles, dihedrals, impropers,
+// bonds, angles, dihedrals, impropers
 // pxx, pyy, pzz, pxy, pxz, pyz
-// fmax, fnorm, nbuild, ndanger,
+// fmax, fnorm, nbuild, ndanger
 // cella, cellb, cellc, cellalpha, cellbeta, cellgamma
 
 // customize a new thermo style by adding a DEFINE to this list
@@ -753,11 +754,6 @@ void Thermo::parse_fields(char *str)
       addfield("TotEng",&Thermo::compute_etotal,FLOAT);
       index_temp = add_compute(id_temp,SCALAR);
       index_pe = add_compute(id_pe,SCALAR);
-    } else if (word == "enthalpy") {
-      addfield("Enthalpy",&Thermo::compute_enthalpy,FLOAT);
-      index_temp = add_compute(id_temp,SCALAR);
-      index_press_scalar = add_compute(id_press,SCALAR);
-      index_pe = add_compute(id_pe,SCALAR);
 
     } else if (word == "evdwl") {
       addfield("E_vdwl",&Thermo::compute_evdwl,FLOAT);
@@ -790,6 +786,19 @@ void Thermo::parse_fields(char *str)
       addfield("E_tail",&Thermo::compute_etail,FLOAT);
       index_pe = add_compute(id_pe,SCALAR);
 
+    } else if (word == "enthalpy") {
+      addfield("Enthalpy",&Thermo::compute_enthalpy,FLOAT);
+      index_temp = add_compute(id_temp,SCALAR);
+      index_press_scalar = add_compute(id_press,SCALAR);
+      index_pe = add_compute(id_pe,SCALAR);
+    } else if (word == "ecouple") {
+      addfield("Ecouple",&Thermo::compute_ecouple,FLOAT);
+      index_pe = add_compute(id_pe,SCALAR);
+    } else if (word == "econserve") {
+      addfield("Econserve",&Thermo::compute_econserve,FLOAT);
+      index_temp = add_compute(id_temp,SCALAR);
+      index_pe = add_compute(id_pe,SCALAR);
+
     } else if (word == "vol") {
       addfield("Volume",&Thermo::compute_vol,FLOAT);
     } else if (word == "density") {
@@ -1234,42 +1243,6 @@ int Thermo::evaluate_keyword(const char *word, double *answer)
     }
     compute_etotal();
 
-  } else if (strcmp(word,"enthalpy") == 0) {
-    if (!pe)
-      error->all(FLERR,
-                 "Thermo keyword in variable requires thermo to use/init pe");
-    if (update->whichflag == 0) {
-      if (pe->invoked_scalar != update->ntimestep)
-        error->all(FLERR,"Compute used in variable thermo keyword between runs "
-                   "is not current");
-    } else {
-      pe->compute_scalar();
-      pe->invoked_flag |= INVOKED_SCALAR;
-    }
-    if (!temperature)
-      error->all(FLERR,"Thermo keyword in variable requires "
-                 "thermo to use/init temp");
-    if (update->whichflag == 0) {
-      if (temperature->invoked_scalar != update->ntimestep)
-        error->all(FLERR,"Compute used in variable thermo keyword between runs "
-                   "is not current");
-    } else if (!(temperature->invoked_flag & INVOKED_SCALAR)) {
-      temperature->compute_scalar();
-      temperature->invoked_flag |= INVOKED_SCALAR;
-    }
-    if (!pressure)
-      error->all(FLERR,"Thermo keyword in variable requires "
-                 "thermo to use/init press");
-    if (update->whichflag == 0) {
-      if (pressure->invoked_scalar != update->ntimestep)
-        error->all(FLERR,"Compute used in variable thermo keyword between runs "
-                   "is not current");
-    } else if (!(pressure->invoked_flag & INVOKED_SCALAR)) {
-      pressure->compute_scalar();
-      pressure->invoked_flag |= INVOKED_SCALAR;
-    }
-    compute_enthalpy();
-
   } else if (strcmp(word,"evdwl") == 0) {
     if (update->eflag_global != update->ntimestep)
       error->all(FLERR,"Energy was not tallied on needed timestep");
@@ -1356,6 +1329,68 @@ int Thermo::evaluate_keyword(const char *word, double *answer)
       error->all(FLERR,"Energy was not tallied on needed timestep");
     compute_etail();
 
+  } else if (strcmp(word,"enthalpy") == 0) {
+    if (!pe)
+      error->all(FLERR,
+                 "Thermo keyword in variable requires thermo to use/init pe");
+    if (update->whichflag == 0) {
+      if (pe->invoked_scalar != update->ntimestep)
+        error->all(FLERR,"Compute used in variable thermo keyword between runs "
+                   "is not current");
+    } else {
+      pe->compute_scalar();
+      pe->invoked_flag |= INVOKED_SCALAR;
+    }
+    if (!temperature)
+      error->all(FLERR,"Thermo keyword in variable requires "
+                 "thermo to use/init temp");
+    if (update->whichflag == 0) {
+      if (temperature->invoked_scalar != update->ntimestep)
+        error->all(FLERR,"Compute used in variable thermo keyword between runs "
+                   "is not current");
+    } else if (!(temperature->invoked_flag & INVOKED_SCALAR)) {
+      temperature->compute_scalar();
+      temperature->invoked_flag |= INVOKED_SCALAR;
+    }
+    if (!pressure)
+      error->all(FLERR,"Thermo keyword in variable requires "
+                 "thermo to use/init press");
+    if (update->whichflag == 0) {
+      if (pressure->invoked_scalar != update->ntimestep)
+        error->all(FLERR,"Compute used in variable thermo keyword between runs "
+                   "is not current");
+    } else if (!(pressure->invoked_flag & INVOKED_SCALAR)) {
+      pressure->compute_scalar();
+      pressure->invoked_flag |= INVOKED_SCALAR;
+    }
+    compute_enthalpy();
+
+  } else if (strcmp(word,"ecouple") == 0) compute_ecouple();
+  } else if (strcmp(word,"econserve") == 0) {
+    if (!pe)
+      error->all(FLERR,
+                 "Thermo keyword in variable requires thermo to use/init pe");
+    if (update->whichflag == 0) {
+      if (pe->invoked_scalar != update->ntimestep)
+        error->all(FLERR,"Compute used in variable thermo keyword between runs "
+                   "is not current");
+    } else {
+      pe->compute_scalar();
+      pe->invoked_flag |= INVOKED_SCALAR;
+    }
+    if (!temperature)
+      error->all(FLERR,"Thermo keyword in variable requires "
+                 "thermo to use/init temp");
+    if (update->whichflag == 0) {
+      if (temperature->invoked_scalar != update->ntimestep)
+        error->all(FLERR,"Compute used in variable thermo keyword between runs "
+                   "is not current");
+    } else if (!(temperature->invoked_flag & INVOKED_SCALAR)) {
+      temperature->compute_scalar();
+      temperature->invoked_flag |= INVOKED_SCALAR;
+    }
+    compute_econserve();
+
   } else if (strcmp(word,"vol") == 0) compute_vol();
   else if (strcmp(word,"density") == 0) compute_density();
   else if (strcmp(word,"lx") == 0) compute_lx();
@@ -1721,27 +1756,28 @@ void Thermo::compute_ke()
 void Thermo::compute_etotal()
 {
   compute_pe();
-  double ke = temperature->scalar;
+  double dvalue_pe = dvalue;
-  ke *= 0.5 * temperature->dof * force->boltz;
+  compute_ke();
-  if (normflag) ke /= natoms;
+  dvalue += dvalue_pe;
-  dvalue += ke;
 }
 
 /* ---------------------------------------------------------------------- */
 
-void Thermo::compute_enthalpy()
+void Thermo::compute_ecouple()
 {
-  compute_etotal();
+  dvalue = modify->ecouple();
-  double etmp = dvalue;
+}
 
-  compute_vol();
+/* ---------------------------------------------------------------------- */
-  double vtmp = dvalue;
-  if (normflag) vtmp /= natoms;
 
-  compute_press();
+void Thermo::compute_econserve()
-  double ptmp = dvalue;
+{
-
+  compute_pe();
-  dvalue = etmp + ptmp*vtmp/(force->nktv2p);
+  double dvalue_pe = dvalue;
+  compute_ke();
+  double dvalue_ke = dvalue;
+  compute_ecouple();
+  dvalue += dvalue_pe + dvalue_ke;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1867,6 +1903,24 @@ void Thermo::compute_etail()
   } else dvalue = 0.0;
 }
 
+/* ---------------------------------------------------------------------- */
+
+void Thermo::compute_enthalpy()
+{
+  compute_etotal();
+  double etmp = dvalue;
+
+  compute_vol();
+  double vtmp = dvalue;
+  if (normflag) vtmp /= natoms;
+
+  compute_press();
+  double ptmp = dvalue;
+
+  dvalue = etmp + ptmp*vtmp/(force->nktv2p);
+}
+
+
 /* ---------------------------------------------------------------------- */
 
 void Thermo::compute_vol()

src/thermo.h

@@ -143,7 +143,6 @@ class Thermo : protected Pointers {
   void compute_pe();
   void compute_ke();
   void compute_etotal();
-  void compute_enthalpy();
 
   void compute_evdwl();
   void compute_ecoul();
@@ -156,6 +155,10 @@ class Thermo : protected Pointers {
   void compute_elong();
   void compute_etail();
 
+  void compute_enthalpy();
+  void compute_ecouple();
+  void compute_econserve();
+
   void compute_vol();
   void compute_density();
   void compute_lx();