From e37ee02eedbfc5ee117caa2e8217f11737270dbf Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Tue, 24 Jul 2018 11:19:00 -0400
Subject: [PATCH] port nh fixes to USER-BOCS package

---
 src/USER-BOCS/fix_bocs.cpp | 45 +++++++++++++++++++++++++-------------
 1 file changed, 30 insertions(+), 15 deletions(-)

diff --git a/src/USER-BOCS/fix_bocs.cpp b/src/USER-BOCS/fix_bocs.cpp
index 37e128f556..7fb8a27110 100644
--- a/src/USER-BOCS/fix_bocs.cpp
+++ b/src/USER-BOCS/fix_bocs.cpp
@@ -846,7 +846,7 @@ void FixBocs::setup(int vflag)
 
   if (pstat_flag) {
     double kt = boltz * t_target;
-    double nkt = atom->natoms * kt;
+    double nkt = (atom->natoms + 1) * kt;
 
     for (int i = 0; i < 3; i++)
       if (p_flag[i])
@@ -1508,7 +1508,7 @@ double FixBocs::compute_scalar()
   double volume;
   double energy;
   double kt = boltz * t_target;
-  double lkt_press = kt;
+  double lkt_press = 0.0;
   int ich;
   if (dimension == 3) volume = domain->xprd * domain->yprd * domain->zprd;
   else volume = domain->xprd * domain->yprd;
@@ -1539,15 +1539,21 @@ double FixBocs::compute_scalar()
   //       sum is over barostatted dimensions
 
   if (pstat_flag) {
-    for (i = 0; i < 3; i++)
-      if (p_flag[i])
+    for (i = 0; i < 3; i++) {
+      if (p_flag[i]) {
         energy += 0.5*omega_dot[i]*omega_dot[i]*omega_mass[i] +
           p_hydro*(volume-vol0) / (pdim*nktv2p);
+        lkt_press += kt;
+      }
+    }
 
     if (pstyle == TRICLINIC) {
-      for (i = 3; i < 6; i++)
-        if (p_flag[i])
+      for (i = 3; i < 6; i++) {
+        if (p_flag[i]) {
           energy += 0.5*omega_dot[i]*omega_dot[i]*omega_mass[i];
+          lkt_press += kt;
+        }
+      }
     }
 
     // extra contributions from thermostat chain for barostat
@@ -1880,15 +1886,14 @@ void FixBocs::nhc_temp_integrate()
 
 void FixBocs::nhc_press_integrate()
 {
-  int ich,i;
+  int ich,i,pdof;
   double expfac,factor_etap,kecurrent;
   double kt = boltz * t_target;
-  double lkt_press = kt;
 
   // Update masses, to preserve initial freq, if flag set
 
   if (omega_mass_flag) {
-    double nkt = atom->natoms * kt;
+    double nkt = (atom->natoms + 1) * kt;
     for (int i = 0; i < 3; i++)
       if (p_flag[i])
         omega_mass[i] = nkt/(p_freq[i]*p_freq[i]);
@@ -1912,14 +1917,24 @@ void FixBocs::nhc_press_integrate()
   }
 
   kecurrent = 0.0;
-  for (i = 0; i < 3; i++)
-    if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
-
-  if (pstyle == TRICLINIC) {
-    for (i = 3; i < 6; i++)
-      if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
+  pdof = 0;
+  for (i = 0; i < 3; i++) {
+    if (p_flag[i]) {
+      kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
+      pdof++;
+    }
   }
 
+  if (pstyle == TRICLINIC) {
+    for (i = 3; i < 6; i++) {
+      if (p_flag[i]) {
+        kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
+        pdof++;
+      }
+    }
+  }
+
+  double lkt_press = pdof * kt;
   etap_dotdot[0] = (kecurrent - lkt_press)/etap_mass[0];
 
   double ncfac = 1.0/nc_pchain;