diff --git a/src/compute_pressure.cpp b/src/compute_pressure.cpp
index 7acf92ea2e..aadf96fdea 100644
--- a/src/compute_pressure.cpp
+++ b/src/compute_pressure.cpp
@@ -215,10 +215,13 @@ void ComputePressure::init()
       vptr[nvirial++] = force->dihedral->virial;
     if (improperflag && force->improper)
       vptr[nvirial++] = force->improper->virial;
-    if (fixflag)
-      for (int i = 0; i < modify->nfix; i++)
-        if (modify->fix[i]->thermo_virial)
+    if (fixflag) {
+      Fix **fix = modify->fix;
+      int nfix = modify->nfix;
+      for (int i = 0; i < nfix; i++)
+        if (fix[i]->virial_global_flag && fix[i]->thermo_virial)
           vptr[nvirial++] = modify->fix[i]->virial;
+    }
   }
 
   // flag Kspace contribution separately, since not summed across procs
diff --git a/src/compute_stress_atom.cpp b/src/compute_stress_atom.cpp
index f918f7f293..aefb58dbb4 100644
--- a/src/compute_stress_atom.cpp
+++ b/src/compute_stress_atom.cpp
@@ -219,8 +219,10 @@ void ComputeStressAtom::compute_peratom()
   //   and fix ave/spatial uses a per-atom stress from this compute as input
 
   if (fixflag) {
-    for (int ifix = 0; ifix < modify->nfix; ifix++)
-      if (modify->fix[ifix]->virial_flag) {
+    Fix **fix = modify->fix;
+    int nfix = modify->nfix;
+    for (int ifix = 0; ifix < nfix; ifix++)
+      if (fix[i]->virial_atom_flag && fix[ifix]->virial_flag) {
         double **vatom = modify->fix[ifix]->vatom;
         if (vatom)
           for (i = 0; i < nlocal; i++)
diff --git a/src/fix_addforce.cpp b/src/fix_addforce.cpp
index 8bd5bc69a3..22c7408996 100644
--- a/src/fix_addforce.cpp
+++ b/src/fix_addforce.cpp
@@ -36,8 +36,8 @@ enum{NONE,CONSTANT,EQUAL,ATOM};
 
 FixAddForce::FixAddForce(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg),
-  xstr(nullptr), ystr(nullptr), zstr(nullptr), estr(nullptr), idregion(nullptr), sforce(nullptr)
-
+  xstr(nullptr), ystr(nullptr), zstr(nullptr), estr(nullptr),
+  idregion(nullptr), sforce(nullptr)
 {
   if (narg < 6) error->all(FLERR,"Illegal fix addforce command");
 
@@ -48,9 +48,10 @@ FixAddForce::FixAddForce(LAMMPS *lmp, int narg, char **arg) :
   global_freq = 1;
   extscalar = 1;
   extvector = 1;
+  energy_global_flag = 1;
+  virial_global_flag = virial_atom_flag = 1;
   respa_level_support = 1;
   ilevel_respa = 0;
-  virial_flag = 1;
 
   xstr = ystr = zstr = nullptr;
 
@@ -138,7 +139,6 @@ int FixAddForce::setmask()
 
   int mask = 0;
   mask |= POST_FORCE;
-  mask |= THERMO_ENERGY;
   mask |= POST_FORCE_RESPA;
   mask |= MIN_POST_FORCE;
   return mask;
@@ -242,10 +242,9 @@ void FixAddForce::post_force(int vflag)
 
   if (update->ntimestep % nevery) return;
 
-  // energy and virial setup
+  // virial setup
 
-  if (vflag) v_setup(vflag);
-  else evflag = 0;
+  v_init(vflag);
 
   if (lmp->kokkos)
     atom->sync_modify(Host, (unsigned int) (F_MASK | MASK_MASK),
diff --git a/src/fix_external.cpp b/src/fix_external.cpp
index 81ab1ec36a..f3ad25d3b3 100644
--- a/src/fix_external.cpp
+++ b/src/fix_external.cpp
@@ -35,9 +35,10 @@ FixExternal::FixExternal(LAMMPS *lmp, int narg, char **arg) :
 
   scalar_flag = 1;
   global_freq = 1;
-  virial_flag = 1;
-  thermo_virial = 1;
   extscalar = 1;
+  energy_global_flag = energy_atom_flag = 1;
+  virial_global_flag = virial_atom_flag = 1;
+  thermo_energy = thermo_virial = 1;
 
   if (strcmp(arg[3],"pf/callback") == 0) {
     if (narg != 6) error->all(FLERR,"Illegal fix external command");
@@ -89,7 +90,6 @@ int FixExternal::setmask()
   if (mode == PF_CALLBACK || mode == PF_ARRAY) {
     mask |= PRE_REVERSE;
     mask |= POST_FORCE;
-    mask |= THERMO_ENERGY;
     mask |= MIN_POST_FORCE;
   }
   return mask;
@@ -195,7 +195,6 @@ void FixExternal::set_energy_global(double caller_energy)
 
 void FixExternal::set_virial_global(double *caller_virial)
 {
-  if (!evflag) return;
   if (!vflag_global) return;
 
   for (int i = 0; i < 6; i++)
@@ -223,7 +222,6 @@ void FixExternal::set_virial_peratom(double **caller_virial)
 {
   int i,j;
 
-  if (!evflag) return;
   if (!vflag_atom) return;
 
   int nlocal = atom->nlocal;
diff --git a/src/fix_gravity.cpp b/src/fix_gravity.cpp
index a96303742c..b4e894316a 100644
--- a/src/fix_gravity.cpp
+++ b/src/fix_gravity.cpp
@@ -46,6 +46,7 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
   scalar_flag = 1;
   global_freq = 1;
   extscalar = 1;
+  energy_global_flag = 1;
   respa_level_support = 1;
   ilevel_respa = 0;
 
@@ -62,7 +63,7 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
     mstyle = CONSTANT;
   }
 
-  int iarg=4;
+  int iarg = 4;
 
   if (strcmp(arg[4],"chute") == 0) {
     if (narg < 6) error->all(FLERR,"Illegal fix gravity command");
@@ -186,7 +187,6 @@ int FixGravity::setmask()
 {
   int mask = 0;
   mask |= POST_FORCE;
-  mask |= THERMO_ENERGY;
   mask |= POST_FORCE_RESPA;
   return mask;
 }
diff --git a/src/fix_indent.cpp b/src/fix_indent.cpp
index cc7a3c8877..b8f2d23312 100644
--- a/src/fix_indent.cpp
+++ b/src/fix_indent.cpp
@@ -47,6 +47,7 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) :
   scalar_flag = 1;
   vector_flag = 1;
   size_vector = 3;
+  energy_global_flag = 1;
   global_freq = 1;
   extscalar = 1;
   extvector = 1;
@@ -107,7 +108,6 @@ int FixIndent::setmask()
 {
   int mask = 0;
   mask |= POST_FORCE;
-  mask |= THERMO_ENERGY;
   mask |= POST_FORCE_RESPA;
   mask |= MIN_POST_FORCE;
   return mask;
diff --git a/src/fix_langevin.cpp b/src/fix_langevin.cpp
index 8bc2aa3976..f798252a50 100644
--- a/src/fix_langevin.cpp
+++ b/src/fix_langevin.cpp
@@ -39,7 +39,6 @@
 #include "error.h"
 #include "group.h"
 
-
 using namespace LAMMPS_NS;
 using namespace FixConst;
 
@@ -62,6 +61,7 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
   scalar_flag = 1;
   global_freq = 1;
   extscalar = 1;
+  ecouple_flag = 1;
   nevery = 1;
 
   if (strstr(arg[3],"v_") == arg[3]) {
@@ -159,7 +159,7 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
   id_temp = nullptr;
   temperature = nullptr;
 
-  energy = 0.0;
+  ecouple = 0.0;
 
   // flangevin is unallocated until first call to setup()
   // compute_scalar checks for this and returns 0.0
@@ -192,7 +192,6 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
       lv[i][2] = 0.0;
     }
   }
-
 }
 
 /* ---------------------------------------------------------------------- */
@@ -224,7 +223,6 @@ int FixLangevin::setmask()
   mask |= POST_FORCE;
   mask |= POST_FORCE_RESPA;
   mask |= END_OF_STEP;
-  mask |= THERMO_ENERGY;
   return mask;
 }
 
@@ -992,7 +990,7 @@ void FixLangevin::end_of_step()
       }
   }
 
-  energy += energy_onestep*update->dt;
+  ecouple += energy_onestep*update->dt;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1070,7 +1068,7 @@ double FixLangevin::compute_scalar()
         if (mask[i] & groupbit)
           energy_onestep += flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] +
                             flangevin[i][2]*v[i][2];
-      energy = 0.5*energy_onestep*update->dt;
+      ecouple = 0.5*energy_onestep*update->dt;
     } else {
       for (int i = 0; i < nlocal; i++)
         if (mask[i] & groupbit) {
@@ -1081,13 +1079,13 @@ double FixLangevin::compute_scalar()
           if (tbiasflag)
             temperature->restore_bias(i, lv[i]);
         }
-      energy = -0.5*energy_onestep*update->dt;
+      ecouple = -0.5*energy_onestep*update->dt;
     }
   }
 
   // convert midstep energy back to previous fullstep energy
 
-  double energy_me = energy - 0.5*energy_onestep*update->dt;
+  double energy_me = ecouple - 0.5*energy_onestep*update->dt;
 
   double energy_all;
   MPI_Allreduce(&energy_me,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp
index 674cda4529..ee2b74ddb9 100644
--- a/src/fix_nh.cpp
+++ b/src/fix_nh.cpp
@@ -66,6 +66,7 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) :
   global_freq = 1;
   extscalar = 1;
   extvector = 0;
+  ecouple_flag = 1;
 
   // default values
 
@@ -632,7 +633,6 @@ int FixNH::setmask()
   int mask = 0;
   mask |= INITIAL_INTEGRATE;
   mask |= FINAL_INTEGRATE;
-  mask |= THERMO_ENERGY;
   mask |= INITIAL_INTEGRATE_RESPA;
   mask |= FINAL_INTEGRATE_RESPA;
   if (pre_exchange_flag) mask |= PRE_EXCHANGE;
diff --git a/src/fix_restrain.cpp b/src/fix_restrain.cpp
index d1a05c2aca..82e91c8289 100644
--- a/src/fix_restrain.cpp
+++ b/src/fix_restrain.cpp
@@ -30,7 +30,6 @@
 #include "memory.h"
 #include "error.h"
 
-
 using namespace LAMMPS_NS;
 using namespace FixConst;
 using namespace MathConst;
@@ -57,6 +56,7 @@ FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) :
   vector_flag = 1;
   size_vector = 3;
   extvector = 1;
+  energy_global_flag = 1;
   respa_level_support = 1;
   ilevel_respa = 0;
 
@@ -176,7 +176,6 @@ int FixRestrain::setmask()
 {
   int mask = 0;
   mask |= POST_FORCE;
-  mask |= THERMO_ENERGY;
   mask |= POST_FORCE_RESPA;
   mask |= MIN_POST_FORCE;
   return mask;
diff --git a/src/fix_spring.cpp b/src/fix_spring.cpp
index d3f89ec97a..4bbaf594bd 100644
--- a/src/fix_spring.cpp
+++ b/src/fix_spring.cpp
@@ -47,6 +47,7 @@ FixSpring::FixSpring(LAMMPS *lmp, int narg, char **arg) :
   global_freq = 1;
   extscalar = 1;
   extvector = 1;
+  energy_global_flag = 1;
   dynamic_group_allow = 1;
   respa_level_support = 1;
   ilevel_respa = 0;
@@ -108,7 +109,6 @@ int FixSpring::setmask()
 {
   int mask = 0;
   mask |= POST_FORCE;
-  mask |= THERMO_ENERGY;
   mask |= POST_FORCE_RESPA;
   mask |= MIN_POST_FORCE;
   return mask;
diff --git a/src/fix_spring_chunk.cpp b/src/fix_spring_chunk.cpp
index 130998dcd5..9ae4ee4632 100644
--- a/src/fix_spring_chunk.cpp
+++ b/src/fix_spring_chunk.cpp
@@ -43,6 +43,7 @@ FixSpringChunk::FixSpringChunk(LAMMPS *lmp, int narg, char **arg) :
   scalar_flag = 1;
   global_freq = 1;
   extscalar = 1;
+  energy_global_flag = 1;
   respa_level_support = 1;
   ilevel_respa = 0;
 
@@ -86,7 +87,6 @@ int FixSpringChunk::setmask()
 {
   int mask = 0;
   mask |= POST_FORCE;
-  mask |= THERMO_ENERGY;
   mask |= POST_FORCE_RESPA;
   mask |= MIN_POST_FORCE;
   return mask;
diff --git a/src/fix_spring_self.cpp b/src/fix_spring_self.cpp
index 50d09acf50..b736271912 100644
--- a/src/fix_spring_self.cpp
+++ b/src/fix_spring_self.cpp
@@ -42,6 +42,7 @@ FixSpringSelf::FixSpringSelf(LAMMPS *lmp, int narg, char **arg) :
   scalar_flag = 1;
   global_freq = 1;
   extscalar = 1;
+  energy_global_flag = 1;
   respa_level_support = 1;
 
   k = utils::numeric(FLERR,arg[3],false,lmp);
@@ -110,7 +111,6 @@ int FixSpringSelf::setmask()
 {
   int mask = 0;
   mask |= POST_FORCE;
-  mask |= THERMO_ENERGY;
   mask |= POST_FORCE_RESPA;
   mask |= MIN_POST_FORCE;
   return mask;
diff --git a/src/fix_temp_berendsen.cpp b/src/fix_temp_berendsen.cpp
index afa41c50e4..61cec81e8f 100644
--- a/src/fix_temp_berendsen.cpp
+++ b/src/fix_temp_berendsen.cpp
@@ -26,7 +26,6 @@
 #include "compute.h"
 #include "error.h"
 
-
 using namespace LAMMPS_NS;
 using namespace FixConst;
 
@@ -45,10 +44,11 @@ FixTempBerendsen::FixTempBerendsen(LAMMPS *lmp, int narg, char **arg) :
 
   restart_global = 1;
   dynamic_group_allow = 1;
-  nevery = 1;
   scalar_flag = 1;
-  global_freq = nevery;
   extscalar = 1;
+  ecouple_flag = 1;
+  nevery = 1;
+  global_freq = nevery;
 
   tstr = nullptr;
   if (strstr(arg[3],"v_") == arg[3]) {
@@ -81,7 +81,7 @@ FixTempBerendsen::FixTempBerendsen(LAMMPS *lmp, int narg, char **arg) :
   modify->add_compute(cmd);
   tflag = 1;
 
-  energy = 0;
+  ecouple = 0;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -102,7 +102,6 @@ int FixTempBerendsen::setmask()
 {
   int mask = 0;
   mask |= END_OF_STEP;
-  mask |= THERMO_ENERGY;
   return mask;
 }
 
@@ -171,7 +170,7 @@ void FixTempBerendsen::end_of_step()
 
   double lamda = sqrt(1.0 + update->dt/t_period*(t_target/t_current - 1.0));
   double efactor = 0.5 * force->boltz * tdof;
-  energy += t_current * (1.0-lamda*lamda) * efactor;
+  ecouple += t_current * (1.0-lamda*lamda) * efactor;
 
   double **v = atom->v;
   int *mask = atom->mask;
@@ -239,7 +238,7 @@ void FixTempBerendsen::reset_target(double t_new)
 
 double FixTempBerendsen::compute_scalar()
 {
-  return energy;
+  return ecouple;
 }
 
 /* ----------------------------------------------------------------------
@@ -250,7 +249,7 @@ void FixTempBerendsen::write_restart(FILE *fp)
 {
   int n = 0;
   double list[1];
-  list[n++] = energy;
+  list[n++] = ecouple;
 
   if (comm->me == 0) {
     int size = n * sizeof(double);
@@ -267,7 +266,7 @@ void FixTempBerendsen::restart(char *buf)
 {
   double *list = (double *) buf;
 
-  energy = list[0];
+  ecouple = list[0];
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/fix_temp_csld.cpp b/src/fix_temp_csld.cpp
index 3b522c185f..da6447b694 100644
--- a/src/fix_temp_csld.cpp
+++ b/src/fix_temp_csld.cpp
@@ -32,7 +32,6 @@
 #include "random_mars.h"
 #include "error.h"
 
-
 using namespace LAMMPS_NS;
 using namespace FixConst;
 
@@ -52,6 +51,7 @@ FixTempCSLD::FixTempCSLD(LAMMPS *lmp, int narg, char **arg) :
   restart_global = 1;
   nevery = 1;
   scalar_flag = 1;
+  ecouple_flag = 1;
   global_freq = nevery;
   dynamic_group_allow = 1;
   extscalar = 1;
@@ -92,7 +92,7 @@ FixTempCSLD::FixTempCSLD(LAMMPS *lmp, int narg, char **arg) :
 
   vhold = nullptr;
   nmax = -1;
-  energy = 0.0;
+  ecouple = 0.0;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -118,7 +118,6 @@ int FixTempCSLD::setmask()
 {
   int mask = 0;
   mask |= END_OF_STEP;
-  mask |= THERMO_ENERGY;
   return mask;
 }
 
@@ -251,7 +250,7 @@ void FixTempCSLD::end_of_step()
   // tally the kinetic energy transferred between heat bath and system
 
   t_current = temperature->compute_scalar();
-  energy +=  ekin_old - t_current * 0.5 * temperature->dof * force->boltz;
+  ecouple +=  ekin_old - t_current * 0.5 * temperature->dof * force->boltz;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -295,10 +294,9 @@ void FixTempCSLD::reset_target(double t_new)
 
 double FixTempCSLD::compute_scalar()
 {
-  return energy;
+  return ecouple;
 }
 
-
 /* ----------------------------------------------------------------------
    pack entire state of Fix into one write
 ------------------------------------------------------------------------- */
@@ -306,11 +304,11 @@ double FixTempCSLD::compute_scalar()
 void FixTempCSLD::write_restart(FILE *fp)
 {
   const int PRNGSIZE = 98+2+3;
-  int nsize = PRNGSIZE*comm->nprocs+2; // pRNG state per proc + nprocs + energy
+  int nsize = PRNGSIZE*comm->nprocs+2; // pRNG state per proc + nprocs + ecouple
   double *list = nullptr;
   if (comm->me == 0) {
     list = new double[nsize];
-    list[0] = energy;
+    list[0] = ecouple;
     list[1] = comm->nprocs;
   }
   double state[PRNGSIZE];
@@ -333,7 +331,7 @@ void FixTempCSLD::restart(char *buf)
 {
   double *list = (double *) buf;
 
-  energy = list[0];
+  ecouple = list[0];
   int nprocs = (int) list[1];
   if (nprocs != comm->nprocs) {
     if (comm->me == 0)
diff --git a/src/fix_temp_csvr.cpp b/src/fix_temp_csvr.cpp
index eedf24c387..8779e7fb5e 100644
--- a/src/fix_temp_csvr.cpp
+++ b/src/fix_temp_csvr.cpp
@@ -33,89 +33,12 @@
 #include "random_mars.h"
 #include "error.h"
 
-
 using namespace LAMMPS_NS;
 using namespace FixConst;
 
 enum{NOBIAS,BIAS};
 enum{CONSTANT,EQUAL};
 
-double FixTempCSVR::gamdev(const int ia)
-{
-  int j;
-  double am,e,s,v1,v2,x,y;
-
-  if (ia < 1) return 0.0;
-  if (ia < 6) {
-    x=1.0;
-    for (j=1; j<=ia; j++)
-      x *= random->uniform();
-
-    // make certain, that -log() doesn't overflow.
-    if (x < 2.2250759805e-308)
-      x = 708.4;
-    else
-      x = -log(x);
-  } else {
-  restart:
-    do {
-      do {
-        do {
-          v1 = random->uniform();
-          v2 = 2.0*random->uniform() - 1.0;
-        } while (v1*v1 + v2*v2 > 1.0);
-
-        y=v2/v1;
-        am=ia-1;
-        s=sqrt(2.0*am+1.0);
-        x=s*y+am;
-      } while (x <= 0.0);
-
-      if (am*log(x/am)-s*y < -700 || v1<0.00001) {
-        goto restart;
-      }
-
-      e=(1.0+y*y)*exp(am*log(x/am)-s*y);
-    } while (random->uniform() > e);
-  }
-  return x;
-}
-
-/* -------------------------------------------------------------------
-  returns the sum of n independent gaussian noises squared
-  (i.e. equivalent to summing the square of the return values of nn
-   calls to gasdev)
----------------------------------------------------------------------- */
-double FixTempCSVR::sumnoises(int nn) {
-  if (nn == 0) {
-    return 0.0;
-  } else if (nn == 1) {
-    const double rr = random->gaussian();
-    return rr*rr;
-  } else if (nn % 2 == 0) {
-    return 2.0 * gamdev(nn / 2);
-  } else {
-    const double rr = random->gaussian();
-    return  2.0 * gamdev((nn-1) / 2) + rr*rr;
-  }
-}
-
-/* -------------------------------------------------------------------
-  returns the scaling factor for velocities to thermalize
-  the system so it samples the canonical ensemble
----------------------------------------------------------------------- */
-
-double FixTempCSVR::resamplekin(double ekin_old, double ekin_new) {
-  const double tdof = temperature->dof;
-  const double c1 = exp(-update->dt/t_period);
-  const double c2 = (1.0-c1)*ekin_new/ekin_old/tdof;
-  const double r1 = random->gaussian();
-  const double r2 = sumnoises(tdof - 1);
-
-  const double scale = c1 + c2*(r1*r1+r2) + 2.0*r1*sqrt(c1*c2);
-  return sqrt(scale);
-}
-
 /* ---------------------------------------------------------------------- */
 
 FixTempCSVR::FixTempCSVR(LAMMPS *lmp, int narg, char **arg) :
@@ -129,6 +52,7 @@ FixTempCSVR::FixTempCSVR(LAMMPS *lmp, int narg, char **arg) :
   restart_global = 1;
   nevery = 1;
   scalar_flag = 1;
+  ecouple_flag = 1;
   global_freq = nevery;
   dynamic_group_allow = 1;
   extscalar = 1;
@@ -168,7 +92,7 @@ FixTempCSVR::FixTempCSVR(LAMMPS *lmp, int narg, char **arg) :
   tflag = 1;
 
   nmax = -1;
-  energy = 0.0;
+  ecouple = 0.0;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -192,7 +116,6 @@ int FixTempCSVR::setmask()
 {
   int mask = 0;
   mask |= END_OF_STEP;
-  mask |= THERMO_ENERGY;
   return mask;
 }
 
@@ -200,7 +123,6 @@ int FixTempCSVR::setmask()
 
 void FixTempCSVR::init()
 {
-
   // check variable
 
   if (tstr) {
@@ -224,7 +146,6 @@ void FixTempCSVR::init()
 
 void FixTempCSVR::end_of_step()
 {
-
   // set current t_target
   // if variable temp, evaluate variable, wrap with clear/add
 
@@ -285,7 +206,7 @@ void FixTempCSVR::end_of_step()
 
   // tally the kinetic energy transferred between heat bath and system
 
-  energy += ekin_old * (1.0 - lamda*lamda);
+  ecouple += ekin_old * (1.0 - lamda*lamda);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -320,6 +241,85 @@ int FixTempCSVR::modify_param(int narg, char **arg)
 
 /* ---------------------------------------------------------------------- */
 
+double FixTempCSVR::gamdev(const int ia)
+{
+  int j;
+  double am,e,s,v1,v2,x,y;
+
+  if (ia < 1) return 0.0;
+  if (ia < 6) {
+    x=1.0;
+    for (j=1; j<=ia; j++)
+      x *= random->uniform();
+
+    // make certain, that -log() doesn't overflow.
+    if (x < 2.2250759805e-308)
+      x = 708.4;
+    else
+      x = -log(x);
+  } else {
+  restart:
+    do {
+      do {
+        do {
+          v1 = random->uniform();
+          v2 = 2.0*random->uniform() - 1.0;
+        } while (v1*v1 + v2*v2 > 1.0);
+
+        y=v2/v1;
+        am=ia-1;
+        s=sqrt(2.0*am+1.0);
+        x=s*y+am;
+      } while (x <= 0.0);
+
+      if (am*log(x/am)-s*y < -700 || v1<0.00001) {
+        goto restart;
+      }
+
+      e=(1.0+y*y)*exp(am*log(x/am)-s*y);
+    } while (random->uniform() > e);
+  }
+  return x;
+}
+
+/* -------------------------------------------------------------------
+  returns the sum of n independent gaussian noises squared
+  (i.e. equivalent to summing the square of the return values of nn
+   calls to gasdev)
+---------------------------------------------------------------------- */
+
+double FixTempCSVR::sumnoises(int nn) {
+  if (nn == 0) {
+    return 0.0;
+  } else if (nn == 1) {
+    const double rr = random->gaussian();
+    return rr*rr;
+  } else if (nn % 2 == 0) {
+    return 2.0 * gamdev(nn / 2);
+  } else {
+    const double rr = random->gaussian();
+    return  2.0 * gamdev((nn-1) / 2) + rr*rr;
+  }
+}
+
+/* -------------------------------------------------------------------
+  returns the scaling factor for velocities to thermalize
+  the system so it samples the canonical ensemble
+---------------------------------------------------------------------- */
+
+double FixTempCSVR::resamplekin(double ekin_old, double ekin_new) {
+  const double tdof = temperature->dof;
+  const double c1 = exp(-update->dt/t_period);
+  const double c2 = (1.0-c1)*ekin_new/ekin_old/tdof;
+  const double r1 = random->gaussian();
+  const double r2 = sumnoises(tdof - 1);
+
+  const double scale = c1 + c2*(r1*r1+r2) + 2.0*r1*sqrt(c1*c2);
+  return sqrt(scale);
+}
+
+/* ---------------------------------------------------------------------- */
+
 void FixTempCSVR::reset_target(double t_new)
 {
   t_target = t_start = t_stop = t_new;
@@ -329,7 +329,7 @@ void FixTempCSVR::reset_target(double t_new)
 
 double FixTempCSVR::compute_scalar()
 {
-  return energy;
+  return ecouple;
 }
 
 /* ----------------------------------------------------------------------
@@ -339,11 +339,11 @@ double FixTempCSVR::compute_scalar()
 void FixTempCSVR::write_restart(FILE *fp)
 {
   const int PRNGSIZE = 98+2+3;
-  int nsize = PRNGSIZE*comm->nprocs+2; // pRNG state per proc + nprocs + energy
+  int nsize = PRNGSIZE*comm->nprocs+2; // pRNG state per proc + nprocs + ecouple
   double *list = nullptr;
   if (comm->me == 0) {
     list = new double[nsize];
-    list[0] = energy;
+    list[0] = ecouple;
     list[1] = comm->nprocs;
   }
   double state[PRNGSIZE];
@@ -366,7 +366,7 @@ void FixTempCSVR::restart(char *buf)
 {
   double *list = (double *) buf;
 
-  energy = list[0];
+  ecouple = list[0];
   int nprocs = (int) list[1];
   if (nprocs != comm->nprocs) {
     if (comm->me == 0)
diff --git a/src/fix_temp_rescale.cpp b/src/fix_temp_rescale.cpp
index 5e6638a56c..80470e3ccb 100644
--- a/src/fix_temp_rescale.cpp
+++ b/src/fix_temp_rescale.cpp
@@ -48,6 +48,7 @@ FixTempRescale::FixTempRescale(LAMMPS *lmp, int narg, char **arg) :
   scalar_flag = 1;
   global_freq = nevery;
   extscalar = 1;
+  ecouple_flag = 1;
   dynamic_group_allow = 1;
 
   tstr = nullptr;
@@ -77,7 +78,7 @@ FixTempRescale::FixTempRescale(LAMMPS *lmp, int narg, char **arg) :
   modify->add_compute(cmd);
   tflag = 1;
 
-  energy = 0.0;
+  ecouple = 0.0;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -98,7 +99,6 @@ int FixTempRescale::setmask()
 {
   int mask = 0;
   mask |= END_OF_STEP;
-  mask |= THERMO_ENERGY;
   return mask;
 }
 
@@ -171,7 +171,7 @@ void FixTempRescale::end_of_step()
     int *mask = atom->mask;
     int nlocal = atom->nlocal;
 
-    energy += (t_current-t_target) * efactor;
+    ecouple += (t_current-t_target) * efactor;
 
     if (which == NOBIAS) {
       for (int i = 0; i < nlocal; i++) {
@@ -236,7 +236,7 @@ void FixTempRescale::reset_target(double t_new)
 
 double FixTempRescale::compute_scalar()
 {
-  return energy;
+  return ecouple;
 }
 
 /* ----------------------------------------------------------------------
@@ -247,7 +247,7 @@ void FixTempRescale::write_restart(FILE *fp)
 {
   int n = 0;
   double list[1];
-  list[n++] = energy;
+  list[n++] = ecouple;
 
   if (comm->me == 0) {
     int size = n * sizeof(double);
@@ -265,7 +265,7 @@ void FixTempRescale::restart(char *buf)
   int n = 0;
   double *list = (double *) buf;
 
-  energy = list[n++];
+  ecouple = list[n++];
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/fix_wall.cpp b/src/fix_wall.cpp
index 31aae965d4..5740045a63 100644
--- a/src/fix_wall.cpp
+++ b/src/fix_wall.cpp
@@ -41,9 +41,10 @@ FixWall::FixWall(LAMMPS *lmp, int narg, char **arg) :
   global_freq = 1;
   extscalar = 1;
   extvector = 1;
+  energy_global_flag = 1;
+  virial_global_flag = virial_atom_flag = 1;
   respa_level_support = 1;
   ilevel_respa = 0;
-  virial_flag = 1;
 
   // parse args
 
@@ -233,7 +234,6 @@ int FixWall::setmask()
   if (fldflag) mask |= PRE_FORCE;
   else mask |= POST_FORCE;
 
-  mask |= THERMO_ENERGY;
   mask |= POST_FORCE_RESPA;
   mask |= MIN_POST_FORCE;
   return mask;
@@ -310,12 +310,12 @@ void FixWall::pre_force(int vflag)
 
 void FixWall::post_force(int vflag)
 {
+  // virial setup
 
-  // energy and virial setup
+  v_init(vflag);
 
-  eflag = 0;
-  if (vflag) v_setup(vflag);
-  else evflag = 0;
+  // energy intialize
+  
   for (int m = 0; m <= nwall; m++) ewall[m] = 0.0;
 
   // coord = current position of wall
diff --git a/src/fix_wall_harmonic.cpp b/src/fix_wall_harmonic.cpp
index 864125d93f..ffd950262b 100644
--- a/src/fix_wall_harmonic.cpp
+++ b/src/fix_wall_harmonic.cpp
@@ -67,7 +67,7 @@ void FixWallHarmonic::wall_particle(int m, int which, double coord)
       if (evflag) {
         if (side < 0) vn = -fwall*delta;
         else vn = fwall*delta;
-        v_tally(dim, i, vn);
+        v_tally(dim,i,vn);
       }
     }
 
diff --git a/src/fix_wall_region.cpp b/src/fix_wall_region.cpp
index 3107f876de..2399dbfe94 100644
--- a/src/fix_wall_region.cpp
+++ b/src/fix_wall_region.cpp
@@ -44,9 +44,10 @@ FixWallRegion::FixWallRegion(LAMMPS *lmp, int narg, char **arg) :
   global_freq = 1;
   extscalar = 1;
   extvector = 1;
+  energy_global_flag = 1;
+  virial_global_flag = virial_atom_flag = 1;
   respa_level_support = 1;
   ilevel_respa = 0;
-  virial_flag = 1;
 
   // parse args
 
@@ -104,7 +105,6 @@ int FixWallRegion::setmask()
 {
   int mask = 0;
   mask |= POST_FORCE;
-  mask |= THERMO_ENERGY;
   mask |= POST_FORCE_RESPA;
   mask |= MIN_POST_FORCE;
   return mask;
@@ -234,11 +234,9 @@ void FixWallRegion::post_force(int vflag)
 
   int onflag = 0;
 
-  // energy and virial setup
+  // virial setup
 
-  eflag = 0;
-  if (vflag) v_setup(vflag);
-  else evflag = 0;
+  v_init(vflag);
 
   // region->match() insures particle is in region or on surface, else error
   // if returned contact dist r = 0, is on surface, also an error