diff --git a/doc/src/fix_bond_react.txt b/doc/src/fix_bond_react.txt index f806da3324..b5ac87327d 100644 --- a/doc/src/fix_bond_react.txt +++ b/doc/src/fix_bond_react.txt @@ -102,19 +102,27 @@ involved in any new reactions. The {xmax} value keyword should typically be set to the maximum distance that non-reacting atoms move during the simulation. -The group-ID set using the {stabilization} keyword should be a -previously unused group-ID. It cannot be specified as 'all'. The fix -bond/react command creates a "dynamic group"_group.html of this name -that includes all non-reacting atoms. This dynamic group-ID should -then be used by a subsequent system-wide time integrator such as nvt, -npt, or nve, as shown in the second example above. It is currently -necessary to place the time integration command after the fix -bond/react command due to the internal dynamic grouping performed by -fix bond/react. +The group-ID set using the {stabilization} keyword can be an existing +static group or a previously-unused group-ID. It cannot be specified +as 'all'. If the group-ID is previously unused, fix bond/react command +creates a "dynamic group"_group.html of this name that is initialized +to include all atoms. If the group-ID is that of an existing static +group, the group is converted into a dynamic group, whose atoms are +limited to those belonging to the original static group. In either +case, this dynamic group-ID should then be used by a subsequent +system-wide time integrator such as nvt, npt, or nve, as shown in the +second example above. The time integration command should be placed +after the fix bond/react command due to the internal dynamic grouping +performed by fix bond/react. By specifying an existing group, you may +thermostat non-reacting parts of your system separately. -NOTE: The internally created group currently applies to all atoms in -the system, i.e. you should generally not have a separate thermostat -which acts on the 'all' group. +NOTE: If the group-ID is an existing static group, react-group-IDs +should also be specified as this group, or a subset. + +NOTE: If the group-ID is previously unused, the internally created +group applies to all atoms in the system, i.e. you should generally +not have a separate thermostat which acts on the 'all' group, or any +other group. The following comments pertain to each {react} argument: diff --git a/src/USER-MISC/fix_bond_react.cpp b/src/USER-MISC/fix_bond_react.cpp index 705bb8ff70..1124d359ce 100644 --- a/src/USER-MISC/fix_bond_react.cpp +++ b/src/USER-MISC/fix_bond_react.cpp @@ -354,6 +354,9 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : id_fix1 = NULL; id_fix2 = NULL; + id_fix3 = NULL; + statted_id = NULL; + custom_exclude_flag = 0; } /* ---------------------------------------------------------------------- */ @@ -426,11 +429,17 @@ FixBondReact::~FixBondReact() // check nfix in case all fixes have already been deleted if (id_fix1 == NULL && modify->nfix) modify->delete_fix(id_fix1); delete [] id_fix1; + + if (custom_exclude_flag == 0) { + if (id_fix3 == NULL && modify->nfix) modify->delete_fix(id_fix3); + delete [] id_fix3; + } } if (id_fix2 == NULL && modify->nfix) modify->delete_fix(id_fix2); delete [] id_fix2; + delete [] statted_id; delete [] guess_branch; delete [] pioneer_count; } @@ -461,21 +470,20 @@ void FixBondReact::post_constructor() int ifix = modify->find_fix(id_fix2); if (ifix == -1) { - char **newarg = new char*[8]; + char **newarg = new char*[7]; newarg[0] = (char *) "bond_react_props_internal"; newarg[1] = (char *) "all"; // group ID is ignored newarg[2] = (char *) "property/atom"; newarg[3] = (char *) "i_limit_tags"; - newarg[4] = (char *) "i_statted_tags"; - newarg[5] = (char *) "i_react_tags"; - newarg[6] = (char *) "ghost"; - newarg[7] = (char *) "yes"; - modify->add_fix(8,newarg); - fix2 = modify->fix[modify->nfix-1]; + newarg[4] = (char *) "i_react_tags"; + newarg[5] = (char *) "ghost"; + newarg[6] = (char *) "yes"; + modify->add_fix(7,newarg); delete [] newarg; } // create master_group if not already existing + // NOTE: limit_tags and react_tags automaticaly intitialized to zero (unless read from restart) if (group->find(master_group) == -1) { group->find_or_create(master_group); char **newarg; @@ -489,33 +497,100 @@ void FixBondReact::post_constructor() delete [] newarg; } - // on to statted_tags (system-wide thermostat) - // intialize per-atom statted_flags to 1 - // (only if not already initialized by restart) - // NOTE: limit_tags and react_tags automaticaly intitialized to zero (unless read from restart) - if (fix2->restart_reset != 1) { - int flag; - int index = atom->find_custom("statted_tags",flag); - int *i_statted_tags = atom->ivector[index]; - - for (int i = 0; i < atom->nlocal; i++) - i_statted_tags[i] = 1; - } - if (stabilization_flag == 1) { - // create exclude_group if not already existing - if (group->find(exclude_group) == -1) { + int igroup = group->find(exclude_group); + // create exclude_group if not already existing, or use as parent group if static + if (igroup == -1 || group->dynamic[igroup] == 0) { + // create stabilization per-atom property + len = strlen("bond_react_stabilization_internal") + 1; + id_fix3 = new char[len]; + strcpy(id_fix3,"bond_react_stabilization_internal"); + + ifix = modify->find_fix(id_fix3); + if (ifix == -1) { + char **newarg = new char*[6]; + newarg[0] = (char *) "bond_react_stabilization_internal"; + newarg[1] = (char *) "all"; // group ID is ignored + newarg[2] = (char *) "property/atom"; + newarg[3] = (char *) "i_statted_tags"; + newarg[4] = (char *) "ghost"; + newarg[5] = (char *) "yes"; + modify->add_fix(6,newarg); + fix2 = modify->fix[modify->nfix-1]; + delete [] newarg; + } + + len = strlen("statted_tags") + 1; + statted_id = new char[len]; + strcpy(statted_id,"statted_tags"); + + // if static group exists, duplicate it, use duplicate as parent group + // original will be converted into dynamic per-atom property group + if (igroup != -1) { + char **newarg; + newarg = new char*[3]; + newarg[0] = (char *) "exclude_PARENT_group"; + newarg[1] = (char *) "union"; + newarg[2] = exclude_group; + group->assign(3,newarg); + delete [] newarg; + } + group->find_or_create(exclude_group); char **newarg; newarg = new char*[5]; newarg[0] = exclude_group; newarg[1] = (char *) "dynamic"; - newarg[2] = (char *) "all"; + if (igroup == -1) newarg[2] = (char *) "all"; + else newarg[2] = (char *) "exclude_PARENT_group"; newarg[3] = (char *) "property"; newarg[4] = (char *) "statted_tags"; group->assign(5,newarg); delete [] newarg; - } + + // on to statted_tags (system-wide thermostat) + // intialize per-atom statted_flags to 1 + // (only if not already initialized by restart) + if (fix2->restart_reset != 1) { + int flag; + int index = atom->find_custom("statted_tags",flag); + int *i_statted_tags = atom->ivector[index]; + + for (int i = 0; i < atom->nlocal; i++) + i_statted_tags[i] = 1; + } + } else { + // sleeping code, for future capabilities + custom_exclude_flag = 1; + // first we have to find correct fix group reference + int n = strlen("GROUP_") + strlen(exclude_group) + 1; + char *fix_group = new char[n]; + strcpy(fix_group,"GROUP_"); + strcat(fix_group,exclude_group); + int ifix = modify->find_fix(fix_group); + Fix *fix = modify->fix[ifix]; + delete [] fix_group; + + // this returns names of corresponding property + int unused; + char * idprop; + idprop = (char *) fix->extract("property",unused); + if (idprop == NULL) + error->all(FLERR,"Exclude group must be a per-atom property group"); + + len = strlen(idprop) + 1; + statted_id = new char[len]; + strcpy(statted_id,idprop); + + // intialize per-atom statted_tags to 1 + // need to correct for smooth restarts + //int flag; + //int index = atom->find_custom(statted_id,flag); + //int *i_statted_tags = atom->ivector[index]; + //for (int i = 0; i < atom->nlocal; i++) + // i_statted_tags[i] = 1; + } + // let's create a new nve/limit fix to limit newly reacted atoms len = strlen("bond_react_MASTER_nve_limit") + 1; @@ -534,40 +609,6 @@ void FixBondReact::post_constructor() fix1 = modify->fix[modify->nfix-1]; delete [] newarg; } - - } - - // currently must redefine dynamic groups so they are updated at proper time - // -> should double check as to why - - int must_redefine_groups = 1; - - if (must_redefine_groups) { - group->find_or_create(master_group); - char **newarg; - newarg = new char*[5]; - newarg[0] = master_group; - newarg[1] = (char *) "dynamic"; - newarg[2] = (char *) "all"; - newarg[3] = (char *) "property"; - newarg[4] = (char *) "limit_tags"; - group->assign(5,newarg); - delete [] newarg; - } - - if (stabilization_flag == 1) { - if (must_redefine_groups) { - group->find_or_create(exclude_group); - char **newarg; - newarg = new char*[5]; - newarg[0] = exclude_group; - newarg[1] = (char *) "dynamic"; - newarg[2] = (char *) "all"; - newarg[3] = (char *) "property"; - newarg[4] = (char *) "statted_tags"; - group->assign(5,newarg); - delete [] newarg; - } } } @@ -1812,7 +1853,7 @@ void FixBondReact::limit_bond(int limit_bond_mode) int index1 = atom->find_custom("limit_tags",flag); int *i_limit_tags = atom->ivector[index1]; - int index2 = atom->find_custom("statted_tags",flag); + int index2 = atom->find_custom(statted_id,flag); int *i_statted_tags = atom->ivector[index2]; int index3 = atom->find_custom("react_tags",flag); @@ -1842,7 +1883,7 @@ void FixBondReact::unlimit_bond() int index1 = atom->find_custom("limit_tags",flag); int *i_limit_tags = atom->ivector[index1]; - int index2 = atom->find_custom("statted_tags",flag); + int index2 = atom->find_custom(statted_id,flag); int *i_statted_tags = atom->ivector[index2]; int index3 = atom->find_custom("react_tags",flag); diff --git a/src/USER-MISC/fix_bond_react.h b/src/USER-MISC/fix_bond_react.h index b960ec8e73..472a02be1a 100644 --- a/src/USER-MISC/fix_bond_react.h +++ b/src/USER-MISC/fix_bond_react.h @@ -57,6 +57,7 @@ class FixBondReact : public Fix { double **cutsq,*fraction; tagint lastcheck; int stabilization_flag; + int custom_exclude_flag; int *stabilize_steps_flag; int *update_edges_flag; int status; @@ -88,6 +89,8 @@ class FixBondReact : public Fix { char *nve_limit_xmax; // indicates max distance allowed to move when relaxing char *id_fix1; // id of internally created fix nve/limit char *id_fix2; // id of internally created fix per-atom properties + char *id_fix3; // id of internally created 'stabilization group' per-atom property fix + char *statted_id; // name of 'stabilization group' per-atom property char *master_group; // group containing relaxing atoms from all fix rxns char *exclude_group; // group for system-wide thermostat