Merge pull request #800 from athomps/fix_gcmc_segfault_fix

Added warning to discourage use of group all and fixed some segfault …

doc/src/fix_gcmc.txt

@@ -95,8 +95,8 @@ which will result in roughly one MC move per atom or molecule
 per MC cycle.
 
 All inserted particles are always added to two groups: the default
-group "all" and the fix group specified in the fix command (which can
-also be "all"). In addition, particles are also added to any groups
+group "all" and the fix group specified in the fix command.
+In addition, particles are also added to any groups
 specified by the {group} and {grouptype} keywords.  If inserted
 particles are individual atoms, they are assigned the atom type given
 by the type argument.  If they are molecules, the type argument has no
@@ -104,6 +104,12 @@ effect and must be set to zero. Instead, the type of each atom in the
 inserted molecule is specified in the file read by the
 "molecule"_molecule.html command.
 
+NOTE: Care should be taken to apply fix gcmc only to
+a group that contains only those atoms and molecules
+that you wish to manipulate using Monte Carlo.
+Hence it is generally not a good idea to specify
+the default group "all" in the fix command, although it is allowed.
+
 This fix cannot be used to perform GCMC insertions of gas atoms or
 molecules other than the exchanged type, but GCMC deletions,
 and MC translations, and rotations can be performed on any atom/molecule in

examples/gcmc/in.gcmc.co2

@@ -58,21 +58,21 @@ timestep        1.0
 
 # rigid constraints with thermostat 
 
-fix             myrigidnvt all rigid/nvt/small molecule temp ${temp} ${temp} 100 mol co2mol
+fix             myrigidnvt co2 rigid/nvt/small molecule temp ${temp} ${temp} 100 mol co2mol
 fix_modify	myrigidnvt dynamic/dof no
 
 # 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 &
+fix             mygcmc co2 gcmc 100 100 0 0 54341 ${temp} ${mu} ${disp} mol &
                 co2mol tfac_insert ${tfac} group co2 rigid myrigidnvt
 
 # atom counts
 
 variable 	carbon atom "type==1"
 variable        oxygen atom "type==2"
-group 		carbon dynamic all var carbon
-group 	        oxygen dynamic all var oxygen
+group 		carbon dynamic co2 var carbon
+group 	        oxygen dynamic co2 var oxygen
 variable        nC equal count(carbon)
 variable        nO equal count(oxygen)
 

examples/gcmc/in.gcmc.lj

@@ -29,14 +29,18 @@ create_box	1 box
 pair_coeff	* * 1.0 1.0
 mass		* 1.0
 
+# we recommend setting up a dedicated group for gcmc
+
+group		gcmcgroup type 1
+
 # gcmc
 
-fix             mygcmc all gcmc 1 100 100 1 29494 ${temp} ${mu} ${disp}
+fix             mygcmc gcmcgroup gcmc 1 100 100 1 29494 ${temp} ${mu} ${disp}
 
 # atom count
 
 variable 	type1 atom "type==1"
-group 		type1 dynamic all var type1
+group 		type1 dynamic gcmcgroup var type1
 variable        n1 equal count(type1)
 
 # averaging

src/MC/fix_gcmc.cpp

@@ -72,7 +72,8 @@ enum{MOVEATOM,MOVEMOL}; // movemode
 FixGCMC::FixGCMC(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg),
   idregion(NULL), full_flag(0), ngroups(0), groupstrings(NULL), ngrouptypes(0), grouptypestrings(NULL),
-  grouptypebits(NULL), grouptypes(NULL), local_gas_list(NULL), molcoords(NULL), random_equal(NULL), random_unequal(NULL),
+  grouptypebits(NULL), grouptypes(NULL), local_gas_list(NULL), molcoords(NULL), molq(NULL), molimage(NULL),
+  random_equal(NULL), random_unequal(NULL),
   fixrigid(NULL), fixshake(NULL), idrigid(NULL), idshake(NULL)
 {
   if (narg < 11) error->all(FLERR,"Illegal fix gcmc command");
@@ -199,8 +200,8 @@ FixGCMC::FixGCMC(LAMMPS *lmp, int narg, char **arg) :
   
   if (exchmode == EXCHATOM) natoms_per_molecule = 1;
   else natoms_per_molecule = onemols[imol]->natoms;
-  nmaxmolcoords = natoms_per_molecule;
-  memory->create(molcoords,nmaxmolcoords,3,"gcmc:molcoords");
+  nmaxmolatoms = natoms_per_molecule;
+  grow_molecule_arrays(nmaxmolatoms);
 
   // compute the number of MC cycles that occur nevery timesteps
 
@@ -513,7 +514,7 @@ void FixGCMC::init()
        "Fix gcmc cannot exchange individual atoms belonging to a molecule");
   }
 
-  // if molecules are exchanged are moves, check for unset mol IDs
+  // if molecules are exchanged or moved, check for unset mol IDs
 
   if (exchmode == EXCHMOL || movemode == MOVEMOL) {
     tagint *molecule = atom->molecule;
@@ -683,6 +684,12 @@ void FixGCMC::init()
   imagezero = ((imageint) IMGMAX << IMG2BITS) |
              ((imageint) IMGMAX << IMGBITS) | IMGMAX;
 
+  // warning if group id is "all"
+
+  if (groupbit & 1)
+    error->warning(FLERR, "Fix gcmc is being applied "
+                   "to the default group all");
+    
   // construct group bitmask for all new atoms
   // aggregated over all group keywords
 
@@ -1153,10 +1160,8 @@ void FixGCMC::attempt_molecule_rotation()
     }
   }
 
-  if (nmolcoords > nmaxmolcoords) {
-    nmaxmolcoords = nmolcoords;
-    molcoords = memory->grow(molcoords,nmaxmolcoords,3,"gcmc:molcoords");
-  }
+  if (nmolcoords > nmaxmolatoms)
+    grow_molecule_arrays(nmolcoords);
   
   double com[3];
   com[0] = com[1] = com[2] = 0.0;
@@ -1816,10 +1821,8 @@ void FixGCMC::attempt_molecule_rotation_full()
     }
   }
 
-  if (nmolcoords > nmaxmolcoords) {
-    nmaxmolcoords = nmolcoords;
-    molcoords = memory->grow(molcoords,nmaxmolcoords,3,"gcmc:molcoords");
-  }
+  if (nmolcoords > nmaxmolatoms)
+    grow_molecule_arrays(nmolcoords);
   
   double com[3];
   com[0] = com[1] = com[2] = 0.0;
@@ -1844,14 +1847,14 @@ void FixGCMC::attempt_molecule_rotation_full()
 
   double **x = atom->x;
   imageint *image = atom->image;
-  imageint image_orig[natoms_per_molecule];
+  
   int n = 0;
   for (int i = 0; i < atom->nlocal; i++) {
     if (mask[i] & molecule_group_bit) {
       molcoords[n][0] = x[i][0];
       molcoords[n][1] = x[i][1];
       molcoords[n][2] = x[i][2];
-      image_orig[n] = image[i];
+      molimage[n] = image[i];
       double xtmp[3];
       domain->unmap(x[i],image[i],xtmp);
       xtmp[0] -= com[0];
@@ -1884,7 +1887,7 @@ void FixGCMC::attempt_molecule_rotation_full()
         x[i][0] = molcoords[n][0];
         x[i][1] = molcoords[n][1];
         x[i][2] = molcoords[n][2];
-        image[i] = image_orig[n];
+        image[i] = molimage[n];
         n++;
       }
     }
@@ -1906,8 +1909,17 @@ void FixGCMC::attempt_molecule_deletion_full()
 
   double energy_before = energy_stored;
 
+  // check nmolq, grow arrays if necessary
+
+  int nmolq = 0;
+  for (int i = 0; i < atom->nlocal; i++)
+    if (atom->molecule[i] == deletion_molecule)
+      if (atom->q_flag) nmolq++;
+  
+  if (nmolq > nmaxmolatoms)
+    grow_molecule_arrays(nmolq);
+  
   int m = 0;
-  double q_tmp[natoms_per_molecule];
   int tmpmask[atom->nlocal];
   for (int i = 0; i < atom->nlocal; i++) {
     if (atom->molecule[i] == deletion_molecule) {
@@ -1915,7 +1927,7 @@ void FixGCMC::attempt_molecule_deletion_full()
       atom->mask[i] = exclusion_group_bit;
       toggle_intramolecular(i);
       if (atom->q_flag) {
-        q_tmp[m] = atom->q[i];
+        molq[m] = atom->q[i];
         m++;
         atom->q[i] = 0.0;
       }
@@ -1948,7 +1960,7 @@ void FixGCMC::attempt_molecule_deletion_full()
         atom->mask[i] = tmpmask[i];
         toggle_intramolecular(i);
         if (atom->q_flag) {
-          atom->q[i] = q_tmp[m];
+          atom->q[i] = molq[m];
           m++;
         }
       }
@@ -2521,3 +2533,10 @@ void FixGCMC::restart(char *buf)
 
   next_reneighbor = static_cast<int> (list[n++]);
 }
+
+void FixGCMC::grow_molecule_arrays(int nmolatoms) {
+    nmaxmolatoms = nmolatoms;
+    molcoords = memory->grow(molcoords,nmaxmolatoms,3,"gcmc:molcoords");
+    molq = memory->grow(molq,nmaxmolatoms,"gcmc:molq");
+    molimage = memory->grow(molimage,nmaxmolatoms,"gcmc:molimage");
+}

src/MC/fix_gcmc.h

@@ -57,7 +57,8 @@ class FixGCMC : public Fix {
   double memory_usage();
   void write_restart(FILE *);
   void restart(char *);
-
+  void grow_molecule_arrays(int);
+  
  private:
   int molecule_group,molecule_group_bit;
   int molecule_group_inversebit;
@@ -78,7 +79,7 @@ class FixGCMC : public Fix {
   bool full_flag;           // true if doing full system energy calculations
 
   int natoms_per_molecule;  // number of atoms in each inserted molecule
-  int nmaxmolcoords;        // number of atoms allocated for molcoords
+  int nmaxmolatoms;         // number of atoms allocated for molecule arrays
 
   int groupbitall;          // group bitmask for inserted atoms
   int ngroups;              // number of group-ids for inserted atoms
@@ -114,6 +115,8 @@ class FixGCMC : public Fix {
   int *local_gas_list;
   double **cutsq;
   double **molcoords;
+  double *molq;
+  imageint *molimage;
   imageint imagezero;
   double overlap_cutoffsq; // square distance cutoff for overlap 
   int overlap_flag;
@@ -224,6 +227,12 @@ W: Energy of old configuration in fix gcmc is > MAXENERGYTEST.
 This probably means that a pair of atoms are closer than the 
 overlap cutoff distance for keyword overlap_cutoff.
 
+W: Fix gcmc is being applied to the default group all
+
+This is allowed, but it will result in Monte Carlo moves
+being performed on all the atoms in the system, which is
+often not what is intended.
+
 E: Invalid atom type in fix gcmc command
 
 The atom type specified in the gcmc command does not exist.