From b1b7f7a248c8c6c13917f98ae750b6e9971d0180 Mon Sep 17 00:00:00 2001
From: Stephen Sanderson <stephen.sanderson@uq.edu.au>
Date: Wed, 17 Mar 2021 11:47:32 +1000
Subject: [PATCH] Fixed treatment of DoF when streaming velocity not subtracted
 in some dimensions

---
 src/compute_temp_profile.cpp | 13 ++++++++++---
 src/compute_temp_profile.h   |  1 +
 2 files changed, 11 insertions(+), 3 deletions(-)

diff --git a/src/compute_temp_profile.cpp b/src/compute_temp_profile.cpp
index 52f2f7de6f..9790faf9b5 100644
--- a/src/compute_temp_profile.cpp
+++ b/src/compute_temp_profile.cpp
@@ -118,6 +118,9 @@ ComputeTempProfile::ComputeTempProfile(LAMMPS *lmp, int narg, char **arg) :
   nbins = nbinx*nbiny*nbinz;
   if (nbins <= 0) error->all(FLERR,"Illegal compute temp/profile command");
 
+  nconstraints = (xflag==0 ? 0 : 1) + (yflag==0 ? 0 : 1) + (zflag==0 ? 0 : 1);
+  nconstraints *= nbins;
+
   memory->create(vbin,nbins,ncount,"temp/profile:vbin");
   memory->create(binave,nbins,ncount,"temp/profile:binave");
 
@@ -196,9 +199,10 @@ void ComputeTempProfile::dof_compute()
   natoms_temp = group->count(igroup);
   dof = domain->dimension * natoms_temp;
 
-  // subtract additional d*Nbins DOF, as in Evans and Morriss paper
+  // subtract additional Nbins DOF for each adjusted direction,
+  // as in Evans and Morriss paper
 
-  dof -= extra_dof + fix_dof + domain->dimension*nbins;
+  dof -= extra_dof + fix_dof + nconstraints;
   if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
   else tfactor = 0.0;
 }
@@ -342,7 +346,7 @@ void ComputeTempProfile::compute_array()
   double dofbin, tfactorbin;
   for (i = 0; i < nbins; i++) {
     if (array[i][0] > 0.0) {
-      dofbin = nper*array[i][0] - domain->dimension;
+      dofbin = nper*array[i][0] - nconstraints/nbins;
       if (dofbin > 0) tfactorbin = force->mvv2e / (dofbin * force->boltz);
       else tfactorbin = 0.0;
       array[i][1] = tfactorbin*tbinall[i];
@@ -582,6 +586,9 @@ void ComputeTempProfile::bin_assign()
 
 void ComputeTempProfile::reset_extra_dof() {
   extra_dof = 0.0;
+  if (xflag == 0) extra_dof++;
+  if (yflag == 0) extra_dof++;
+  if (zflag == 0 && domain->dimension == 3) extra_dof++;
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/compute_temp_profile.h b/src/compute_temp_profile.h
index 92e76f9095..5e62d923b7 100644
--- a/src/compute_temp_profile.h
+++ b/src/compute_temp_profile.h
@@ -48,6 +48,7 @@ class ComputeTempProfile : public Compute {
   int nbinx,nbiny,nbinz,nbins;
   int ivx,ivy,ivz;
   double tfactor;
+  double nconstraints;
 
   int box_change,triclinic;
   int *periodicity;