From 6f33ce9dba55ff55254c19bc1a4606452d194dfc Mon Sep 17 00:00:00 2001
From: Aidan Thompson <athomps@sandia.gov>
Date: Wed, 12 Dec 2018 17:31:42 -0700
Subject: [PATCH] Made two changes: 	-recomputed up-to-date pressure tensor
 (fixes energy conservation problem with aniso) 	-changed ndof for iso
 (fixes volume fluctuation problem with iso)

---
 src/compute_pressure.cpp | 54 ++++++++++++++++++++++++++++++++++++++++
 src/compute_pressure.h   |  1 +
 src/fix_nh.cpp           | 19 ++++++++++----
 3 files changed, 69 insertions(+), 5 deletions(-)

diff --git a/src/compute_pressure.cpp b/src/compute_pressure.cpp
index dde02a5aed..33d5a8f349 100644
--- a/src/compute_pressure.cpp
+++ b/src/compute_pressure.cpp
@@ -261,6 +261,60 @@ void ComputePressure::compute_vector()
   }
 }
 
+/* ----------------------------------------------------------------------
+   compute pressure tensor
+   assume KE tensor has already been computed
+------------------------------------------------------------------------- */
+
+void ComputePressure::compute_vector_ke_scalar()
+{
+  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 t;
+  if (keflag) {
+    if (temperature->invoked_scalar != update->ntimestep)
+      t = temperature->compute_scalar();
+    else t = temperature->scalar;
+  }
+
+  if (dimension == 3) {
+    inv_volume = 1.0 / (domain->xprd * domain->yprd * domain->zprd);
+    virial_compute(6,3);
+    if (keflag) {
+      double kescalar = temperature->dof * boltz * t / 3.0;
+      for (int i = 0; i < 3; i++)
+        vector[i] = (kescalar + virial[i]) * inv_volume * nktv2p;
+      for (int i = 3; i < 6; i++)
+        vector[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) {
+      double kescalar = temperature->dof * boltz * t / 2.0;
+      vector[0] = (kescalar + virial[0]) * inv_volume * nktv2p;
+      vector[1] = (kescalar + virial[1]) * inv_volume * nktv2p;
+      vector[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;
+    }
+  }
+}
+
 /* ---------------------------------------------------------------------- */
 
 void ComputePressure::virial_compute(int n, int ndiag)
diff --git a/src/compute_pressure.h b/src/compute_pressure.h
index a59a64e634..508aa45fe5 100644
--- a/src/compute_pressure.h
+++ b/src/compute_pressure.h
@@ -31,6 +31,7 @@ class ComputePressure : public Compute {
   virtual void init();
   virtual double compute_scalar();
   virtual void compute_vector();
+  void compute_vector_ke_scalar();
   void reset_extra_compute_fix(const char *);
 
  protected:
diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp
index a1a562f2bb..b65554dbe1 100644
--- a/src/fix_nh.cpp
+++ b/src/fix_nh.cpp
@@ -29,6 +29,7 @@
 #include "modify.h"
 #include "fix_deform.h"
 #include "compute.h"
+#include "compute_pressure.h"
 #include "kspace.h"
 #include "update.h"
 #include "respa.h"
@@ -777,7 +778,7 @@ void FixNH::setup(int /*vflag*/)
 
   if (pstat_flag) {
     if (pstyle == ISO) pressure->compute_scalar();
-    else pressure->compute_vector();
+    else ((ComputePressure *)pressure)->compute_vector_ke_scalar();
     couple();
     pressure->addstep(update->ntimestep+1);
   }
@@ -850,7 +851,7 @@ void FixNH::initial_integrate(int /*vflag*/)
       pressure->compute_scalar();
     } else {
       temperature->compute_vector();
-      pressure->compute_vector();
+      ((ComputePressure *)pressure)->compute_vector_ke_scalar();
     }
     couple();
     pressure->addstep(update->ntimestep+1);
@@ -904,9 +905,16 @@ void FixNH::final_integrate()
   t_current = temperature->compute_scalar();
   tdof = temperature->dof;
 
+  // need to recompute pressure to account for change in KE
+  // t_current is up-to-date, but compute_temperature is not
+  // compute appropriately coupled elements of mvv_current
+
   if (pstat_flag) {
     if (pstyle == ISO) pressure->compute_scalar();
-    else pressure->compute_vector();
+    else {
+      temperature->compute_vector();
+      ((ComputePressure *)pressure)->compute_vector_ke_scalar();
+    }
     couple();
     pressure->addstep(update->ntimestep+1);
   }
@@ -957,7 +965,7 @@ void FixNH::initial_integrate_respa(int /*vflag*/, int ilevel, int /*iloop*/)
         pressure->compute_scalar();
       } else {
         temperature->compute_vector();
-        pressure->compute_vector();
+        ((ComputePressure *)pressure)->compute_vector_ke_scalar();
       }
       couple();
       pressure->addstep(update->ntimestep+1);
@@ -1871,7 +1879,8 @@ void FixNH::nhc_press_integrate()
       }
   }
 
-  lkt_press = pdof * kt;
+  if (pstyle == ISO) lkt_press = kt;
+  else lkt_press = pdof * kt;
   etap_dotdot[0] = (kecurrent - lkt_press)/etap_mass[0];
 
   double ncfac = 1.0/nc_pchain;