diff --git a/doc/src/fix_gcmc.txt b/doc/src/fix_gcmc.txt
index c78142a99991a47fa5033f1f384f591effc312c0..38f0fb95ce0530d425689b0daa9052e6a6995074 100644
--- a/doc/src/fix_gcmc.txt
+++ b/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
diff --git a/examples/gcmc/in.gcmc.co2 b/examples/gcmc/in.gcmc.co2
index 128f05b489885885fdb68c9953572d71e6d44f76..bb6916fc4847df10707f2712c10955d921878601 100644
--- a/examples/gcmc/in.gcmc.co2
+++ b/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)
diff --git a/examples/gcmc/in.gcmc.lj b/examples/gcmc/in.gcmc.lj
index 3fe78efb252d2274a03b812f0bd972db2666952c..a1c7c6eb1002c5bbb7af59e364be39fa9cbc27ef 100644
--- a/examples/gcmc/in.gcmc.lj
+++ b/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
diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp
index a50ca0fd21e7004ed27fa021deb172536fcbe33a..f6843eb0a3a420a0d18f0e15fe1036990c565a50 100644
--- a/src/MC/fix_gcmc.cpp
+++ b/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");
+}
diff --git a/src/MC/fix_gcmc.h b/src/MC/fix_gcmc.h
index fc6021efea32fdd6b01e7f9161577ec318ff4a96..77dda17fdaeccbf5cd16f63b3e86e83baa7c51da 100644
--- a/src/MC/fix_gcmc.h
+++ b/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.