From d1356c1d7633f639ed0bd58a564ccbd2453d60eb Mon Sep 17 00:00:00 2001
From: Aidan Thompson <athomps@sandia.gov>
Date: Tue, 3 Apr 2018 15:07:04 -0600
Subject: [PATCH] Made ndof dynamic for temperature fixes and computes

---
 examples/gcmc/in.gcmc.co2 |  7 ++++++-
 examples/gcmc/in.gcmc.h2o | 15 +++++++++------
 src/MC/fix_gcmc.cpp       |  8 ++++++++
 3 files changed, 23 insertions(+), 7 deletions(-)

diff --git a/examples/gcmc/in.gcmc.co2 b/examples/gcmc/in.gcmc.co2
index bb6916fc48..b5e11d212d 100644
--- a/examples/gcmc/in.gcmc.co2
+++ b/examples/gcmc/in.gcmc.co2
@@ -59,7 +59,9 @@ timestep        1.0
 # rigid constraints with thermostat 
 
 fix             myrigidnvt co2 rigid/nvt/small molecule temp ${temp} ${temp} 100 mol co2mol
-fix_modify	myrigidnvt dynamic/dof no
+
+# dynamically update  fix rigid/nvt/small temperature ndof
+fix_modify	myrigidnvt dynamic/dof yes
 
 # gcmc
 
@@ -82,7 +84,10 @@ variable	tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1)
 variable	iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
 variable	dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)
 variable	racc equal f_mygcmc[8]/(f_mygcmc[7]+0.1)
+
+# dynamically update default temperature ndof
 compute_modify  thermo_temp dynamic/dof yes
+
 thermo_style    custom step temp press pe ke density atoms v_iacc v_dacc v_tacc v_racc v_nC v_nO
 thermo          1000
 
diff --git a/examples/gcmc/in.gcmc.h2o b/examples/gcmc/in.gcmc.h2o
index 2c03b1ab78..8e408c135b 100644
--- a/examples/gcmc/in.gcmc.h2o
+++ b/examples/gcmc/in.gcmc.h2o
@@ -30,8 +30,6 @@ create_box      2 box                       &
                 extra/angle/per/atom 1      &
                 extra/special/per/atom 2
 
-# we can load multiple molecule templates, but don't have to use them all
-molecule        co2mol CO2.txt
 molecule        h2omol H2O.txt
 create_atoms   	0 box mol h2omol 464563 units box
                         
@@ -58,15 +56,21 @@ timestep        1.0
 
 minimize 0.0 0.0 100 1000
 reset_timestep 0
+
 # rigid constraints with thermostat 
 
-fix             mynvt all nvt temp ${temp} ${temp} 100
-fix             wshake  all shake 0.0001 50 0 b 1 a 1 mol h2omol
-# gcmc
+fix             mynvt h2o nvt temp ${temp} ${temp} 100
+fix             wshake h2o shake 0.0001 50 0 b 1 a 1 mol h2omol
 
+# important to make temperature dofs dynamic
 
+compute_modify  thermo_temp dynamic/dof yes
+compute_modify  mynvt_temp dynamic/dof yes
 
 run 1000
+reset_timestep 0
+
+# gcmc
 
 variable        tfac equal 5.0/3.0 # (3 trans + 2 rot)/(3 trans)
 fix             mygcmc all gcmc 100 100 0 0 54341 ${temp} ${mu} ${disp} mol &
@@ -87,7 +91,6 @@ variable	tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1)
 variable	iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
 variable	dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)
 variable	racc equal f_mygcmc[8]/(f_mygcmc[7]+0.1)
-compute_modify  thermo_temp dynamic/dof yes
 thermo_style    custom step temp press pe ke density atoms v_iacc v_dacc v_tacc v_racc v_nO v_nH
 thermo          1000
 
diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp
index 7ee9cd028e..0c9e776b7d 100644
--- a/src/MC/fix_gcmc.cpp
+++ b/src/MC/fix_gcmc.cpp
@@ -1253,6 +1253,10 @@ void FixGCMC::attempt_molecule_deletion()
 
   if (ngas == 0) return;
 
+  // work-around to avoid n=0 problem with fix rigid/nvt/small
+
+  if (ngas == natoms_per_molecule) return;
+
   tagint deletion_molecule = pick_random_gas_molecule();
   if (deletion_molecule == -1) return;
 
@@ -1910,6 +1914,10 @@ void FixGCMC::attempt_molecule_deletion_full()
 
   if (ngas == 0) return;
 
+  // work-around to avoid n=0 problem with fix rigid/nvt/small
+
+  if (ngas == natoms_per_molecule) return;
+
   tagint deletion_molecule = pick_random_gas_molecule();
   if (deletion_molecule == -1) return;
 
-- 
GitLab