lijiext
/
lammps
mirror of https://github.com/lammps/lammps.git
1
0
Fork 1
Code Issues Packages Projects Releases Wiki Activity

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

This commit is contained in:
Aidan Thompson 2018-02-05 13:29:14 -07:00
parent fa4c7fc664
commit 1d403b2aa3
5 changed files with 66 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
doc/src/fix_gcmc.txt
View File

@ -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

8
examples/gcmc/in.gcmc.co2
View File

@ -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)

8
examples/gcmc/in.gcmc.lj
View File

@ -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

55
src/MC/fix_gcmc.cpp
View File

@ -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");
}

13
src/MC/fix_gcmc.h
View File

@ -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.