diff --git a/doc/src/fix_bond_break.txt b/doc/src/fix_bond_break.txt index 59fea8f45b..5927ceb4e5 100644 --- a/doc/src/fix_bond_break.txt +++ b/doc/src/fix_bond_break.txt @@ -137,8 +137,8 @@ doc page for more info. [Related commands:] "fix bond/create"_fix_bond_create.html, "fix -bond/swap"_fix_bond_swap.html, "dump local"_dump.html, -"special_bonds"_special_bonds.html +bond/react"_fix_bond_react.html, "fix bond/swap"_fix_bond_swap.html, +"dump local"_dump.html, "special_bonds"_special_bonds.html [Default:] diff --git a/doc/src/fix_bond_create.txt b/doc/src/fix_bond_create.txt index 02655577fd..de1a9f93f1 100644 --- a/doc/src/fix_bond_create.txt +++ b/doc/src/fix_bond_create.txt @@ -232,8 +232,8 @@ doc page for more info. [Related commands:] "fix bond/break"_fix_bond_break.html, "fix -bond/swap"_fix_bond_swap.html, "dump local"_dump.html, -"special_bonds"_special_bonds.html +bond/react"_fix_bond_react.html, "fix bond/swap"_fix_bond_swap.html, +"dump local"_dump.html, "special_bonds"_special_bonds.html [Default:] diff --git a/doc/src/fix_bond_react.txt b/doc/src/fix_bond_react.txt index f806da3324..fc08a9b14e 100644 --- a/doc/src/fix_bond_react.txt +++ b/doc/src/fix_bond_react.txt @@ -24,11 +24,11 @@ common_keyword = {stabilization} :l {stabilization} values = {no} or {yes} {group-ID} {xmax} {no} = no reaction site stabilization {yes} = perform reaction site stabilization - {group-ID} = user-assigned ID for all non-reacting atoms (group created internally) + {group-ID} = user-assigned prefix for the dynamic group of non-reacting atoms {xmax} = xmax value that is used by an internally created "nve/limit"_fix_nve_limit.html integrator :pre react = mandatory argument indicating new reaction specification :l react-ID = user-assigned name for the reaction :l - react-group-ID = only atoms in this group are available for the reaction :l + react-group-ID = only atoms in this group are considered for the reaction :l Nevery = attempt reaction every this many steps :l Rmin = bonding pair atoms must be separated by more than Rmin to initiate reaction (distance units) :l Rmax = bonding pair atoms must be separated by less than Rmax to initiate reaction (distance units) :l @@ -51,7 +51,7 @@ react = mandatory argument indicating new reaction specification :l molecule mol1 pre_reacted_topology.txt molecule mol2 post_reacted_topology.txt -fix 5 all bond/react stabilization no react myrxn1 all 1 0 3.25 mol1 mol2 map_file.txt :pre +fix 5 all bond/react react myrxn1 all 1 0 3.25 mol1 mol2 map_file.txt :pre molecule mol1 pre_reacted_rxn1.txt molecule mol2 post_reacted_rxn1.txt @@ -60,7 +60,7 @@ molecule mol4 post_reacted_rxn2.txt fix 5 all bond/react stabilization yes nvt_grp .03 & react myrxn1 all 1 0 3.25 mol1 mol2 map_file_rxn1.txt prob 0.50 12345 & react myrxn2 all 1 0 2.75 mol3 mol4 map_file_rxn2.txt prob 0.25 12345 -fix 6 nvt_grp nvt temp 300 300 100 # set thermostat after bond/react :pre +fix 6 nvt_grp_REACT nvt temp 300 300 100 # set thermostat after bond/react :pre [Description:] @@ -102,19 +102,29 @@ 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, the fix bond/react +command creates a "dynamic group"_group.html that is initialized to +include all atoms. If the group-ID is that of an existing static +group, the group is used as the parent group of new, +internally-created dynamic group. In both cases, this new dynamic +group is named by appending '_REACT' to the group-ID, e.g. +nvt_grp_REACT. By specifying an existing group, you may thermostat +constant-topology parts of your system separately. The dynamic group +contains only non-reacting atoms at a given timestep, and therefore +should 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. -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 static 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/doc/src/group.txt b/doc/src/group.txt index 8472677372..256c4ed7f5 100644 --- a/doc/src/group.txt +++ b/doc/src/group.txt @@ -225,8 +225,7 @@ atomfile-style variable. The variable is evaluated and atoms whose per-atom values are 0.0, are removed from the dynamic group. If the {property} keyword is used, the per-atom property name must be a previously defined per-atom property. The per-atom property is evaluated and atoms whose -values are 0.0 are removed from the dynamic group, otherwise they -are added to the group. +values are 0.0 are removed from the dynamic group. The assignment of atoms to a dynamic group is done at the beginning of each run and on every timestep that is a multiple of {N}, which is the diff --git a/examples/USER/misc/bond_react/nylon,6-6_melt/in.large_nylon_melt b/examples/USER/misc/bond_react/nylon,6-6_melt/in.large_nylon_melt index f2dc506dde..3982da799d 100644 --- a/examples/USER/misc/bond_react/nylon,6-6_melt/in.large_nylon_melt +++ b/examples/USER/misc/bond_react/nylon,6-6_melt/in.large_nylon_melt @@ -36,7 +36,7 @@ fix myrxns all bond/react stabilization yes statted_grp .03 & react rxn2 all 1 0.0 5.0 mol3 mol4 rxn1_stp2_map # stable at 800K -fix 1 statted_grp nvt temp 800 800 100 +fix 1 statted_grp_REACT nvt temp 800 800 100 # in order to customize behavior of reacting atoms, # you can use the internally created 'bond_react_MASTER_group', like so: diff --git a/examples/USER/misc/bond_react/tiny_nylon/in.tiny_nylon b/examples/USER/misc/bond_react/tiny_nylon/in.tiny_nylon index 1f7e9c42b7..73cb2eb270 100644 --- a/examples/USER/misc/bond_react/tiny_nylon/in.tiny_nylon +++ b/examples/USER/misc/bond_react/tiny_nylon/in.tiny_nylon @@ -36,7 +36,7 @@ fix myrxns all bond/react stabilization yes statted_grp .03 & react rxn1 all 1 0.0 2.9 mol1 mol2 rxn1_stp1_map & react rxn2 all 1 0.0 5.0 mol3 mol4 rxn1_stp2_map -fix 1 statted_grp nvt temp 300 300 100 +fix 1 statted_grp_REACT nvt temp 300 300 100 fix 4 bond_react_MASTER_group temp/rescale 1 300 300 10 1 diff --git a/src/SPIN/fix_nve_spin.cpp b/src/SPIN/fix_nve_spin.cpp index 898eb95396..7d96aca6d4 100644 --- a/src/SPIN/fix_nve_spin.cpp +++ b/src/SPIN/fix_nve_spin.cpp @@ -60,8 +60,11 @@ static const char cite_fix_nve_spin[] = "dynamics and molecular dynamics},\n" "author={Tranchida, J and Plimpton, SJ and Thibaudeau, P and Thompson, AP},\n" "journal={Journal of Computational Physics},\n" + "volume={372},\n" + "pages={406-425},\n" "year={2018},\n" "publisher={Elsevier}\n" + "doi={10.1016/j.jcp.2018.06.042}\n" "}\n\n"; enum{NONE}; diff --git a/src/USER-MISC/fix_bond_react.cpp b/src/USER-MISC/fix_bond_react.cpp index 705bb8ff70..a7b85e92d3 100644 --- a/src/USER-MISC/fix_bond_react.cpp +++ b/src/USER-MISC/fix_bond_react.cpp @@ -78,6 +78,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : fix1 = NULL; fix2 = NULL; + fix3 = NULL; if (narg < 8) error->all(FLERR,"Illegal fix bond/react command: " "too few arguments"); @@ -279,6 +280,9 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : open(files[i]); onemol = atom->molecules[unreacted_mol[i]]; twomol = atom->molecules[reacted_mol[i]]; + if (onemol->natoms != twomol->natoms) + error->all(FLERR,"Post-reacted template must contain the same " + "number of atoms as the pre-reacted template"); get_molxspecials(); read(i); fclose(fp); @@ -354,6 +358,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 +433,15 @@ 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 (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,61 +472,126 @@ 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 - if (group->find(master_group) == -1) { - 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; - } - - // 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; - } + 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) { - // 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 *) id_fix3; + 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); + fix3 = 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, use as parent group + // also, rename dynamic exclude_group by appending '_REACT' + char *exclude_PARENT_group; + int n = strlen(exclude_group) + 1; + exclude_PARENT_group = new char[n]; + strcpy(exclude_PARENT_group,exclude_group); + n += strlen("_REACT"); + delete [] exclude_group; + exclude_group = new char[n]; + strcpy(exclude_group,exclude_PARENT_group); + strcat(exclude_group,"_REACT"); + 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; - } + delete [] exclude_PARENT_group; + + // on to statted_tags (system-wide thermostat) + // intialize per-atom statted_flags to 1 + // (only if not already initialized by restart) + if (fix3->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 +610,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 +1854,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 +1884,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..ca1f3bd20c 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; @@ -77,9 +78,9 @@ class FixBondReact : public Fix { class Molecule *onemol; // pre-reacted molecule template class Molecule *twomol; // post-reacted molecule template Fix *fix1; // nve/limit used to relax reaction sites - Fix *fix2; // properties/atom used to indicate 1) indicate relaxing atoms - // 2) system-wide thermostat - // 3) to which 'react' atom belongs + Fix *fix2; // properties/atom used to indicate 1) relaxing atoms + // 2) to which 'react' atom belongs + Fix *fix3; // property/atom used for system-wide thermostat class RanMars **random; class NeighList *list; @@ -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