Merge branch 'master' into user-cgdna

This commit is contained in:
Oliver Henrich 2018-11-06 19:48:06 +00:00
commit 4e6253254c
9 changed files with 155 additions and 98 deletions

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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