From 05a0275791d304d45f84bd5331dbfe12b7cb16b5 Mon Sep 17 00:00:00 2001
From: pscrozi <pscrozi@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 23 Jan 2015 17:28:13 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@12986
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/MC/fix_gcmc.cpp | 1073 +++++++++++++++++++++++++++++++++++++------
 src/MC/fix_gcmc.h   |   35 +-
 2 files changed, 945 insertions(+), 163 deletions(-)

diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp
index 8f7329ef71..f68430337e 100644
--- a/src/MC/fix_gcmc.cpp
+++ b/src/MC/fix_gcmc.cpp
@@ -27,15 +27,25 @@
 #include "modify.h"
 #include "fix.h"
 #include "comm.h"
+#include "compute.h"
 #include "group.h"
 #include "domain.h"
 #include "region.h"
 #include "random_park.h"
 #include "force.h"
 #include "pair.h"
+#include "bond.h"
+#include "angle.h"
+#include "dihedral.h"
+#include "improper.h"
+#include "kspace.h"
+#include "math_extra.h"
 #include "math_const.h"
 #include "memory.h"
 #include "error.h"
+#include "thermo.h"
+#include "output.h"
+#include "neighbor.h"
 #include <iostream>
 
 using namespace std;
@@ -55,6 +65,8 @@ FixGCMC::FixGCMC(LAMMPS *lmp, int narg, char **arg) :
   if (atom->molecular == 2) 
     error->all(FLERR,"Fix gcmc does not (yet) work with atom_style template");
 
+  dynamic_group_allow = 1;
+    
   vector_flag = 1;
   size_vector = 8;
   global_freq = 1;
@@ -84,18 +96,6 @@ FixGCMC::FixGCMC(LAMMPS *lmp, int narg, char **arg) :
 
   options(narg-11,&arg[11]);
 
-  // only one GCMC fix may handle a molecule
-  if (mode == MOLECULE) {
-    for (int i = 0; i < modify->nfix; i++) {
-      if (modify->fix[i] == this) continue;
-      if (strcmp(modify->fix[i]->style,"gcmc") == 0) {
-         FixGCMC *f = (FixGCMC *) modify->fix[i];
-         if (f->mode == MOLECULE)
-            error->all(FLERR,"Only one fix gcmc with 'molecule yes' allowed");
-      }
-    }
-  }
-
   // random number generator, same for all procs
 
   random_equal = new RanPark(lmp,seed);
@@ -153,8 +153,7 @@ FixGCMC::FixGCMC(LAMMPS *lmp, int narg, char **arg) :
       error->all(FLERR,"Fix gcmc molecule must have coordinates");
     if (onemols[imol]->typeflag == 0)
       error->all(FLERR,"Fix gcmc molecule must have atom types");
-    if (ngcmc_type+onemols[imol]->ntypes <= 0 || 
-        ngcmc_type+onemols[imol]->ntypes > atom->ntypes)
+    if (ngcmc_type+onemols[imol]->ntypes <= 0 || ngcmc_type+onemols[imol]->ntypes > atom->ntypes)
       error->all(FLERR,"Invalid atom type in fix gcmc mol command");
 
     if (atom->molecular == 2 && onemols != atom->avec->onemols)
@@ -214,13 +213,19 @@ void FixGCMC::options(int narg, char **arg)
   iregion = -1; 
   region_volume = 0;
   max_region_attempts = 1000; 
-  rotation_group = 0;
-  rotation_groupbit = 0;
-  rotation_inversegroupbit = 0;
+  molecule_group = 0;
+  molecule_group_bit = 0;
+  molecule_group_inversebit = 0;
+  exclusion_group = 0;
+  exclusion_group_bit = 0;
+  exclusion_group_inversebit = 0;
   pressure_flag = false;
   pressure = 0.0;
   fugacity_coeff = 1.0;
   shakeflag = 0;
+  charge = 0.0;
+  charge_flag = false;
+  full_flag = false;
   idshake = NULL;
   
   int iarg = 0;
@@ -261,6 +266,22 @@ void FixGCMC::options(int narg, char **arg)
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix gcmc command");
       fugacity_coeff = force->numeric(FLERR,arg[iarg+1]);
       iarg += 2;
+    } else if (strcmp(arg[iarg],"charge") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix gcmc command");
+      charge = force->numeric(FLERR,arg[iarg+1]);
+      charge_flag = true;
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"shake") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix gcmc command");
+      int n = strlen(arg[iarg+1]) + 1;
+      delete [] idshake;
+      idshake = new char[n];
+      strcpy(idshake,arg[iarg+1]);
+      shakeflag = 1;
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"full_energy") == 0) {
+      full_flag = true;
+      iarg += 1;
     } else error->all(FLERR,"Illegal fix gcmc command");
   }
 }
@@ -273,12 +294,20 @@ FixGCMC::~FixGCMC()
   delete random_equal;
   delete random_unequal;
 
-  // remove rotation group this fix defined
-  // only do it if the rotation group exists and group itself exists
+  // remove groups this fix defined
+  // only do it if the groups exists and group itself exists
 
-  if (rotation_group && (strcmp(group->names[0],"all") == 0)) {
+  if (molecule_group && (strcmp(group->names[0],"all") == 0)) {
     char **group_arg = new char*[2];
-    group_arg[0] = group->names[rotation_group];
+    group_arg[0] = group->names[molecule_group];
+    group_arg[1] = (char *) "delete";
+    group->assign(2,group_arg);
+    delete [] group_arg;
+  } 
+  
+  if (exclusion_group && (strcmp(group->names[0],"all") == 0)) {
+    char **group_arg = new char*[2];
+    group_arg[0] = group->names[exclusion_group];
     group_arg[1] = (char *) "delete";
     group->assign(2,group_arg);
     delete [] group_arg;
@@ -286,10 +315,10 @@ FixGCMC::~FixGCMC()
 
   memory->destroy(local_gas_list);
   memory->destroy(atom_coord);
-
-  delete [] idshake;
   memory->destroy(coords);
   memory->destroy(imageflags);
+  
+  delete [] idshake;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -304,7 +333,28 @@ int FixGCMC::setmask()
 /* ---------------------------------------------------------------------- */
 
 void FixGCMC::init()
-{
+{ 
+  // decide whether to switch to the full_energy option
+
+  if (!full_flag) {
+    if ((force->kspace) || 
+        (force->pair == NULL) ||
+        (force->pair->single_enable == 0) ||
+        (force->pair_match("hybrid",0)) ||
+        (force->pair_match("eam",0)) ||
+        (domain->triclinic == 1)) {
+      full_flag = true;
+      if (comm->me == 0) 
+        error->warning(FLERR,"fix gcmc using full_energy option");
+    }
+  }
+  
+  if (full_flag) {
+    char *id_pe = (char *) "thermo_pe";
+    int ipe = modify->find_compute(id_pe);
+    c_pe = modify->compute[ipe];
+  }
+
   int *type = atom->type;
 
   if (mode == ATOM) {
@@ -362,17 +412,42 @@ void FixGCMC::init()
       error->all(FLERR,"Fix gcmc and fix shake not using "
                  "same molecule template ID");
   }
-               
-  if (force->pair->single_enable == 0)
-    error->all(FLERR,"Fix gcmc incompatible with given pair_style");
 
   if (domain->dimension == 2)
     error->all(FLERR,"Cannot use fix gcmc in a 2d simulation");
 
-  if (domain->triclinic == 1)
-    error->all(FLERR,"Cannot use fix gcmc with a triclinic box");
+  // create a new group for interaction exclusions
 
-  // create a new group for rotation molecules
+  if (full_flag) {
+    char **group_arg = new char*[4];
+    // create unique group name for atoms to be excluded
+    int len = strlen(id) + 30;
+    group_arg[0] = new char[len];
+    sprintf(group_arg[0],"FixGCMC:gcmc_exclusion_group:%s",id);
+    group_arg[1] = (char *) "subtract";
+    group_arg[2] = (char *) "all";
+    group_arg[3] = (char *) "all";
+    group->assign(4,group_arg);
+    exclusion_group = group->find(group_arg[0]);
+    if (exclusion_group == -1) 
+      error->all(FLERR,"Could not find fix gcmc exclusion group ID");
+    exclusion_group_bit = group->bitmask[exclusion_group];
+    exclusion_group_inversebit = exclusion_group_bit ^ ~0;
+    
+    // neighbor list exclusion setup
+    // turn off interactions between group all and the exclusion group
+    
+    int narg = 4;  
+    char **arg = new char*[narg];;
+    arg[0] = (char *) "exclude";
+    arg[1] = (char *) "group";
+    arg[2] = group_arg[0];
+    arg[3] = (char *) "all";
+    neighbor->modify_params(narg,arg);
+    delete [] group_arg[0];
+  }
+    
+  // create a new group for temporary use with selected molecules
 
   if (mode == MOLECULE) {
     char **group_arg = new char*[3];
@@ -385,11 +460,11 @@ void FixGCMC::init()
     sprintf(digits,"%d",ngcmc_type);
     group_arg[2] = digits;
     group->assign(3,group_arg);
-    rotation_group = group->find(group_arg[0]);
-    if (rotation_group == -1) 
+    molecule_group = group->find(group_arg[0]);
+    if (molecule_group == -1) 
       error->all(FLERR,"Could not find fix gcmc rotation group ID");
-    rotation_groupbit = group->bitmask[rotation_group];
-    rotation_inversegroupbit = rotation_groupbit ^ ~0;
+    molecule_group_bit = group->bitmask[molecule_group];
+    molecule_group_inversebit = molecule_group_bit ^ ~0;
     delete [] group_arg[0];
   }
     
@@ -474,31 +549,65 @@ void FixGCMC::pre_exchange()
   comm->borders();
   update_gas_atoms_list();
 
-  if (mode == MOLECULE) {
-    for (int i = 0; i < ncycles; i++) {
-      int random_int_fraction =
-        static_cast<int>(random_equal->uniform()*ncycles) + 1;
-      if (random_int_fraction <= nmcmoves) {
-        if (random_equal->uniform() < 0.5) attempt_molecule_translation();
-        else attempt_molecule_rotation();
-      } else {
-        if (random_equal->uniform() < 0.5) attempt_molecule_deletion();
-        else attempt_molecule_insertion();
+  if (full_flag) {
+    energy_stored = energy_full();
+
+    if (mode == MOLECULE) {
+      for (int i = 0; i < ncycles; i++) {
+        int random_int_fraction =
+          static_cast<int>(random_equal->uniform()*ncycles) + 1;
+        if (random_int_fraction <= nmcmoves) {
+          if (random_equal->uniform() < 0.5) attempt_molecule_translation_full();
+          else attempt_molecule_rotation_full();
+        } else {
+          if (random_equal->uniform() < 0.5) attempt_molecule_deletion_full();
+          else attempt_molecule_insertion_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();
+        }
       }
     }
+    domain->pbc();
+    comm->exchange();
+    atom->nghost = 0;
+    comm->borders();
+    
   } 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();
+    
+    if (mode == MOLECULE) {
+      for (int i = 0; i < ncycles; i++) {
+        int random_int_fraction =
+          static_cast<int>(random_equal->uniform()*ncycles) + 1;
+        if (random_int_fraction <= nmcmoves) {
+          if (random_equal->uniform() < 0.5) attempt_molecule_translation();
+          else attempt_molecule_rotation();
+        } else {
+          if (random_equal->uniform() < 0.5) attempt_molecule_deletion();
+          else attempt_molecule_insertion();
+        }
+      }
+    } 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();
+        }
       }
     }
   }
-
   next_reneighbor = update->ntimestep + nevery;
 }
 
@@ -520,19 +629,33 @@ void FixGCMC::attempt_atomic_translation()
     double rsq = 1.1;
     double rx,ry,rz;
     rx = ry = rz = 0.0;
+    double coord[3];
     while (rsq > 1.0) {
       rx = 2*random_unequal->uniform() - 1.0;
       ry = 2*random_unequal->uniform() - 1.0;
       rz = 2*random_unequal->uniform() - 1.0;
       rsq = rx*rx + ry*ry + rz*rz;
-    }
-    double coord[3];
+    }    
     coord[0] = x[i][0] + displace*rx;
     coord[1] = x[i][1] + displace*ry;
     coord[2] = x[i][2] + displace*rz;
+    if (regionflag) {
+      while (domain->regions[iregion]->match(coord[0],coord[1],coord[2]) == 0) {
+        while (rsq > 1.0) {
+          rx = 2*random_unequal->uniform() - 1.0;
+          ry = 2*random_unequal->uniform() - 1.0;
+          rz = 2*random_unequal->uniform() - 1.0;
+          rsq = rx*rx + ry*ry + rz*rz;
+        }    
+        coord[0] = x[i][0] + displace*rx;
+        coord[1] = x[i][1] + displace*ry;
+        coord[2] = x[i][2] + displace*rz;
+      }
+    }
+
     double energy_after = energy(i,ngcmc_type,-1,coord);
     if (random_unequal->uniform() < 
-        exp(-beta*(energy_after - energy_before))) {
+        exp(beta*(energy_before - energy_after))) {
       x[i][0] = coord[0];
       x[i][1] = coord[1];
       x[i][2] = coord[2];
@@ -685,9 +808,41 @@ void FixGCMC::attempt_molecule_translation()
   com_displace[0] = displace*rx;
   com_displace[1] = displace*ry;
   com_displace[2] = displace*rz;
+  
+  int nlocal = atom->nlocal;
+  if (regionflag) {
+    int *mask = atom->mask;
+    for (int i = 0; i < nlocal; i++) {
+      if (atom->molecule[i] == translation_molecule) {
+        mask[i] |= molecule_group_bit;
+      } else {
+        mask[i] &= molecule_group_inversebit;
+      }
+    }
+    double com[3];
+    com[0] = com[1] = com[2] = 0.0;
+    group->xcm(molecule_group,gas_mass,com);
+    coord[0] = com[0] + displace*rx;
+    coord[1] = com[1] + displace*ry;
+    coord[2] = com[2] + displace*rz;
+    while (domain->regions[iregion]->match(coord[0],coord[1],coord[2]) == 0) {
+      while (rsq > 1.0) {
+        rx = 2*random_equal->uniform() - 1.0;
+        ry = 2*random_equal->uniform() - 1.0;
+        rz = 2*random_equal->uniform() - 1.0;
+        rsq = rx*rx + ry*ry + rz*rz;
+      }    
+      coord[0] = com[0] + displace*rx;
+      coord[1] = com[1] + displace*ry;
+      coord[2] = com[2] + displace*rz;
+    }
+    com_displace[0] = displace*rx;
+    com_displace[1] = displace*ry;
+    com_displace[2] = displace*rz;
+  }
 
   double energy_after = 0.0;
-  for (int i = 0; i < atom->nlocal; i++) {
+  for (int i = 0; i < nlocal; i++) {
     if (atom->molecule[i] == translation_molecule) {
       coord[0] = x[i][0] + com_displace[0];
       coord[1] = x[i][1] + com_displace[1];
@@ -700,8 +855,8 @@ void FixGCMC::attempt_molecule_translation()
   MPI_Allreduce(&energy_after,&energy_after_sum,1,MPI_DOUBLE,MPI_SUM,world);
 
   if (random_equal->uniform() < 
-      exp(-beta*(energy_after_sum - energy_before_sum))) {
-    for (int i = 0; i < atom->nlocal; i++) {
+      exp(beta*(energy_before_sum - energy_after_sum))) {
+    for (int i = 0; i < nlocal; i++) {
       if (atom->molecule[i] == translation_molecule) {
         x[i][0] += com_displace[0];
         x[i][1] += com_displace[1];
@@ -736,36 +891,41 @@ void FixGCMC::attempt_molecule_rotation()
   int *mask = atom->mask;
   for (int i = 0; i < nlocal; i++) {
     if (atom->molecule[i] == rotation_molecule) {
-      mask[i] |= rotation_groupbit;
+      mask[i] |= molecule_group_bit;
     } else {
-      mask[i] &= rotation_inversegroupbit;
+      mask[i] &= molecule_group_inversebit;
     }
   }
 
   double com[3];
   com[0] = com[1] = com[2] = 0.0;
-  group->xcm(rotation_group,gas_mass,com);
+  group->xcm(molecule_group,gas_mass,com);
 
-  double rot[9];
-  get_rotation_matrix(max_rotation_angle,&rot[0]);
+  double r[3],rotmat[3][3],quat[4];
+  r[0] = random_equal->uniform() - 0.5;
+  r[1] = random_equal->uniform() - 0.5;
+  r[2] = random_equal->uniform() - 0.5;
+
+  double theta = random_equal->uniform() * max_rotation_angle;
+  MathExtra::norm3(r);
+  MathExtra::axisangle_to_quat(r,theta,quat);
+  MathExtra::quat_to_mat(quat,rotmat);
 
   double **x = atom->x;
   imageint *image = atom->image;
   double energy_after = 0.0;
   int n = 0;
   for (int i = 0; i < nlocal; i++) {
-    if (mask[i] & rotation_groupbit) {
+    if (mask[i] & molecule_group_bit) {
       double xtmp[3];
       domain->unmap(x[i],image[i],xtmp);
       xtmp[0] -= com[0];
       xtmp[1] -= com[1];
       xtmp[2] -= com[2];
-      atom_coord[n][0] = 
-        rot[0]*xtmp[0] + rot[1]*xtmp[1] + rot[2]*xtmp[2] + com[0];
-      atom_coord[n][1] = 
-        rot[3]*xtmp[0] + rot[4]*xtmp[1] + rot[5]*xtmp[2] + com[1];
-      atom_coord[n][2] = 
-        rot[6]*xtmp[0] + rot[7]*xtmp[1] + rot[8]*xtmp[2] + com[2];
+      MathExtra::matvec(rotmat,xtmp,atom_coord[n]);    
+      atom_coord[n][0] += com[0];
+      atom_coord[n][1] += com[1];
+      atom_coord[n][2] += com[2];
       xtmp[0] = atom_coord[n][0];
       xtmp[1] = atom_coord[n][1];
       xtmp[2] = atom_coord[n][2];
@@ -779,10 +939,10 @@ void FixGCMC::attempt_molecule_rotation()
   MPI_Allreduce(&energy_after,&energy_after_sum,1,MPI_DOUBLE,MPI_SUM,world);
 
   if (random_equal->uniform() < 
-      exp(-beta*(energy_after_sum - energy_before_sum))) {
+      exp(beta*(energy_before_sum - energy_after_sum))) {
     int n = 0;
     for (int i = 0; i < nlocal; i++) {
-      if (mask[i] & rotation_groupbit) {
+      if (mask[i] & molecule_group_bit) {
         image[i] = imagetmp;
         x[i][0] = atom_coord[n][0];
         x[i][1] = atom_coord[n][1];
@@ -865,19 +1025,24 @@ void FixGCMC::attempt_molecule_insertion()
     com_coord[2] = zlo + random_equal->uniform() * (zhi-zlo);
   }
 
-  double rot[9];
-  get_rotation_matrix(MY_2PI,&rot[0]);
-  
+  double r[3],rotmat[3][3],quat[4];
+  r[0] = random_equal->uniform() - 0.5;
+  r[1] = random_equal->uniform() - 0.5;
+  r[2] = random_equal->uniform() - 0.5;
+
+  double theta = random_equal->uniform() * MY_2PI;
+  MathExtra::norm3(r);
+  MathExtra::axisangle_to_quat(r,theta,quat);
+  MathExtra::quat_to_mat(quat,rotmat);
+    
   double insertion_energy = 0.0;
   bool procflag[natoms_per_molecule];
 
   for (int i = 0; i < natoms_per_molecule; i++) {
-    atom_coord[i][0] = rot[0]*onemols[imol]->x[i][0] + 
-      rot[1]*onemols[imol]->x[i][1] + rot[2]*onemols[imol]->x[i][2] + com_coord[0];
-    atom_coord[i][1] = rot[3]*onemols[imol]->x[i][0] + 
-      rot[4]*onemols[imol]->x[i][1] + rot[5]*onemols[imol]->x[i][2] + com_coord[1];
-    atom_coord[i][2] = rot[6]*onemols[imol]->x[i][0] + 
-      rot[7]*onemols[imol]->x[i][1] + rot[8]*onemols[imol]->x[i][2] + com_coord[2];
+    MathExtra::matvec(rotmat,onemols[imol]->x[i],atom_coord[i]);
+    atom_coord[i][0] += com_coord[0];
+    atom_coord[i][1] += com_coord[1];
+    atom_coord[i][2] += com_coord[2];
 
     double xtmp[3];
     xtmp[0] = atom_coord[i][0];
@@ -899,7 +1064,7 @@ void FixGCMC::attempt_molecule_insertion()
                 MPI_DOUBLE,MPI_SUM,world);
 
   if (random_equal->uniform() < zz*volume*natoms_per_molecule*
-      exp(-beta*insertion_energy_sum)/(ngas+1)) {  
+      exp(-beta*insertion_energy_sum)/(ngas + natoms_per_molecule)) {  
 
     tagint maxmol = 0;
     for (int i = 0; i < atom->nlocal; i++) maxmol = MAX(maxmol,atom->molecule[i]);
@@ -919,6 +1084,11 @@ void FixGCMC::attempt_molecule_insertion()
 
     int nlocalprev = atom->nlocal;
     
+    double vnew[3];
+    vnew[0] = random_unequal->gaussian()*sigma;
+    vnew[1] = random_unequal->gaussian()*sigma;
+    vnew[2] = random_unequal->gaussian()*sigma;
+    
     for (int i = 0; i < natoms_per_molecule; i++) {
       if (procflag[i]) {
         atom->avec->create_atom(ngcmc_type+onemols[imol]->type[i],atom_coord[i]);
@@ -930,9 +1100,9 @@ void FixGCMC::attempt_molecule_insertion()
         if (maxtag_all+i+1 >= MAXTAGINT)
           error->all(FLERR,"Fix gcmc ran out of available atom IDs");
         atom->tag[m] = maxtag_all + i + 1;
-        atom->v[m][0] = random_unequal->gaussian()*sigma;
-        atom->v[m][1] = random_unequal->gaussian()*sigma;
-        atom->v[m][2] = random_unequal->gaussian()*sigma;
+        atom->v[m][0] = vnew[0];
+        atom->v[m][1] = vnew[1];
+        atom->v[m][2] = vnew[2];
         
         atom->add_molecule_atom(onemols[imol],i,m,maxtag_all);
         for (int j = 0; j < nfix; j++)
@@ -942,7 +1112,7 @@ void FixGCMC::attempt_molecule_insertion()
     }
 
     if (shakeflag) 
-      fixshake->set_molecule(nlocalprev,maxtag_all,imol,NULL,NULL,NULL);
+      fixshake->set_molecule(nlocalprev,maxtag_all,imol,com_coord,vnew,quat);
 
     atom->natoms += natoms_per_molecule;
     if (atom->natoms < 0 || atom->natoms > MAXBIGINT)
@@ -951,7 +1121,7 @@ void FixGCMC::attempt_molecule_insertion()
     atom->nangles += onemols[imol]->nangles;
     atom->ndihedrals += onemols[imol]->ndihedrals;
     atom->nimpropers += onemols[imol]->nimpropers;
-    atom->map_init();
+    if (atom->map_style) atom->map_init();
     atom->nghost = 0;
     comm->borders();
     update_gas_atoms_list();
@@ -999,6 +1169,621 @@ double FixGCMC::energy(int i, int itype, tagint imolecule, double *coord)
   return total_energy;
 }
 
+/* ----------------------------------------------------------------------
+   compute the energy of the given gas molecule in its current position 
+   sum across all procs that own atoms of the given molecule
+------------------------------------------------------------------------- */
+
+double FixGCMC::molecule_energy(tagint gas_molecule_id)
+{
+  double mol_energy = 0.0;
+  for (int i = 0; i < atom->nlocal; i++)
+    if (atom->molecule[i] == gas_molecule_id) {
+      mol_energy += energy(i,atom->type[i],gas_molecule_id,atom->x[i]);
+    }
+
+  double mol_energy_sum = 0.0;
+  MPI_Allreduce(&mol_energy,&mol_energy_sum,1,MPI_DOUBLE,MPI_SUM,world);
+  
+  return mol_energy_sum;
+}
+
+/* ----------------------------------------------------------------------
+------------------------------------------------------------------------- */
+
+void FixGCMC::attempt_atomic_translation_full()
+{
+  ntranslation_attempts += 1.0;
+  
+  if (ngas == 0) return;
+
+  double energy_before = energy_stored;
+  
+  int i = pick_random_gas_atom();
+
+  double **x = atom->x;
+  double xtmp[3];
+    
+  if (i >= 0) {
+  
+    double rsq = 1.1;
+    double rx,ry,rz;
+    rx = ry = rz = 0.0;
+    double coord[3];
+    while (rsq > 1.0) {
+      rx = 2*random_unequal->uniform() - 1.0;
+      ry = 2*random_unequal->uniform() - 1.0;
+      rz = 2*random_unequal->uniform() - 1.0;
+      rsq = rx*rx + ry*ry + rz*rz;
+    }    
+    coord[0] = x[i][0] + displace*rx;
+    coord[1] = x[i][1] + displace*ry;
+    coord[2] = x[i][2] + displace*rz;
+    if (regionflag) {
+      while (domain->regions[iregion]->match(coord[0],coord[1],coord[2]) == 0) {
+        while (rsq > 1.0) {
+          rx = 2*random_unequal->uniform() - 1.0;
+          ry = 2*random_unequal->uniform() - 1.0;
+          rz = 2*random_unequal->uniform() - 1.0;
+          rsq = rx*rx + ry*ry + rz*rz;
+        }    
+        coord[0] = x[i][0] + displace*rx;
+        coord[1] = x[i][1] + displace*ry;
+        coord[2] = x[i][2] + displace*rz;
+      }
+    }
+    xtmp[0] = x[i][0];
+    xtmp[1] = x[i][1];
+    xtmp[2] = x[i][2];
+    x[i][0] = coord[0];
+    x[i][1] = coord[1];
+    x[i][2] = coord[2];
+  }
+  
+  double energy_after = energy_full();
+  
+  if (random_equal->uniform() < 
+      exp(beta*(energy_before - energy_after))) {
+    energy_stored = energy_after;
+    ntranslation_successes += 1.0;
+  } else {
+    if (i >= 0) {
+      x[i][0] = xtmp[0];
+      x[i][1] = xtmp[1];
+      x[i][2] = xtmp[2];
+    }
+    energy_stored = energy_before;
+  } 
+  update_gas_atoms_list();
+}
+
+/* ----------------------------------------------------------------------
+------------------------------------------------------------------------- */
+
+void FixGCMC::attempt_atomic_deletion_full()
+{
+  double q_tmp;
+  
+  ndeletion_attempts += 1.0;
+
+  if (ngas == 0) return;
+  
+  double energy_before = energy_stored;
+
+  int i = pick_random_gas_atom();
+
+  if (i >= 0) {
+    atom->mask[i] |= exclusion_group_bit;
+    if (atom->q_flag) {
+      q_tmp = atom->q[i];
+      atom->q[i] = 0.0;
+    }
+  }
+  double energy_after = energy_full();
+
+  if (random_equal->uniform() < 
+      ngas*exp(beta*(energy_before - energy_after))/(zz*volume)) {
+    if (i >= 0) {
+      atom->avec->copy(atom->nlocal-1,i,1);
+      atom->nlocal--;  
+    }
+    atom->natoms--;
+    if (atom->map_style) atom->map_init();    
+    ndeletion_successes += 1.0;   
+    energy_stored = energy_after;    
+  } else {
+    if (i >= 0) {
+      atom->mask[i] &= exclusion_group_inversebit;
+      if (atom->q_flag) atom->q[i] = q_tmp;
+    }
+    energy_stored = energy_before;
+  }
+  update_gas_atoms_list();
+}
+
+/* ----------------------------------------------------------------------
+------------------------------------------------------------------------- */
+
+void FixGCMC::attempt_atomic_insertion_full()
+{
+  ninsertion_attempts += 1.0;
+
+  double energy_before = energy_stored;
+  
+  double coord[3];
+  if (regionflag) {
+    int region_attempt = 0;
+    coord[0] = region_xlo + random_equal->uniform() * (region_xhi-region_xlo);
+    coord[1] = region_ylo + random_equal->uniform() * (region_yhi-region_ylo);
+    coord[2] = region_zlo + random_equal->uniform() * (region_zhi-region_zlo);
+    while (domain->regions[iregion]->match(coord[0],coord[1],coord[2]) == 0) {
+      coord[0] = region_xlo + random_equal->uniform() * (region_xhi-region_xlo);
+      coord[1] = region_ylo + random_equal->uniform() * (region_yhi-region_ylo);
+      coord[2] = region_zlo + random_equal->uniform() * (region_zhi-region_zlo);
+      region_attempt++;
+      if (region_attempt >= max_region_attempts) return;
+    }
+  } else {
+    coord[0] = xlo + random_equal->uniform() * (xhi-xlo);
+    coord[1] = ylo + random_equal->uniform() * (yhi-ylo);
+    coord[2] = zlo + random_equal->uniform() * (zhi-zlo);
+  }
+  
+  int proc_flag = 0;
+  if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
+      coord[1] >= sublo[1] && coord[1] < subhi[1] &&
+      coord[2] >= sublo[2] && coord[2] < subhi[2]) {
+    proc_flag = 1;
+    atom->avec->create_atom(ngcmc_type,coord);
+    int m = atom->nlocal - 1;
+    atom->mask[m] = 1 | groupbit;
+    atom->v[m][0] = random_unequal->gaussian()*sigma;
+    atom->v[m][1] = random_unequal->gaussian()*sigma;
+    atom->v[m][2] = random_unequal->gaussian()*sigma;
+    if (charge_flag) atom->q[m] = charge;
+    
+    int nfix = modify->nfix;
+    Fix **fix = modify->fix;
+    for (int j = 0; j < nfix; j++)
+      if (fix[j]->create_attribute) fix[j]->set_arrays(m);
+  }
+
+  atom->natoms++;
+  if (atom->tag_enable) {
+    atom->tag_extend();
+    if (atom->map_style) atom->map_init();
+  }
+  atom->nghost = 0;
+  comm->borders();
+    
+  double energy_after = energy_full();
+  
+  if (random_equal->uniform() <
+      zz*volume*exp(beta*(energy_before - energy_after))/(ngas+1)) {
+
+    ninsertion_successes += 1.0;
+    energy_stored = energy_after;
+  } else {
+    atom->natoms--;
+    if (proc_flag) atom->nlocal--;
+    energy_stored = energy_before;
+  }
+  update_gas_atoms_list();
+}
+
+/* ----------------------------------------------------------------------
+------------------------------------------------------------------------- */
+
+void FixGCMC::attempt_molecule_translation_full()
+{
+  ntranslation_attempts += 1.0;
+
+  if (ngas == 0) return;
+
+  tagint translation_molecule = pick_random_gas_molecule();
+  if (translation_molecule == -1) return;
+
+  double energy_before = energy_stored;
+  
+  double **x = atom->x;
+  double rx,ry,rz;
+  double com_displace[3],coord[3];
+  double rsq = 1.1;
+  while (rsq > 1.0) {
+    rx = 2*random_equal->uniform() - 1.0;
+    ry = 2*random_equal->uniform() - 1.0;
+    rz = 2*random_equal->uniform() - 1.0;
+    rsq = rx*rx + ry*ry + rz*rz;
+  }
+  com_displace[0] = displace*rx;
+  com_displace[1] = displace*ry;
+  com_displace[2] = displace*rz;
+  
+  int nlocal = atom->nlocal;
+  if (regionflag) {
+    int *mask = atom->mask;
+    for (int i = 0; i < nlocal; i++) {
+      if (atom->molecule[i] == translation_molecule) {
+        mask[i] |= molecule_group_bit;
+      } else {
+        mask[i] &= molecule_group_inversebit;
+      }
+    }
+    double com[3];
+    com[0] = com[1] = com[2] = 0.0;
+    group->xcm(molecule_group,gas_mass,com);
+    coord[0] = com[0] + displace*rx;
+    coord[1] = com[1] + displace*ry;
+    coord[2] = com[2] + displace*rz;
+    while (domain->regions[iregion]->match(coord[0],coord[1],coord[2]) == 0) {
+      while (rsq > 1.0) {
+        rx = 2*random_equal->uniform() - 1.0;
+        ry = 2*random_equal->uniform() - 1.0;
+        rz = 2*random_equal->uniform() - 1.0;
+        rsq = rx*rx + ry*ry + rz*rz;
+      }    
+      coord[0] = com[0] + displace*rx;
+      coord[1] = com[1] + displace*ry;
+      coord[2] = com[2] + displace*rz;
+    }
+    com_displace[0] = displace*rx;
+    com_displace[1] = displace*ry;
+    com_displace[2] = displace*rz;
+  }
+  
+  for (int i = 0; i < nlocal; i++) {
+    if (atom->molecule[i] == translation_molecule) {
+      x[i][0] += com_displace[0];
+      x[i][1] += com_displace[1];
+      x[i][2] += com_displace[2];
+    }
+  }
+
+  double energy_after = energy_full();
+
+  if (random_equal->uniform() < 
+      exp(beta*(energy_before - energy_after))) {
+    ntranslation_successes += 1.0;
+    energy_stored = energy_after;
+  } else {
+    energy_stored = energy_before;
+    for (int i = 0; i < nlocal; i++) {
+      if (atom->molecule[i] == translation_molecule) {
+        x[i][0] -= com_displace[0];
+        x[i][1] -= com_displace[1];
+        x[i][2] -= com_displace[2];
+      }
+    }
+  }
+  update_gas_atoms_list();
+}
+
+/* ----------------------------------------------------------------------
+------------------------------------------------------------------------- */
+
+void FixGCMC::attempt_molecule_rotation_full()
+{
+  nrotation_attempts += 1.0;
+
+  if (ngas == 0) return;
+  
+  tagint rotation_molecule = pick_random_gas_molecule();
+  if (rotation_molecule == -1) return;
+  
+  double energy_before = energy_stored;
+
+  int nlocal = atom->nlocal;
+  int *mask = atom->mask;
+  for (int i = 0; i < nlocal; i++) {
+    if (atom->molecule[i] == rotation_molecule) {
+      mask[i] |= molecule_group_bit;
+    } else {
+      mask[i] &= molecule_group_inversebit;
+    }
+  }
+
+  double com[3];
+  com[0] = com[1] = com[2] = 0.0;
+  group->xcm(molecule_group,gas_mass,com);
+
+  double r[3],rotmat[3][3],quat[4];
+  r[0] = random_equal->uniform() - 0.5;
+  r[1] = random_equal->uniform() - 0.5;
+  r[2] = random_equal->uniform() - 0.5;
+
+  double theta = random_equal->uniform() * max_rotation_angle;
+  MathExtra::norm3(r);
+  MathExtra::axisangle_to_quat(r,theta,quat);
+  MathExtra::quat_to_mat(quat,rotmat);
+
+  double **x = atom->x;
+  imageint *image = atom->image;
+  imageint image_orig[natoms_per_molecule];
+  int n = 0;
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & molecule_group_bit) {
+      atom_coord[n][0] = x[i][0];
+      atom_coord[n][1] = x[i][1];
+      atom_coord[n][2] = x[i][2]; 
+      image_orig[n] = image[i];      
+      double xtmp[3];
+      domain->unmap(x[i],image[i],xtmp);
+      xtmp[0] -= com[0];
+      xtmp[1] -= com[1];
+      xtmp[2] -= com[2];
+      MathExtra::matvec(rotmat,xtmp,x[i]); 
+      x[i][0] += com[0];
+      x[i][1] += com[1];
+      x[i][2] += com[2];
+      image[i] = imagetmp;
+      domain->remap(x[i],image[i]);
+      n++;
+    }
+  }
+  
+  double energy_after = energy_full();
+
+  if (random_equal->uniform() < 
+      exp(beta*(energy_before - energy_after))) {
+    nrotation_successes += 1.0;
+    energy_stored = energy_after;    
+  } else {
+    energy_stored = energy_before;
+    int n = 0;
+    for (int i = 0; i < nlocal; i++) {
+      if (mask[i] & molecule_group_bit) {
+        x[i][0] = atom_coord[n][0];
+        x[i][1] = atom_coord[n][1];
+        x[i][2] = atom_coord[n][2];
+        image[i] = image_orig[n];
+        n++;
+      }
+    }
+  }
+  update_gas_atoms_list();
+}
+
+/* ----------------------------------------------------------------------
+------------------------------------------------------------------------- */
+
+void FixGCMC::attempt_molecule_deletion_full()
+{
+  ndeletion_attempts += 1.0;
+
+  if (ngas == 0) return;
+  
+  tagint deletion_molecule = pick_random_gas_molecule();
+  if (deletion_molecule == -1) return;
+
+  double energy_before = energy_stored;
+  
+  int m = 0;
+  double q_tmp[natoms_per_molecule];
+  for (int i = 0; i < atom->nlocal; i++) {
+    if (atom->molecule[i] == deletion_molecule) {
+      atom->mask[i] |= exclusion_group_bit;
+      if (atom->q_flag) {
+        q_tmp[m] = atom->q[i];
+        m++;
+        atom->q[i] = 0.0;
+      }
+      
+    }
+  }
+  
+  double energy_after = energy_full();  
+
+  if (random_equal->uniform() < 
+      ngas*exp(beta*(energy_before - energy_after))/(zz*volume*natoms_per_molecule)) {
+    int i = 0;
+    while (i < atom->nlocal) {
+      if (atom->molecule[i] == deletion_molecule) {
+        atom->avec->copy(atom->nlocal-1,i,1);
+        atom->nlocal--;
+      } else i++;
+    }
+    atom->natoms -= natoms_per_molecule;
+    if (atom->map_style) atom->map_init();
+    ndeletion_successes += 1.0;
+    energy_stored = energy_after;
+  } else {
+    energy_stored = energy_before;
+    int m = 0;
+    for (int i = 0; i < atom->nlocal; i++) {
+      if (atom->molecule[i] == deletion_molecule) {
+        atom->mask[i] &= exclusion_group_inversebit;
+        if (atom->q_flag) {
+          atom->q[i] = q_tmp[m];
+          m++;
+        }
+      }    
+    }
+  }
+  update_gas_atoms_list();
+}
+
+/* ----------------------------------------------------------------------
+------------------------------------------------------------------------- */
+
+void FixGCMC::attempt_molecule_insertion_full()
+{
+  ninsertion_attempts += 1.0;
+  
+  double energy_before = energy_stored;
+  
+  tagint maxmol = 0;
+  for (int i = 0; i < atom->nlocal; i++) maxmol = MAX(maxmol,atom->molecule[i]);
+  tagint maxmol_all;
+  MPI_Allreduce(&maxmol,&maxmol_all,1,MPI_LMP_TAGINT,MPI_MAX,world);
+  maxmol_all++;
+  if (maxmol_all >= MAXTAGINT) 
+    error->all(FLERR,"Fix gcmc ran out of available molecule IDs");
+  int insertion_molecule = maxmol_all;
+    
+  tagint maxtag = 0;
+  for (int i = 0; i < atom->nlocal; i++) maxtag = MAX(maxtag,atom->tag[i]);
+  tagint maxtag_all;
+  MPI_Allreduce(&maxtag,&maxtag_all,1,MPI_LMP_TAGINT,MPI_MAX,world);
+
+  int nfix = modify->nfix;
+  Fix **fix = modify->fix;
+
+  int nlocalprev = atom->nlocal;
+  
+  double com_coord[3];
+  if (regionflag) {
+    int region_attempt = 0;
+    com_coord[0] = region_xlo + random_equal->uniform() * 
+      (region_xhi-region_xlo);
+    com_coord[1] = region_ylo + random_equal->uniform() * 
+      (region_yhi-region_ylo);
+    com_coord[2] = region_zlo + random_equal->uniform() * 
+      (region_zhi-region_zlo);
+    while (domain->regions[iregion]->match(com_coord[0],com_coord[1],
+                                           com_coord[2]) == 0) {
+      com_coord[0] = region_xlo + random_equal->uniform() * 
+        (region_xhi-region_xlo);
+      com_coord[1] = region_ylo + random_equal->uniform() * 
+        (region_yhi-region_ylo);
+      com_coord[2] = region_zlo + random_equal->uniform() * 
+        (region_zhi-region_zlo);
+      region_attempt++;
+      if (region_attempt >= max_region_attempts) return;
+    }
+  } else {
+    com_coord[0] = xlo + random_equal->uniform() * (xhi-xlo);
+    com_coord[1] = ylo + random_equal->uniform() * (yhi-ylo);
+    com_coord[2] = zlo + random_equal->uniform() * (zhi-zlo);
+  }
+  
+  double r[3],rotmat[3][3],quat[4];
+  r[0] = random_equal->uniform() - 0.5;
+  r[1] = random_equal->uniform() - 0.5;
+  r[2] = random_equal->uniform() - 0.5;
+
+  double theta = random_equal->uniform() * MY_2PI;
+  MathExtra::norm3(r);
+  MathExtra::axisangle_to_quat(r,theta,quat);
+  MathExtra::quat_to_mat(quat,rotmat);
+
+  double vnew[3];
+  vnew[0] = random_unequal->gaussian()*sigma;
+  vnew[1] = random_unequal->gaussian()*sigma;
+  vnew[2] = random_unequal->gaussian()*sigma;
+    
+  for (int i = 0; i < natoms_per_molecule; i++) {
+    double xtmp[3];
+    MathExtra::matvec(rotmat,onemols[imol]->x[i],xtmp);
+    xtmp[0] += com_coord[0];
+    xtmp[1] += com_coord[1];
+    xtmp[2] += com_coord[2];
+    
+    domain->remap(xtmp);
+
+    if (xtmp[0] >= sublo[0] && xtmp[0] < subhi[0] &&
+        xtmp[1] >= sublo[1] && xtmp[1] < subhi[1] &&
+        xtmp[2] >= sublo[2] && xtmp[2] < subhi[2]) {
+
+      atom->avec->create_atom(ngcmc_type+onemols[imol]->type[i],xtmp);
+      int m = atom->nlocal - 1;
+      atom->mask[m] = 1 | groupbit;
+      atom->image[m] = imagetmp;
+      domain->remap(atom->x[m],atom->image[m]);
+      atom->molecule[m] = insertion_molecule;
+      if (maxtag_all+i+1 >= MAXTAGINT)
+        error->all(FLERR,"Fix gcmc ran out of available atom IDs");
+      atom->tag[m] = maxtag_all + i + 1;
+      atom->v[m][0] = vnew[0];
+      atom->v[m][1] = vnew[1];
+      atom->v[m][2] = vnew[2];
+
+      atom->add_molecule_atom(onemols[imol],i,m,maxtag_all);
+
+      for (int j = 0; j < nfix; j++)
+        if (fix[j]->create_attribute) fix[j]->set_arrays(m);
+    }
+  }
+
+  if (shakeflag) 
+    fixshake->set_molecule(nlocalprev,maxtag_all,imol,com_coord,vnew,quat);
+
+  atom->natoms += natoms_per_molecule;
+  if (atom->natoms < 0 || atom->natoms > MAXBIGINT)
+    error->all(FLERR,"Too many total atoms");
+  atom->nbonds += onemols[imol]->nbonds;
+  atom->nangles += onemols[imol]->nangles;
+  atom->ndihedrals += onemols[imol]->ndihedrals;
+  atom->nimpropers += onemols[imol]->nimpropers;
+  if (atom->map_style) atom->map_init();
+  atom->nghost = 0;
+  comm->borders();
+  double energy_after = energy_full();
+
+  if (random_equal->uniform() < zz*volume*natoms_per_molecule*
+      exp(beta*(energy_before - energy_after)/(ngas + natoms_per_molecule))) {  
+   
+    ninsertion_successes += 1.0;
+    energy_stored = energy_after;
+
+  } else {
+
+    atom->nbonds -= onemols[imol]->nbonds;
+    atom->nangles -= onemols[imol]->nangles;
+    atom->ndihedrals -= onemols[imol]->ndihedrals;
+    atom->nimpropers -= onemols[imol]->nimpropers;
+    atom->natoms -= natoms_per_molecule;
+  
+    energy_stored = energy_before;
+    int i = 0;
+    while (i < atom->nlocal) {
+      if (atom->molecule[i] == insertion_molecule) {
+        atom->avec->copy(atom->nlocal-1,i,1);
+        atom->nlocal--;
+      } else i++;
+    }
+  }
+  update_gas_atoms_list();
+}
+
+/* ----------------------------------------------------------------------
+   compute system potential energy
+------------------------------------------------------------------------- */
+
+double FixGCMC::energy_full()
+{ 
+  if (domain->triclinic) domain->x2lamda(atom->nlocal);
+  domain->pbc();
+  comm->exchange();
+  comm->borders();
+  if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
+  if (modify->n_pre_neighbor) modify->pre_neighbor();
+  neighbor->build();
+  int eflag = 1;
+  int vflag = 0;
+  
+  if (modify->n_pre_force) modify->pre_force(vflag);
+  
+  if (force->pair) force->pair->compute(eflag,vflag);
+
+  if (atom->molecular) {
+    if (force->bond) force->bond->compute(eflag,vflag);
+    if (force->angle) force->angle->compute(eflag,vflag);
+    if (force->dihedral) force->dihedral->compute(eflag,vflag);
+    if (force->improper) force->improper->compute(eflag,vflag);
+  }
+
+  if (force->kspace) force->kspace->compute(eflag,vflag);
+ 
+  if (modify->n_post_force) modify->post_force(vflag);
+  if (modify->n_end_of_step) modify->end_of_step();
+ 
+  update->eflag_global = update->ntimestep;
+  double total_energy = c_pe->compute_scalar();
+  if (output->thermo->normflag) total_energy *= atom->natoms;
+
+  return total_energy;
+}
+
 /* ----------------------------------------------------------------------
 ------------------------------------------------------------------------- */
 
@@ -1036,63 +1821,18 @@ tagint FixGCMC::pick_random_gas_molecule()
   return gas_molecule_id_all;
 }
 
-/* ----------------------------------------------------------------------
-   compute the energy of the given gas molecule in its current position 
-   sum across all procs that own atoms of the given molecule
-------------------------------------------------------------------------- */
-
-double FixGCMC::molecule_energy(tagint gas_molecule_id)
-{
-  double mol_energy = 0.0;
-  for (int i = 0; i < atom->nlocal; i++)
-    if (atom->molecule[i] == gas_molecule_id) {
-      mol_energy += energy(i,atom->type[i],gas_molecule_id,atom->x[i]);
-    }
-
-  double mol_energy_sum = 0.0;
-  MPI_Allreduce(&mol_energy,&mol_energy_sum,1,MPI_DOUBLE,MPI_SUM,world);
-  
-  return mol_energy_sum;
-}
-
-/* ----------------------------------------------------------------------
-   compute a 3x3 rotation matrix using 3 random Euler angles, 
-   each with a random maximum value supplied by the caller
-------------------------------------------------------------------------- */
-
-void FixGCMC::get_rotation_matrix(double max_angle, double *rot)
-{
-  double angle_x = max_angle*random_equal->uniform();
-  double angle_y = max_angle*random_equal->uniform();
-  double angle_z = max_angle*random_equal->uniform();
-  
-  double a = cos(angle_x);
-  double b = sin(angle_x);
-  double c = cos(angle_y);
-  double d = sin(angle_y);
-  double e = cos(angle_z);
-  double f = sin(angle_z);
-  double ad = a*d;
-  double bd = b*d;
-
-  rot[0] = c*e;
-  rot[1] = -c*f;
-  rot[2] = -d;
-  rot[3] = -bd*e + a*f;
-  rot[4] = bd*f + a*e;
-  rot[5] = -b*c;
-  rot[6] = ad*e + b*f;
-  rot[7] = -ad*f + b*e;
-  rot[8] = a*c;
-}
-
 /* ----------------------------------------------------------------------
    update the list of gas atoms
 ------------------------------------------------------------------------- */
 
 void FixGCMC::update_gas_atoms_list()
 {
-  if (atom->nlocal > gcmc_nmax) {
+  int nlocal = atom->nlocal;
+  int *mask = atom->mask;
+  tagint *molecule = atom->molecule;
+  double **x = atom->x;
+      
+  if (nlocal > gcmc_nmax) {
     memory->sfree(local_gas_list);
     gcmc_nmax = atom->nmax;
     local_gas_list = (int *) memory->smalloc(gcmc_nmax*sizeof(int),
@@ -1100,19 +1840,58 @@ void FixGCMC::update_gas_atoms_list()
   }
 
   ngas_local = 0;
+  
   if (regionflag) {
-    for (int i = 0; i < atom->nlocal; i++) {
-      if (atom->mask[i] & groupbit) {
-        double **x = atom->x;
-        if (domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2]) == 1) {
-          local_gas_list[ngas_local] = i;
-          ngas_local++;
+  
+    if (mode == MOLECULE) {
+    
+      tagint maxmol = 0;
+      for (int i = 0; i < nlocal; i++) maxmol = MAX(maxmol,molecule[i]);
+      tagint maxmol_all;
+      MPI_Allreduce(&maxmol,&maxmol_all,1,MPI_LMP_TAGINT,MPI_MAX,world);
+      double comx[maxmol_all];
+      double comy[maxmol_all];
+      double comz[maxmol_all];
+      for (int imolecule = 0; imolecule < maxmol_all; imolecule++) {
+        for (int i = 0; i < nlocal; i++) {
+          if (molecule[i] == imolecule) {
+            mask[i] |= molecule_group_bit;
+          } else {
+            mask[i] &= molecule_group_inversebit;
+          }
+        }
+        double com[3];
+        com[0] = com[1] = com[2] = 0.0;
+        group->xcm(molecule_group,gas_mass,com);
+        comx[imolecule] = com[0];
+        comy[imolecule] = com[1];
+        comz[imolecule] = com[2];
+      }
+    
+      for (int i = 0; i < nlocal; i++) {
+        if (mask[i] & groupbit) {
+          if (domain->regions[iregion]->match(comx[molecule[i]],
+             comy[molecule[i]],comz[molecule[i]]) == 1) {
+            local_gas_list[ngas_local] = i;
+            ngas_local++;
+          }
+        }
+      }
+      
+    } else {
+      for (int i = 0; i < nlocal; i++) {
+        if (mask[i] & groupbit) {
+          if (domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2]) == 1) {
+            local_gas_list[ngas_local] = i;
+            ngas_local++;
+          }
         }
       }
     }
+
   } else {
-    for (int i = 0; i < atom->nlocal; i++) {
-      if (atom->mask[i] & groupbit) {
+    for (int i = 0; i < nlocal; i++) {
+      if (mask[i] & groupbit) {
         local_gas_list[ngas_local] = i;
         ngas_local++;
       }
diff --git a/src/MC/fix_gcmc.h b/src/MC/fix_gcmc.h
index c0e7b4280a..57ad0c4296 100644
--- a/src/MC/fix_gcmc.h
+++ b/src/MC/fix_gcmc.h
@@ -1,4 +1,4 @@
-/* -*- c++ -*- ----------------------------------------------------------
+/* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
@@ -40,10 +40,17 @@ class FixGCMC : public Fix {
   void attempt_molecule_deletion();
   void attempt_molecule_insertion();
   double energy(int, int, tagint, double *);
+  void attempt_atomic_translation_full();
+  void attempt_atomic_deletion_full();
+  void attempt_atomic_insertion_full();
+  void attempt_molecule_translation_full();
+  void attempt_molecule_rotation_full();
+  void attempt_molecule_deletion_full();
+  void attempt_molecule_insertion_full();
+  double energy_full();
   int pick_random_gas_atom();
   tagint pick_random_gas_molecule();
   double molecule_energy(tagint);
-  void get_rotation_matrix(double, double *);
   void update_gas_atoms_list();
   double compute_vector(int);
   double memory_usage();
@@ -51,8 +58,10 @@ class FixGCMC : public Fix {
   void restart(char *);
 
  private:
-  int rotation_group,rotation_groupbit;
-  int rotation_inversegroupbit;
+  int molecule_group,molecule_group_bit;
+  int molecule_group_inversebit;
+  int exclusion_group,exclusion_group_bit;
+  int exclusion_group_inversebit;
   int ngcmc_type,nevery,seed;
   int ncycles,nexchanges,nmcmoves;
   int ngas;                 // # of gas atoms on all procs
@@ -63,7 +72,8 @@ class FixGCMC : public Fix {
   int iregion;              // GCMC region
   char *idregion;           // GCMC region id
   bool pressure_flag;       // true if user specified reservoir pressure
-                            // else false
+  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
 
@@ -84,10 +94,11 @@ class FixGCMC : public Fix {
   double displace;
   double max_rotation_angle;
   double beta,zz,sigma,volume;
-  double pressure,fugacity_coeff;
+  double pressure,fugacity_coeff,charge;
   double xlo,xhi,ylo,yhi,zlo,zhi;
   double region_xlo,region_xhi,region_ylo,region_yhi,region_zlo,region_zhi;
   double region_volume;
+  double energy_stored;
   double *sublo,*subhi;
   int *local_gas_list;
   double **cutsq;
@@ -108,6 +119,8 @@ class FixGCMC : public Fix {
   class Fix *fixshake;
   int shakeflag;
   char *idshake;
+  
+  class Compute *c_pe;
 
   void options(int, char **);
 };
@@ -169,21 +182,11 @@ Should not choose the GCMC molecule feature if no molecules are being
 simulated. The general molecule flag is off, but GCMC's molecule flag
 is on.
 
-E: Fix gcmc incompatible with given pair_style
-
-Some pair_styles do not provide single-atom energies, which are needed
-by fix gcmc.
-
 E: Cannot use fix gcmc in a 2d simulation
 
 Fix gcmc is set up to run in 3d only. No 2d simulations with fix gcmc
 are allowed.
 
-E: Cannot use fix gcmc with a triclinic box
-
-Fix gcmc is set up to run with othogonal boxes only. Simulations with
-triclinic boxes and fix gcmc are not allowed.
-
 E: Could not find fix gcmc rotation group ID
 
 Self-explanatory.