Merge pull request #1764 from jwood13/stable

Add a max and min option to fix_gcmc

doc/src/fix_gcmc.rst

@@ -51,6 +51,8 @@ Syntax
        *intra_energy* value = intramolecular energy (energy units)
        *tfac_insert* value = scale up/down temperature of inserted atoms (unitless)
        *overlap_cutoff* value = maximum pair distance for overlap rejection (distance units)
+       *max* value = Maximum number of molecules allowed in the system
+       *min* value = Minimum number of molecules allowed in the system
 
 
 
@@ -385,6 +387,12 @@ assigning an infinite positive energy to all new configurations that
 place any pair of atoms closer than the specified overlap cutoff
 distance.
 
+The *max* and *min* keywords allow for the restriction of the number
+of atoms in the simulation. They automatically reject all insertion 
+or deletion moves that would take the system beyond the set boundaries.
+Should the system already be beyond the boundary, only moves that bring 
+the system closer to the bounds may be accepted.
+
 The *group* keyword adds all inserted atoms to the
 :doc:`group <group>` of the group-ID value. The *grouptype* keyword
 adds all inserted atoms of the specified type to the

doc/txt/fix_gcmc.txt

@@ -48,7 +48,9 @@ keyword = {mol}, {region}, {maxangle}, {pressure}, {fugacity_coeff}, {full_energ
     group-ID = group-ID for inserted atoms (string)
   {intra_energy} value = intramolecular energy (energy units)
   {tfac_insert} value = scale up/down temperature of inserted atoms (unitless)
-  {overlap_cutoff} value = maximum pair distance for overlap rejection (distance units) :pre
+  {overlap_cutoff} value = maximum pair distance for overlap rejection (distance units)
+  {max} value = Maximum number of molecules allowed in the system
+  {min} value = Minimum number of molecules allowed in the system :pre
 :ule
 
 [Examples:]
@@ -364,6 +366,12 @@ assigning an infinite positive energy to all new configurations that
 place any pair of atoms closer than the specified overlap cutoff
 distance.
 
+The {max} and {min} keywords allow for the restriction of the number
+of atoms in the simulation. They automatically reject all insertion 
+or deletion moves that would take the system beyond the set boundaries.
+Should the system already be beyond the boundary, only moves that bring 
+the system closer to the bounds may be accepted.
+
 The {group} keyword adds all inserted atoms to the
 "group"_group.html of the group-ID value. The {grouptype} keyword
 adds all inserted atoms of the specified type to the

src/MC/fix_gcmc.cpp

@@ -268,6 +268,8 @@ void FixGCMC::options(int narg, char **arg)
   tfac_insert = 1.0;
   overlap_cutoffsq = 0.0;
   overlap_flag = 0;
+  min_ngas = -1;
+  max_ngas = INT_MAX;
 
   int iarg = 0;
   while (iarg < narg) {
@@ -387,6 +389,14 @@ void FixGCMC::options(int narg, char **arg)
       overlap_cutoffsq = rtmp*rtmp;
       overlap_flag = 1;
       iarg += 2;
+    } else if (strcmp(arg[iarg],"min") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix gcmc command");
+      min_ngas = force->numeric(FLERR,arg[iarg+1]);
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"max") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix gcmc command");
+      max_ngas = force->numeric(FLERR,arg[iarg+1]);
+      iarg += 2;
     } else error->all(FLERR,"Illegal fix gcmc command");
   }
 }
@@ -893,7 +903,7 @@ void FixGCMC::attempt_atomic_deletion()
 {
   ndeletion_attempts += 1.0;
 
-  if (ngas == 0) return;
+  if (ngas == 0 || ngas <= min_ngas) return;
 
   int i = pick_random_gas_atom();
 
@@ -934,6 +944,8 @@ void FixGCMC::attempt_atomic_insertion()
 
   ninsertion_attempts += 1.0;
 
+  if (ngas >= max_ngas) return;
+  
   // pick coordinates for insertion point
 
   double coord[3];
@@ -1248,7 +1260,7 @@ void FixGCMC::attempt_molecule_deletion()
 {
   ndeletion_attempts += 1.0;
 
-  if (ngas == 0) return;
+  if (ngas == 0 || ngas <= min_ngas) return;
 
   // work-around to avoid n=0 problem with fix rigid/nvt/small
 
@@ -1287,6 +1299,8 @@ void FixGCMC::attempt_molecule_insertion()
   double lamda[3];
   ninsertion_attempts += 1.0;
 
+  if (ngas >= max_ngas) return;
+  
   double com_coord[3];
   if (regionflag) {
     int region_attempt = 0;
@@ -1570,7 +1584,7 @@ void FixGCMC::attempt_atomic_deletion_full()
 
   ndeletion_attempts += 1.0;
 
-  if (ngas == 0) return;
+  if (ngas == 0 || ngas <= min_ngas) return;
 
   double energy_before = energy_stored;
 
@@ -1619,6 +1633,8 @@ void FixGCMC::attempt_atomic_insertion_full()
   double lamda[3];
   ninsertion_attempts += 1.0;
 
+  if (ngas >= max_ngas) return;
+  
   double energy_before = energy_stored;
 
   double coord[3];
@@ -1914,7 +1930,7 @@ void FixGCMC::attempt_molecule_deletion_full()
 {
   ndeletion_attempts += 1.0;
 
-  if (ngas == 0) return;
+  if (ngas == 0 || ngas <= min_ngas) return;
 
   // work-around to avoid n=0 problem with fix rigid/nvt/small
 
@@ -1997,6 +2013,8 @@ void FixGCMC::attempt_molecule_insertion_full()
   double lamda[3];
   ninsertion_attempts += 1.0;
 
+  if (ngas >= max_ngas) return;
+
   double energy_before = energy_stored;
 
   tagint maxmol = 0;

src/MC/fix_gcmc.h

@@ -119,6 +119,8 @@ class FixGCMC : public Fix {
   imageint imagezero;
   double overlap_cutoffsq; // square distance cutoff for overlap
   int overlap_flag;
+  int max_ngas;
+  int min_ngas;
 
   double energy_intra;