diff --git a/doc/src/fix_gcmc.txt b/doc/src/fix_gcmc.txt index 405738276f..c78142a999 100644 --- a/doc/src/fix_gcmc.txt +++ b/doc/src/fix_gcmc.txt @@ -26,16 +26,20 @@ zero or more keyword/value pairs may be appended to args :l keyword = {mol}, {region}, {maxangle}, {pressure}, {fugacity_coeff}, {full_energy}, {charge}, {group}, {grouptype}, {intra_energy}, {tfac_insert}, or {overlap_cutoff} {mol} value = template-ID template-ID = ID of molecule template specified in a separate "molecule"_molecule.html command + {mcmoves} values = Patomtrans Pmoltrans Pmolrotate + Patomtrans = proportion of atom translation MC moves + Pmoltrans = proportion of molecule translation MC moves + Pmolrotate = proportion of molecule rotation MC moves {rigid} value = fix-ID fix-ID = ID of "fix rigid/small"_fix_rigid.html command {shake} value = fix-ID fix-ID = ID of "fix shake"_fix_shake.html command {region} value = region-ID - region-ID = ID of region where MC moves are allowed + region-ID = ID of region where GCMC exchanges and MC moves are allowed {maxangle} value = maximum molecular rotation angle (degrees) {pressure} value = pressure of the gas reservoir (pressure units) {fugacity_coeff} value = fugacity coefficient of the gas reservoir (unitless) - {full_energy} = compute the entire system energy when performing MC moves + {full_energy} = compute the entire system energy when performing GCMC exchanges and MC moves {charge} value = charge of inserted atoms (charge units) {group} value = group-ID group-ID = group-ID for inserted atoms (string) @@ -56,34 +60,42 @@ fix 4 my_gas gcmc 1 10 10 1 123456543 300.0 -12.5 1.0 region disk :pre [Description:] This fix performs grand canonical Monte Carlo (GCMC) exchanges of -atoms or molecules of the given type with an imaginary ideal gas +atoms or molecules with an imaginary ideal gas reservoir at the specified T and chemical potential (mu) as discussed -in "(Frenkel)"_#Frenkel. If used with the "fix nvt"_fix_nh.html +in "(Frenkel)"_#Frenkel. It also +attempts Monte Carlo (MC) moves (translations and molecule +rotations) within the simulation cell or +region. If used with the "fix nvt"_fix_nh.html command, simulations in the grand canonical ensemble (muVT, constant chemical potential, constant volume, and constant temperature) can be performed. Specific uses include computing isotherms in microporous materials, or computing vapor-liquid coexistence curves. -Every N timesteps the fix attempts a number of GCMC exchanges -(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. diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp index 48ba8ed5cf..a50ca0fd21 100644 --- a/src/MC/fix_gcmc.cpp +++ b/src/MC/fix_gcmc.cpp @@ -64,15 +64,16 @@ using namespace MathConst; #define MAXENERGYTEST 1.0e50 -enum{ATOM,MOLECULE}; +enum{EXCHATOM,EXCHMOL}; // exchmode +enum{MOVEATOM,MOVEMOL}; // movemode /* ---------------------------------------------------------------------- */ FixGCMC::FixGCMC(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), idregion(NULL), full_flag(0), ngroups(0), groupstrings(NULL), ngrouptypes(0), grouptypestrings(NULL), - grouptypebits(NULL), grouptypes(NULL), local_gas_list(NULL), atom_coord(NULL), random_equal(NULL), random_unequal(NULL), - 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 (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(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(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(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(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(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(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]); diff --git a/src/MC/fix_gcmc.h b/src/MC/fix_gcmc.h index 3656a1df58..fc6021efea 100644 --- a/src/MC/fix_gcmc.h +++ b/src/MC/fix_gcmc.h @@ -64,10 +64,12 @@ class FixGCMC : public Fix { int exclusion_group,exclusion_group_bit; int ngcmc_type,nevery,seed; int ncycles,nexchanges,nmcmoves; + double patomtrans, pmoltrans, pmolrotate, pmctot; int ngas; // # of gas atoms on all procs int ngas_local; // # of gas atoms on this proc int ngas_before; // # of gas atoms on procs < this proc - int mode; // ATOM or MOLECULE + int exchmode; // exchange ATOM or MOLECULE + int movemode; // move ATOM or MOLECULE int regionflag; // 0 = anywhere in box, 1 = specific region int iregion; // gcmc region char *idregion; // gcmc region id @@ -75,7 +77,8 @@ class FixGCMC : public Fix { bool charge_flag; // true if user specified atomic charge bool full_flag; // true if doing full system energy calculations - int natoms_per_molecule; // number of atoms in each gas molecule + int natoms_per_molecule; // number of atoms in each inserted molecule + int nmaxmolcoords; // number of atoms allocated for molcoords int groupbitall; // group bitmask for inserted atoms int ngroups; // number of group-ids for inserted atoms @@ -110,7 +113,7 @@ class FixGCMC : public Fix { double *sublo,*subhi; int *local_gas_list; double **cutsq; - double **atom_coord; + double **molcoords; imageint imagezero; double overlap_cutoffsq; // square distance cutoff for overlap int overlap_flag; @@ -126,8 +129,6 @@ class FixGCMC : public Fix { class Molecule **onemols; int imol,nmol; - double **coords; - imageint *imageflags; class Fix *fixrigid, *fixshake; int rigidflag, shakeflag; char *idrigid, *idshake;