Fixed recent segfault in fix gcmc and added mcmoves keyword

doc/src/fix_gcmc.txt

@@ -26,16 +26,20 @@ zero or more keyword/value pairs may be appended to args :l
 keyword = {mol}, {region}, {maxangle}, {pressure}, {fugacity_coeff}, {full_energy}, {charge}, {group}, {grouptype}, {intra_energy}, {tfac_insert}, or {overlap_cutoff}
   {mol} value = template-ID
     template-ID = ID of molecule template specified in a separate "molecule"_molecule.html command
+  {mcmoves} values = Patomtrans Pmoltrans Pmolrotate
+    Patomtrans = proportion of atom translation MC moves
+    Pmoltrans = proportion of molecule translation MC moves
+    Pmolrotate = proportion of molecule rotation MC moves
   {rigid} value = fix-ID
     fix-ID = ID of "fix rigid/small"_fix_rigid.html command
   {shake} value = fix-ID
     fix-ID = ID of "fix shake"_fix_shake.html command
   {region} value = region-ID
-    region-ID = ID of region where MC moves are allowed
+    region-ID = ID of region where GCMC exchanges and MC moves are allowed
   {maxangle} value = maximum molecular rotation angle (degrees)
   {pressure} value = pressure of the gas reservoir (pressure units)
   {fugacity_coeff} value = fugacity coefficient of the gas reservoir (unitless)
-  {full_energy} = compute the entire system energy when performing MC moves
+  {full_energy} = compute the entire system energy when performing GCMC exchanges and MC moves
   {charge} value = charge of inserted atoms (charge units)
   {group} value = group-ID
     group-ID = group-ID for inserted atoms (string)
@@ -56,34 +60,42 @@ fix 4 my_gas gcmc 1 10 10 1 123456543 300.0 -12.5 1.0 region disk :pre
 [Description:]
 
 This fix performs grand canonical Monte Carlo (GCMC) exchanges of
-atoms or molecules of the given type with an imaginary ideal gas
+atoms or molecules with an imaginary ideal gas
 reservoir at the specified T and chemical potential (mu) as discussed
-in "(Frenkel)"_#Frenkel. If used with the "fix nvt"_fix_nh.html
+in "(Frenkel)"_#Frenkel. It also
+attempts  Monte Carlo (MC) moves (translations and molecule
+rotations) within the simulation cell or
+region. If used with the "fix nvt"_fix_nh.html
 command, simulations in the grand canonical ensemble (muVT, constant
 chemical potential, constant volume, and constant temperature) can be
 performed.  Specific uses include computing isotherms in microporous
 materials, or computing vapor-liquid coexistence curves.
 
-Every N timesteps the fix attempts a number of GCMC exchanges
+Every N timesteps the fix attempts both GCMC exchanges
-(insertions or deletions) of gas atoms or molecules of the given type
+(insertions or deletions) and MC moves of gas atoms or molecules.
-between the simulation cell and the imaginary reservoir. It also
+On those timesteps, the average number of attempted GCMC exchanges is X,
-attempts a number of Monte Carlo moves (translations and molecule
+while the average number of attempted MC moves is M. 
-rotations) of gas of the given type within the simulation cell or
+For GCMC exchanges of either molecular or atomic gasses,
-region.  The average number of attempted GCMC exchanges is X. The
+these exchanges can be either deletions or insertions,
-average number of attempted MC moves is M.  M should typically be
+with equal probability.
+
+The possible choices for MC moves are translation of an atom,
+translation of a molecule, and rotation of a molecule.
+The relative amounts of each are determined by the optional
+{mcmoves} keyword (see below).
+The default behavior is as follows.
+If the {mol} keyword is used, only molecule translations
+and molecule rotations are performed with equal probability.
+Conversely, if the {mol} keyword is not used, only atom
+translations are performed.
+M should typically be
 chosen to be approximately equal to the expected number of gas atoms
 or molecules of the given type within the simulation cell or region,
-which will result in roughly one MC translation per atom or molecule
+which will result in roughly one MC move per atom or molecule
 per MC cycle.
 
-For MC moves of molecular gasses, rotations and translations are each
+All inserted particles are always added to two groups: the default
-attempted with 50% probability. For MC moves of atomic gasses,
+group "all" and the fix group specified in the fix command (which can
-translations are attempted 100% of the time. For MC exchanges of
-either molecular or atomic gasses, deletions and insertions are each
-attempted with 50% probability.
-
-All inserted particles are always assigned to two groups: the default
-group "all" and the group specified in the fix gcmc command (which can
 also be "all"). 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
@@ -92,9 +104,9 @@ 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.
 
-This fix cannot be used to perform MC insertions of gas atoms or
+This fix cannot be used to perform GCMC insertions of gas atoms or
-molecules other than the exchanged type, but MC deletions,
+molecules other than the exchanged type, but GCMC deletions,
-translations, and rotations can be performed on any atom/molecule in
+and MC translations, and rotations can be performed on any atom/molecule in
 the fix group.  All atoms in the simulation cell can be moved using
 regular time integration translations, e.g. via "fix nvt"_fix_nh.html,
 resulting in a hybrid GCMC+MD simulation. A smaller-than-usual
@@ -150,8 +162,8 @@ about this point.
 Individual atoms are inserted, unless the {mol} keyword is used.  It
 specifies a {template-ID} previously defined using the
 "molecule"_molecule.html command, which reads a file that defines the
-molecule.  The coordinates, atom types, charges, etc, as well as any
+molecule.  The coordinates, atom types, charges, etc., as well as any
-bond/angle/etc and special neighbor information for the molecule can
+bonding and special neighbor information for the molecule can
 be specified in the molecule file.  See the "molecule"_molecule.html
 command for details.  The only settings required to be in this file
 are the coordinates and types of atoms in the molecule.
@@ -162,7 +174,7 @@ error when it tries to find bonded neighbors.  LAMMPS will warn you if
 any of the atoms eligible for deletion have a non-zero molecule ID,
 but does not check for this at the time of deletion.
 
-If you wish to insert molecules via the {mol} keyword, that will be
+If you wish to insert molecules using the {mol} keyword that will be
 treated as rigid bodies, use the {rigid} keyword, specifying as its
 value the ID of a separate "fix rigid/small"_fix_rigid.html command
 which also appears in your input script.
@@ -178,6 +190,15 @@ their bonds or angles constrained via SHAKE, use the {shake} keyword,
 specifying as its value the ID of a separate "fix
 shake"_fix_shake.html command which also appears in your input script.
 
+Optionally, users may specify the relative amounts of different MC
+moves using the {mcmoves} keyword. The values {Patomtrans},
+{Pmoltrans}, {Pmolrotate} specify the average proportion of
+atom translations, molecule translations, and molecule rotations,
+respectively. The values must be non-negative integers or real
+numbers, with at least one non-zero value. For example, (10,30,0)
+would result in 25% of the MC moves being atomic translations, 75%
+molecular translations, and no molecular rotations. 
+
 Optionally, users may specify the maximum rotation angle for molecular
 rotations using the {maxangle} keyword and specifying the angle in
 degrees. Rotations are performed by generating a random point on the
@@ -188,19 +209,19 @@ to the unit vector defined by the point on the unit sphere.  The same
 procedure is used for randomly rotating molecules when they are
 inserted, except that the maximum angle is 360 degrees.
 
-Note that fix GCMC does not use configurational bias MC or any other
+Note that fix gcmc does not use configurational bias MC or any other
 kind of sampling of intramolecular degrees of freedom.  Inserted
 molecules can have different orientations, but they will all have the
 same intramolecular configuration, which was specified in the molecule
 command input.
 
 For atomic gasses, inserted atoms have the specified atom type, but
-deleted atoms are any atoms that have been inserted or that belong to
+deleted atoms are any atoms that have been inserted or that already
-the user-specified fix group. For molecular gasses, exchanged
+belong to the fix group. For molecular gasses, exchanged
 molecules use the same atom types as in the template molecule supplied
 by the user.  In both cases, exchanged atoms/molecules are assigned to
-two groups: the default group "all" and the group specified in the fix
+two groups: the default group "all" and the fix group
-gcmc command (which can also be "all").
+(which can also be "all").
 
 The chemical potential is a user-specified input parameter defined
 as:
@@ -241,13 +262,16 @@ which case the user-specified chemical potential is ignored. The user
 may also specify the fugacity coefficient phi using the
 {fugacity_coeff} keyword, which defaults to unity.
 
-The {full_energy} option means that fix GCMC will compute the total
+The {full_energy} option means that the fix calculates the total
-potential energy of the entire simulated system. The total system
+potential energy of the entire simulated system, instead of just
-energy before and after the proposed GCMC move is then used in the
+the energy of the part that is changed. The total system
+energy before and after the proposed GCMC exchange or MC move
+is then used in the
 Metropolis criterion to determine whether or not to accept the
-proposed GCMC move. By default, this option is off, in which case only
+proposed change. By default, this option is off,
-partial energies are computed to determine the difference in energy
+in which case only
-that would be caused by the proposed GCMC move.
+partial energies are computed to determine the energy difference
+due to the proposed change.
 
 The {full_energy} option is needed for systems with complicated
 potential energy calculations, including the following:
@@ -263,10 +287,11 @@ In these cases, LAMMPS will automatically apply the {full_energy}
 keyword and issue a warning message.
 
 When the {mol} keyword is used, the {full_energy} option also includes
-the intramolecular energy of inserted and deleted molecules. If this
+the intramolecular energy of inserted and deleted molecules, whereas
+this energy is not included when {full_energy} is not used. If this
 is not desired, the {intra_energy} keyword can be used to define an
 amount of energy that is subtracted from the final energy when a
-molecule is inserted, and added to the initial energy when a molecule
+molecule is inserted, and subtracted from the initial energy when a molecule
 is deleted. For molecules that have a non-zero intramolecular energy,
 this will ensure roughly the same behavior whether or not the
 {full_energy} option is used.
@@ -291,7 +316,8 @@ include: "efield"_fix_efield.html, "gravity"_fix_gravity.html,
 "temp/berendsen"_fix_temp_berendsen.html,
 "temp/rescale"_fix_temp_rescale.html, and "wall fixes"_fix_wall.html.
 For that energy to be included in the total potential energy of the
-system (the quantity used when performing GCMC moves), you MUST enable
+system (the quantity used when performing GCMC exchange and MC moves),
+you MUST enable
 the "fix_modify"_fix_modify.html {energy} option for that fix.  The
 doc pages for individual "fix"_fix.html commands specify if this
 should be done.
@@ -305,9 +331,14 @@ about simulating non-neutral systems with kspace on.
 
 Use of this fix typically will cause the number of atoms to fluctuate,
 therefore, you will want to use the
-"compute_modify"_compute_modify.html command to insure that the
+"compute_modify dynamic/dof"_compute_modify.html command to insure that the
 current number of atoms is used as a normalizing factor each time
-temperature is computed.  Here is the necessary command:
+temperature is computed. A simple example of this is:
+
+compute_modify thermo_temp dynamic yes :pre
+
+A more complicated example is listed earlier on this page
+in the context of NVT dynamics.
 
 NOTE: If the density of the cell is initially very small or zero, and
 increases to a much larger density after a period of equilibration,
@@ -327,17 +358,9 @@ assigning an infinite positive energy to all new configurations that
 place any pair of atoms closer than the specified overlap cutoff
 distance.
 
-compute_modify thermo_temp dynamic yes :pre
+The {group} keyword adds all inserted atoms to the
-
-If LJ units are used, note that a value of 0.18292026 is used by this
-fix as the reduced value for Planck's constant.  This value was
-derived from LJ parameters for argon, where h* = h/sqrt(sigma^2 *
-epsilon * mass), sigma = 3.429 angstroms, epsilon/k = 121.85 K, and
-mass = 39.948 amu.
-
-The {group} keyword assigns all inserted atoms to the
 "group"_group.html of the group-ID value. The {grouptype} keyword
-assigns all inserted atoms of the specified type to the
+adds all inserted atoms of the specified type to the
 "group"_group.html of the group-ID value.
 
 [Restart, fix_modify, output, run start/stop, minimize info:]
@@ -384,7 +407,8 @@ Can be run in parallel, but aspects of the GCMC part will not scale
 well in parallel. Only usable for 3D simulations.
 
 When using fix gcmc in combination with fix shake or fix rigid,
-only gcmc exchange moves are supported.
+only GCMC exchange moves are supported, so the argument
+{M} must be zero.
 
 Note that very lengthy simulations involving insertions/deletions of
 billions of gas molecules may run out of atom or molecule IDs and
@@ -409,7 +433,9 @@ the user for each subsequent fix gcmc command.
 [Default:]
 
 The option defaults are mol = no, maxangle = 10, overlap_cutoff = 0.0,
-fugacity_coeff = 1, and full_energy = no,
+fugacity_coeff = 1.0, intra_energy = 0.0, tfac_insert = 1.0.
+(Patomtrans, Pmoltrans, Pmolrotate) = (1, 0, 0) for mol = no and
+(0, 1, 1) for mol = yes. full_energy = no,
 except for the situations where full_energy is required, as
 listed above.
 

src/MC/fix_gcmc.cpp

@@ -64,15 +64,16 @@ using namespace MathConst;
 
 #define MAXENERGYTEST 1.0e50
 
-enum{ATOM,MOLECULE};
+enum{EXCHATOM,EXCHMOL}; // exchmode
+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), atom_coord(NULL), random_equal(NULL), random_unequal(NULL),
+  grouptypebits(NULL), grouptypes(NULL), local_gas_list(NULL), molcoords(NULL), random_equal(NULL), random_unequal(NULL),
-  coords(NULL), imageflags(NULL), fixrigid(NULL), fixshake(NULL), idrigid(NULL), idshake(NULL)
+  fixrigid(NULL), fixshake(NULL), idrigid(NULL), idshake(NULL)
 {
   if (narg < 11) error->all(FLERR,"Illegal fix gcmc command");
 
@@ -161,9 +162,9 @@ FixGCMC::FixGCMC(LAMMPS *lmp, int narg, char **arg) :
      static_cast<double> (attempts);
   }
 
-  // error check and further setup for mode = MOLECULE
+  // error check and further setup for exchmode = EXCHMOL
 
-  if (mode == MOLECULE) {
+  if (exchmode == EXCHMOL) {
     if (onemols[imol]->xflag == 0)
       error->all(FLERR,"Fix gcmc molecule must have coordinates");
     if (onemols[imol]->typeflag == 0)
@@ -182,9 +183,9 @@ FixGCMC::FixGCMC(LAMMPS *lmp, int narg, char **arg) :
   if (charge_flag && atom->q == NULL)
     error->all(FLERR,"Fix gcmc atom has charge, but atom style does not");
 
-  if (rigidflag && mode == ATOM)
+  if (rigidflag && exchmode == EXCHATOM)
     error->all(FLERR,"Cannot use fix gcmc rigid and not molecule");
-  if (shakeflag && mode == ATOM)
+  if (shakeflag && exchmode == EXCHATOM)
     error->all(FLERR,"Cannot use fix gcmc shake and not molecule");
   if (rigidflag && shakeflag)
     error->all(FLERR,"Cannot use fix gcmc rigid and shake");
@@ -193,13 +194,13 @@ FixGCMC::FixGCMC(LAMMPS *lmp, int narg, char **arg) :
   if (shakeflag && (nmcmoves > 0))
     error->all(FLERR,"Cannot use fix gcmc shake with MC moves");
 
-  // setup of coords and imageflags array
+  // setup of array of coordinates for molecule insertion
-
+  // also used by rotation moves for any molecule
-  if (mode == ATOM) natoms_per_molecule = 1;
+  
+  if (exchmode == EXCHATOM) natoms_per_molecule = 1;
   else natoms_per_molecule = onemols[imol]->natoms;
-  memory->create(coords,natoms_per_molecule,3,"gcmc:coords");
+  nmaxmolcoords = natoms_per_molecule;
-  memory->create(imageflags,natoms_per_molecule,"gcmc:imageflags");
+  memory->create(molcoords,nmaxmolcoords,3,"gcmc:molcoords");
-  memory->create(atom_coord,natoms_per_molecule,3,"gcmc:atom_coord");
 
   // compute the number of MC cycles that occur nevery timesteps
 
@@ -235,7 +236,12 @@ void FixGCMC::options(int narg, char **arg)
 
   // defaults
 
-  mode = ATOM;
+  exchmode = EXCHATOM;
+  movemode = MOVEATOM;
+  patomtrans = 0.0;
+  pmoltrans = 0.0;
+  pmolrotate = 0.0;
+  pmctot = 0.0;
   max_rotation_angle = 10*MY_PI/180;
   regionflag = 0;
   iregion = -1;
@@ -277,10 +283,21 @@ void FixGCMC::options(int narg, char **arg)
       if (atom->molecules[imol]->nset > 1 && comm->me == 0)
         error->warning(FLERR,"Molecule template for "
                        "fix gcmc has multiple molecules");
-      mode = MOLECULE;
+      exchmode = EXCHMOL;
       onemols = atom->molecules;
       nmol = onemols[imol]->nset;
       iarg += 2;
+  } else if (strcmp(arg[iarg],"mcmoves") == 0) {
+      if (iarg+4 > narg) error->all(FLERR,"Illegal fix gcmc command");
+      patomtrans = force->numeric(FLERR,arg[iarg+1]);
+      pmoltrans = force->numeric(FLERR,arg[iarg+2]);
+      pmolrotate = force->numeric(FLERR,arg[iarg+3]);
+      if (patomtrans < 0 || pmoltrans < 0 || pmolrotate < 0)
+        error->all(FLERR,"Illegal fix gcmc command");
+      pmctot = patomtrans + pmoltrans + pmolrotate;
+      if (pmctot <= 0)
+        error->all(FLERR,"Illegal fix gcmc command");
+      iarg += 4;
     } else if (strcmp(arg[iarg],"region") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix gcmc command");
       iregion = domain->find_region(arg[iarg+1]);
@@ -387,9 +404,7 @@ FixGCMC::~FixGCMC()
   delete random_unequal;
 
   memory->destroy(local_gas_list);
-  memory->destroy(atom_coord);
+  memory->destroy(molcoords);
-  memory->destroy(coords);
-  memory->destroy(imageflags);
 
   delete [] idrigid;
   delete [] idshake;
@@ -430,6 +445,30 @@ void FixGCMC::init()
 
   triclinic = domain->triclinic;
 
+  // set probabilities for MC moves
+
+  if (pmctot == 0.0)
+    if (exchmode == EXCHATOM) {
+      movemode = MOVEATOM;
+      patomtrans = 1.0;
+      pmoltrans = 0.0;
+      pmolrotate = 0.0;
+    } else {
+      movemode = MOVEMOL;
+      patomtrans = 0.0;
+      pmoltrans = 0.5;
+      pmolrotate = 0.5;
+   }
+  else {
+    if (pmoltrans == 0.0 && pmolrotate == 0.0)
+      movemode = MOVEATOM;
+    else
+      movemode = MOVEMOL;
+    patomtrans /= pmctot;
+    pmoltrans /= pmctot;
+    pmolrotate /= pmctot;
+  }
+  
   // decide whether to switch to the full_energy option
 
   if (!full_flag) {
@@ -454,14 +493,14 @@ void FixGCMC::init()
 
   int *type = atom->type;
 
-  if (mode == ATOM) {
+  if (exchmode == EXCHATOM) {
     if (ngcmc_type <= 0 || ngcmc_type > atom->ntypes)
       error->all(FLERR,"Invalid atom type in fix gcmc command");
   }
 
-  // if mode == ATOM, warn if any deletable atom has a mol ID
+  // if atoms are exchanged, warn if any deletable atom has a mol ID
 
-  if ((mode == ATOM) && atom->molecule_flag) {
+  if ((exchmode == EXCHATOM) && atom->molecule_flag) {
     tagint *molecule = atom->molecule;
     int flag = 0;
     for (int i = 0; i < atom->nlocal; i++)
@@ -474,9 +513,9 @@ void FixGCMC::init()
        "Fix gcmc cannot exchange individual atoms belonging to a molecule");
   }
 
-  // if mode == MOLECULE, check for unset mol IDs
+  // if molecules are exchanged are moves, check for unset mol IDs
 
-  if (mode == MOLECULE) {
+  if (exchmode == EXCHMOL || movemode == MOVEMOL) {
     tagint *molecule = atom->molecule;
     int *mask = atom->mask;
     int flag = 0;
@@ -490,11 +529,11 @@ void FixGCMC::init()
        "All mol IDs should be set for fix gcmc group atoms");
   }
 
-  if (((mode == MOLECULE) && (atom->molecule_flag == 0)) ||
+  if (exchmode == EXCHMOL || movemode == MOVEMOL)
-      ((mode == MOLECULE) && (!atom->tag_enable || !atom->map_style)))
+    if (atom->molecule_flag == 0 || !atom->tag_enable || !atom->map_style)
-    error->all(FLERR,
+      error->all(FLERR,
-               "Fix gcmc molecule command requires that "
+       "Fix gcmc molecule command requires that "
-               "atoms have molecule attributes");
+       "atoms have molecule attributes");
 
   // if rigidflag defined, check for rigid/small fix
   // its molecule template must be same as this one
@@ -566,7 +605,7 @@ void FixGCMC::init()
 
   // create a new group for temporary use with selected molecules
 
-  if (mode == MOLECULE) {
+  if (exchmode == EXCHMOL || movemode == MOVEMOL) {
     char **group_arg = new char*[3];
     // create unique group name for atoms to be rotated
     int len = strlen(id) + 30;
@@ -586,10 +625,10 @@ void FixGCMC::init()
     delete [] group_arg;
   }
 
-  // get all of the needed molecule data if mode == MOLECULE,
+  // get all of the needed molecule data if exchanging
-  // otherwise just get the gas mass
+  // or moving molecules, otherwise just get the gas mass
 
-  if (mode == MOLECULE) {
+  if (exchmode == EXCHMOL || movemode == MOVEMOL) {
 
     onemols[imol]->compute_mass();
     onemols[imol]->compute_com();
@@ -713,27 +752,21 @@ void FixGCMC::pre_exchange()
         error->warning(FLERR,"Energy of old configuration in "
                        "fix gcmc is > MAXENERGYTEST.");
 
-    if (mode == MOLECULE) {
+    for (int i = 0; i < ncycles; i++) {
-      for (int i = 0; i < ncycles; i++) {
+      int ixm = static_cast<int>(random_equal->uniform()*ncycles) + 1;
-        int random_int_fraction =
+      if (ixm <= nmcmoves) {
-          static_cast<int>(random_equal->uniform()*ncycles) + 1;
+        double xmcmove = random_equal->uniform();
-        if (random_int_fraction <= nmcmoves) {
+        if (xmcmove < patomtrans) attempt_atomic_translation_full();
-          if (random_equal->uniform() < 0.5) attempt_molecule_translation_full();
+        else if (xmcmove < patomtrans+pmoltrans) attempt_molecule_translation_full();
-          else attempt_molecule_rotation_full();
+        else attempt_molecule_rotation_full();
-        } else {
+      } else {
-          if (random_equal->uniform() < 0.5) attempt_molecule_deletion_full();
+        double xgcmc = random_equal->uniform();
-          else attempt_molecule_insertion_full();
+        if (exchmode == EXCHATOM) {
-        }
+          if (xgcmc < 0.5) attempt_atomic_deletion_full();
-      }
-    } else {
-      for (int i = 0; i < ncycles; i++) {
-        int random_int_fraction =
-          static_cast<int>(random_equal->uniform()*ncycles) + 1;
-        if (random_int_fraction <= nmcmoves) {
-          attempt_atomic_translation_full();
-        } else {
-          if (random_equal->uniform() < 0.5) attempt_atomic_deletion_full();
           else attempt_atomic_insertion_full();
+        } else {
+          if (xgcmc < 0.5) attempt_molecule_deletion_full();
+          else attempt_molecule_insertion_full();
         }
       }
     }
@@ -746,27 +779,21 @@ void FixGCMC::pre_exchange()
 
   } else {
 
-    if (mode == MOLECULE) {
+    for (int i = 0; i < ncycles; i++) {
-      for (int i = 0; i < ncycles; i++) {
+      int ixm = static_cast<int>(random_equal->uniform()*ncycles) + 1;
-        int random_int_fraction =
+      if (ixm <= nmcmoves) {
-          static_cast<int>(random_equal->uniform()*ncycles) + 1;
+        double xmcmove = random_equal->uniform();
-        if (random_int_fraction <= nmcmoves) {
+        if (xmcmove < patomtrans) attempt_atomic_translation();
-          if (random_equal->uniform() < 0.5) attempt_molecule_translation();
+        else if (xmcmove < patomtrans+pmoltrans) attempt_molecule_translation();
-          else attempt_molecule_rotation();
+        else attempt_molecule_rotation();
-        } else {
+      } else {
-          if (random_equal->uniform() < 0.5) attempt_molecule_deletion();
+        double xgcmc = random_equal->uniform();
-          else attempt_molecule_insertion();
+        if (exchmode == EXCHATOM) {
-        }
+          if (xgcmc < 0.5) attempt_atomic_deletion();
-      }
-    } else {
-      for (int i = 0; i < ncycles; i++) {
-        int random_int_fraction =
-          static_cast<int>(random_equal->uniform()*ncycles) + 1;
-        if (random_int_fraction <= nmcmoves) {
-          attempt_atomic_translation();
-        } else {
-          if (random_equal->uniform() < 0.5) attempt_atomic_deletion();
           else attempt_atomic_insertion();
+        } else {
+          if (xgcmc < 0.5) attempt_molecule_deletion();
+          else attempt_molecule_insertion();
         }
       }
     }
@@ -1116,14 +1143,21 @@ void FixGCMC::attempt_molecule_rotation()
                    "fix gcmc is > MAXENERGYTEST.");
 
   int *mask = atom->mask;
+  int nmolcoords = 0;
   for (int i = 0; i < atom->nlocal; i++) {
     if (atom->molecule[i] == rotation_molecule) {
       mask[i] |= molecule_group_bit;
+      nmolcoords++;
     } else {
       mask[i] &= molecule_group_inversebit;
     }
   }
 
+  if (nmolcoords > nmaxmolcoords) {
+    nmaxmolcoords = nmolcoords;
+    molcoords = memory->grow(molcoords,nmaxmolcoords,3,"gcmc:molcoords");
+  }
+  
   double com[3];
   com[0] = com[1] = com[2] = 0.0;
   group->xcm(molecule_group,gas_mass,com);
@@ -1156,13 +1190,13 @@ void FixGCMC::attempt_molecule_rotation()
       xtmp[0] -= com[0];
       xtmp[1] -= com[1];
       xtmp[2] -= com[2];
-      MathExtra::matvec(rotmat,xtmp,atom_coord[n]);
+      MathExtra::matvec(rotmat,xtmp,molcoords[n]);
-      atom_coord[n][0] += com[0];
+      molcoords[n][0] += com[0];
-      atom_coord[n][1] += com[1];
+      molcoords[n][1] += com[1];
-      atom_coord[n][2] += com[2];
+      molcoords[n][2] += com[2];
-      xtmp[0] = atom_coord[n][0];
+      xtmp[0] = molcoords[n][0];
-      xtmp[1] = atom_coord[n][1];
+      xtmp[1] = molcoords[n][1];
-      xtmp[2] = atom_coord[n][2];
+      xtmp[2] = molcoords[n][2];
       domain->remap(xtmp);
       if (!domain->inside(xtmp))
         error->one(FLERR,"Fix gcmc put atom outside box");
@@ -1181,9 +1215,9 @@ void FixGCMC::attempt_molecule_rotation()
     for (int i = 0; i < atom->nlocal; i++) {
       if (mask[i] & molecule_group_bit) {
         image[i] = imagezero;
-        x[i][0] = atom_coord[n][0];
+        x[i][0] = molcoords[n][0];
-        x[i][1] = atom_coord[n][1];
+        x[i][1] = molcoords[n][1];
-        x[i][2] = atom_coord[n][2];
+        x[i][2] = molcoords[n][2];
         domain->remap(x[i],image[i]);
         n++;
       }
@@ -1303,18 +1337,18 @@ void FixGCMC::attempt_molecule_insertion()
   bool procflag[natoms_per_molecule];
 
   for (int i = 0; i < natoms_per_molecule; i++) {
-    MathExtra::matvec(rotmat,onemols[imol]->x[i],atom_coord[i]);
+    MathExtra::matvec(rotmat,onemols[imol]->x[i],molcoords[i]);
-    atom_coord[i][0] += com_coord[0];
+    molcoords[i][0] += com_coord[0];
-    atom_coord[i][1] += com_coord[1];
+    molcoords[i][1] += com_coord[1];
-    atom_coord[i][2] += com_coord[2];
+    molcoords[i][2] += com_coord[2];
 
     // use temporary variable for remapped position
-    // so unmapped position is preserved in atom_coord
+    // so unmapped position is preserved in molcoords
 
     double xtmp[3];
-    xtmp[0] = atom_coord[i][0];
+    xtmp[0] = molcoords[i][0];
-    xtmp[1] = atom_coord[i][1];
+    xtmp[1] = molcoords[i][1];
-    xtmp[2] = atom_coord[i][2];
+    xtmp[2] = molcoords[i][2];
     domain->remap(xtmp);
     if (!domain->inside(xtmp))
       error->one(FLERR,"Fix gcmc put atom outside box");
@@ -1372,7 +1406,7 @@ void FixGCMC::attempt_molecule_insertion()
 
     for (int i = 0; i < natoms_per_molecule; i++) {
       if (procflag[i]) {
-        atom->avec->create_atom(onemols[imol]->type[i],atom_coord[i]);
+        atom->avec->create_atom(onemols[imol]->type[i],molcoords[i]);
         int m = atom->nlocal - 1;
 
         // add to groups
@@ -1772,14 +1806,21 @@ void FixGCMC::attempt_molecule_rotation_full()
   double energy_before = energy_stored;
 
   int *mask = atom->mask;
+  int nmolcoords = 0;
   for (int i = 0; i < atom->nlocal; i++) {
     if (atom->molecule[i] == rotation_molecule) {
       mask[i] |= molecule_group_bit;
+      nmolcoords++;
     } else {
       mask[i] &= molecule_group_inversebit;
     }
   }
 
+  if (nmolcoords > nmaxmolcoords) {
+    nmaxmolcoords = nmolcoords;
+    molcoords = memory->grow(molcoords,nmaxmolcoords,3,"gcmc:molcoords");
+  }
+  
   double com[3];
   com[0] = com[1] = com[2] = 0.0;
   group->xcm(molecule_group,gas_mass,com);
@@ -1807,9 +1848,9 @@ void FixGCMC::attempt_molecule_rotation_full()
   int n = 0;
   for (int i = 0; i < atom->nlocal; i++) {
     if (mask[i] & molecule_group_bit) {
-      atom_coord[n][0] = x[i][0];
+      molcoords[n][0] = x[i][0];
-      atom_coord[n][1] = x[i][1];
+      molcoords[n][1] = x[i][1];
-      atom_coord[n][2] = x[i][2];
+      molcoords[n][2] = x[i][2];
       image_orig[n] = image[i];
       double xtmp[3];
       domain->unmap(x[i],image[i],xtmp);
@@ -1840,9 +1881,9 @@ void FixGCMC::attempt_molecule_rotation_full()
     int n = 0;
     for (int i = 0; i < atom->nlocal; i++) {
       if (mask[i] & molecule_group_bit) {
-        x[i][0] = atom_coord[n][0];
+        x[i][0] = molcoords[n][0];
-        x[i][1] = atom_coord[n][1];
+        x[i][1] = molcoords[n][1];
-        x[i][2] = atom_coord[n][2];
+        x[i][2] = molcoords[n][2];
         image[i] = image_orig[n];
         n++;
       }
@@ -2140,7 +2181,7 @@ double FixGCMC::energy(int i, int itype, tagint imolecule, double *coord)
   for (int j = 0; j < nall; j++) {
 
     if (i == j) continue;
-    if (mode == MOLECULE)
+    if (exchmode == EXCHMOL || movemode == MOVEMOL)
       if (imolecule == molecule[j]) continue;
 
     delx = coord[0] - x[j][0];
@@ -2212,9 +2253,10 @@ double FixGCMC::energy_full()
     tagint *molecule = atom->molecule;
     int nall = atom->nlocal + atom->nghost;
     for (int i = 0; i < atom->nlocal; i++) {
-      if (mode == MOLECULE) imolecule = molecule[i];
+      if (exchmode == EXCHMOL || movemode == MOVEMOL)
+        imolecule = molecule[i];
       for (int j = i+1; j < nall; j++) {
-        if (mode == MOLECULE)
+        if (exchmode == EXCHMOL || movemode == MOVEMOL)
           if (imolecule == molecule[j]) continue;
 
         delx = x[i][0] - x[j][0];
@@ -2352,7 +2394,7 @@ void FixGCMC::update_gas_atoms_list()
 
   if (regionflag) {
 
-    if (mode == MOLECULE) {
+    if (exchmode == EXCHMOL || movemode == MOVEMOL) {
 
       tagint maxmol = 0;
       for (int i = 0; i < nlocal; i++) maxmol = MAX(maxmol,molecule[i]);

src/MC/fix_gcmc.h

@@ -64,10 +64,12 @@ class FixGCMC : public Fix {
   int exclusion_group,exclusion_group_bit;
   int ngcmc_type,nevery,seed;
   int ncycles,nexchanges,nmcmoves;
+  double patomtrans, pmoltrans, pmolrotate, pmctot;
   int ngas;                 // # of gas atoms on all procs
   int ngas_local;           // # of gas atoms on this proc
   int ngas_before;          // # of gas atoms on procs < this proc
-  int mode;                 // ATOM or MOLECULE
+  int exchmode;             // exchange ATOM or MOLECULE
+  int movemode;             // move ATOM or MOLECULE
   int regionflag;           // 0 = anywhere in box, 1 = specific region
   int iregion;              // gcmc region
   char *idregion;           // gcmc region id
@@ -75,7 +77,8 @@ class FixGCMC : public Fix {
   bool charge_flag;         // true if user specified atomic charge
   bool full_flag;           // true if doing full system energy calculations
 
-  int natoms_per_molecule;  // number of atoms in each gas molecule
+  int natoms_per_molecule;  // number of atoms in each inserted molecule
+  int nmaxmolcoords;        // number of atoms allocated for molcoords
 
   int groupbitall;          // group bitmask for inserted atoms
   int ngroups;              // number of group-ids for inserted atoms
@@ -110,7 +113,7 @@ class FixGCMC : public Fix {
   double *sublo,*subhi;
   int *local_gas_list;
   double **cutsq;
-  double **atom_coord;
+  double **molcoords;
   imageint imagezero;
   double overlap_cutoffsq; // square distance cutoff for overlap 
   int overlap_flag;
@@ -126,8 +129,6 @@ class FixGCMC : public Fix {
 
   class Molecule **onemols;
   int imol,nmol;
-  double **coords;
-  imageint *imageflags;
   class Fix *fixrigid, *fixshake;
   int rigidflag, shakeflag;
   char *idrigid, *idshake;