forked from lijiext/lammps
Merge pull request #778 from athomps/fix_gcmc_segfault_fix
Fixed recent segfault in fix gcmc and added mcmoves keyword
This commit is contained in:
Steve Plimpton
GitHub
commit
553b3ff69a
|
|
@ -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
|
||||
(insertions or deletions) of gas atoms or molecules of the given type
|
||||
between the simulation cell and the imaginary reservoir. It also
|
||||
attempts a number of Monte Carlo moves (translations and molecule
|
||||
rotations) of gas of the given type within the simulation cell or
|
||||
region. The average number of attempted GCMC exchanges is X. The
|
||||
average number of attempted MC moves is M. M should typically be
|
||||
Every N timesteps the fix attempts both GCMC exchanges
|
||||
(insertions or deletions) and MC moves of gas atoms or molecules.
|
||||
On those timesteps, the average number of attempted GCMC exchanges is X,
|
||||
while the average number of attempted MC moves is M.
|
||||
For GCMC exchanges of either molecular or atomic gasses,
|
||||
these exchanges can be either deletions or insertions,
|
||||
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
|
||||
attempted with 50% probability. For MC moves of atomic gasses,
|
||||
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
|
||||
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
|
||||
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
|
||||
molecules other than the exchanged type, but MC deletions,
|
||||
translations, and rotations can be performed on any atom/molecule in
|
||||
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
|
||||
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
|
||||
bond/angle/etc and special neighbor information for the molecule can
|
||||
molecule. The coordinates, atom types, charges, etc., as well as any
|
||||
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
|
||||
the user-specified fix group. For molecular gasses, exchanged
|
||||
deleted atoms are any atoms that have been inserted or that already
|
||||
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
|
||||
gcmc command (which can also be "all").
|
||||
two groups: the default group "all" and the fix group
|
||||
(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
|
||||
potential energy of the entire simulated system. The total system
|
||||
energy before and after the proposed GCMC move is then used in the
|
||||
The {full_energy} option means that the fix calculates the total
|
||||
potential energy of the entire simulated system, instead of just
|
||||
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
|
||||
partial energies are computed to determine the difference in energy
|
||||
that would be caused by the proposed GCMC move.
|
||||
proposed change. By default, this option is off,
|
||||
in which case only
|
||||
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
|
||||
|
||||
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
|
||||
The {group} keyword adds 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.
|
||||
|
||||
|
|
|
|||
|
|
@ -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),
|
||||
coords(NULL), imageflags(NULL), fixrigid(NULL), fixshake(NULL), idrigid(NULL), idshake(NULL)
|
||||
grouptypebits(NULL), grouptypes(NULL), local_gas_list(NULL), molcoords(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");
|
||||
|
||||
|
|
@ -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
|
||||
|
||||
if (mode == ATOM) natoms_per_molecule = 1;
|
||||
// setup of array of coordinates for molecule insertion
|
||||
// also used by rotation moves for any molecule
|
||||
|
||||
if (exchmode == EXCHATOM) natoms_per_molecule = 1;
|
||||
else natoms_per_molecule = onemols[imol]->natoms;
|
||||
memory->create(coords,natoms_per_molecule,3,"gcmc:coords");
|
||||
memory->create(imageflags,natoms_per_molecule,"gcmc:imageflags");
|
||||
memory->create(atom_coord,natoms_per_molecule,3,"gcmc:atom_coord");
|
||||
nmaxmolcoords = natoms_per_molecule;
|
||||
memory->create(molcoords,nmaxmolcoords,3,"gcmc:molcoords");
|
||||
|
||||
// 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(coords);
|
||||
memory->destroy(imageflags);
|
||||
memory->destroy(molcoords);
|
||||
|
||||
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)) ||
|
||||
((mode == MOLECULE) && (!atom->tag_enable || !atom->map_style)))
|
||||
error->all(FLERR,
|
||||
"Fix gcmc molecule command requires that "
|
||||
"atoms have molecule attributes");
|
||||
if (exchmode == EXCHMOL || movemode == MOVEMOL)
|
||||
if (atom->molecule_flag == 0 || !atom->tag_enable || !atom->map_style)
|
||||
error->all(FLERR,
|
||||
"Fix gcmc molecule command requires that "
|
||||
"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,
|
||||
// otherwise just get the gas mass
|
||||
// get all of the needed molecule data if exchanging
|
||||
// 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++) {
|
||||
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();
|
||||
for (int i = 0; i < ncycles; i++) {
|
||||
int ixm = static_cast<int>(random_equal->uniform()*ncycles) + 1;
|
||||
if (ixm <= nmcmoves) {
|
||||
double xmcmove = random_equal->uniform();
|
||||
if (xmcmove < patomtrans) attempt_atomic_translation_full();
|
||||
else if (xmcmove < patomtrans+pmoltrans) attempt_molecule_translation_full();
|
||||
else attempt_molecule_rotation_full();
|
||||
} else {
|
||||
double xgcmc = random_equal->uniform();
|
||||
if (exchmode == EXCHATOM) {
|
||||
if (xgcmc < 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++) {
|
||||
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();
|
||||
for (int i = 0; i < ncycles; i++) {
|
||||
int ixm = static_cast<int>(random_equal->uniform()*ncycles) + 1;
|
||||
if (ixm <= nmcmoves) {
|
||||
double xmcmove = random_equal->uniform();
|
||||
if (xmcmove < patomtrans) attempt_atomic_translation();
|
||||
else if (xmcmove < patomtrans+pmoltrans) attempt_molecule_translation();
|
||||
else attempt_molecule_rotation();
|
||||
} else {
|
||||
double xgcmc = random_equal->uniform();
|
||||
if (exchmode == EXCHATOM) {
|
||||
if (xgcmc < 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]);
|
||||
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];
|
||||
MathExtra::matvec(rotmat,xtmp,molcoords[n]);
|
||||
molcoords[n][0] += com[0];
|
||||
molcoords[n][1] += com[1];
|
||||
molcoords[n][2] += com[2];
|
||||
xtmp[0] = molcoords[n][0];
|
||||
xtmp[1] = molcoords[n][1];
|
||||
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][1] = atom_coord[n][1];
|
||||
x[i][2] = atom_coord[n][2];
|
||||
x[i][0] = molcoords[n][0];
|
||||
x[i][1] = molcoords[n][1];
|
||||
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]);
|
||||
atom_coord[i][0] += com_coord[0];
|
||||
atom_coord[i][1] += com_coord[1];
|
||||
atom_coord[i][2] += com_coord[2];
|
||||
MathExtra::matvec(rotmat,onemols[imol]->x[i],molcoords[i]);
|
||||
molcoords[i][0] += com_coord[0];
|
||||
molcoords[i][1] += com_coord[1];
|
||||
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[1] = atom_coord[i][1];
|
||||
xtmp[2] = atom_coord[i][2];
|
||||
xtmp[0] = molcoords[i][0];
|
||||
xtmp[1] = molcoords[i][1];
|
||||
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];
|
||||
atom_coord[n][1] = x[i][1];
|
||||
atom_coord[n][2] = x[i][2];
|
||||
molcoords[n][0] = x[i][0];
|
||||
molcoords[n][1] = x[i][1];
|
||||
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][1] = atom_coord[n][1];
|
||||
x[i][2] = atom_coord[n][2];
|
||||
x[i][0] = molcoords[n][0];
|
||||
x[i][1] = molcoords[n][1];
|
||||
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]);
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
Loading…
Reference in New Issue