diff --git a/doc/src/fix_bond_react.txt b/doc/src/fix_bond_react.txt index 9c82e2f0ff..4eced39207 100644 --- a/doc/src/fix_bond_react.txt +++ b/doc/src/fix_bond_react.txt @@ -144,7 +144,7 @@ modified to match the post-reaction template. A bonding atom pair will be identified if several conditions are met. First, a pair of atoms I,J within the specified react-group-ID of type -itype and jtype must separated by a distance between {Rmin} and +itype and jtype must be separated by a distance between {Rmin} and {Rmax}. It is possible that multiple bonding atom pairs are identified: if the bonding atoms in the pre-reacted template are not 1-2, 1-3, or 1-4 neighbors, the closest bonding atom partner is set as @@ -213,9 +213,10 @@ mandatory keyword is 'equivalences' and the optional keywords are N {equivalences} = # of atoms N in the reaction molecule templates N {edgeIDs} = # of edge atoms N in the pre-reacted molecule template N {deleteIDs} = # of atoms N that are specified for deletion -N {customIDs} = # of atoms N that are specified for a custom update :pre +N {customIDs} = # of atoms N that are specified for a custom update +N {constraints} = # of specified reaction constraints N :pre -The body of the map file contains two mandatory sections and three +The body of the map file contains two mandatory sections and four optional sections. The first mandatory section begins with the keyword 'BondingIDs' and lists the atom IDs of the bonding atom pair in the pre-reacted molecule template. The second mandatory section begins @@ -232,7 +233,10 @@ Edges' and allows for forcing the update of a specific atom's atomic charge. The first column is the ID of an atom near the edge of the pre-reacted molecule template, and the value of the second column is either 'none' or 'charges.' Further details are provided in the -discussion of the 'update_edges' keyword. +discussion of the 'update_edges' keyword. The fourth optional section +begins with the keyword 'Constraints' and lists additional criteria +that must be satisfied in order for the reaction to occur. Currently, +there is one type of constraint available, as discussed below. A sample map file is given below: @@ -265,6 +269,16 @@ Equivalences :pre :line +Any number of additional constraints may be specified in the +Constraints section of the map file. Currently there is one type of +additional constraint, of type 'distance', whose syntax is as follows: + +distance {ID1} {ID2} {rmin} {rmax} :pre + +where 'distance' is the required keyword, {ID1} and {ID2} are +pre-reaction atom IDs, and these two atoms must be separated by a +distance between {rmin} and {rmax} for the reaction to occur. + Once a reaction site has been successfully identified, data structures within LAMMPS that store bond topology are updated to reflect the post-reacted molecule template. All force fields with fixed bonds, diff --git a/src/USER-MISC/fix_bond_react.cpp b/src/USER-MISC/fix_bond_react.cpp index 4b15f50d7d..b9c357f363 100644 --- a/src/USER-MISC/fix_bond_react.cpp +++ b/src/USER-MISC/fix_bond_react.cpp @@ -58,7 +58,8 @@ static const char cite_fix_bond_react[] = #define BIG 1.0e20 #define DELTA 16 #define MAXLINE 256 -#define MAXGUESS 20 +#define MAXGUESS 20 // max # of guesses allowed by superimpose algorithm +#define MAXCONARGS 5 // max # of arguments for any type of constraint // various statuses of superimpose algorithm: // ACCEPT: site successfully matched to pre-reacted template @@ -168,6 +169,8 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : memory->create(limit_duration,nreacts,"bond/react:limit_duration"); memory->create(stabilize_steps_flag,nreacts,"bond/react:stabilize_steps_flag"); memory->create(update_edges_flag,nreacts,"bond/react:update_edges_flag"); + memory->create(nconstraints,nreacts,"bond/react:nconstraints"); + memory->create(constraints,nreacts,MAXCONARGS,"bond/react:constraints"); memory->create(iatomtype,nreacts,"bond/react:iatomtype"); memory->create(jatomtype,nreacts,"bond/react:jatomtype"); memory->create(ibonding,nreacts,"bond/react:ibonding"); @@ -185,6 +188,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : max_rxn[i] = INT_MAX; stabilize_steps_flag[i] = 0; update_edges_flag[i] = 0; + nconstraints[i] = 0; // set default limit duration to 60 timesteps limit_duration[i] = 60; reaction_count[i] = 0; @@ -958,6 +962,7 @@ void FixBondReact::far_partner() delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; + domain->minimum_image(delx,dely,delz); // ghost location fix rsq = delx*delx + dely*dely + delz*delz; if (rsq >= cutsq[rxnID][1] || rsq <= cutsq[rxnID][0]) { @@ -1139,7 +1144,7 @@ void FixBondReact::superimpose_algorithm() } } - if (status == ACCEPT) { // reaction site found successfully! + if (status == ACCEPT && check_constraints()) { // reaction site found successfully! glove_ghostcheck(); } hang_catch++; @@ -1572,6 +1577,32 @@ void FixBondReact::ring_check() } } +/* ---------------------------------------------------------------------- +evaluate constraints: return 0 if any aren't satisfied +------------------------------------------------------------------------- */ + +int FixBondReact::check_constraints() +{ + tagint atom1,atom2; + double delx,dely,delz,rsq; + + double **x = atom->x; + + for (int i = 0; i < nconstraints[rxnID]; i++) { + if (constraints[rxnID][0] == 0) { // 'distance' type + atom1 = atom->map(glove[(int) constraints[rxnID][1]-1][1]); + atom2 = atom->map(glove[(int) constraints[rxnID][2]-1][1]); + delx = x[atom1][0] - x[atom2][0]; + dely = x[atom1][1] - x[atom2][1]; + delz = x[atom1][2] - x[atom2][2]; + domain->minimum_image(delx,dely,delz); // ghost location fix + rsq = delx*delx + dely*dely + delz*delz; + if (rsq < constraints[rxnID][3] || rsq > constraints[rxnID][4]) return 0; + } + } + return 1; +} + /* ---------------------------------------------------------------------- Get xspecials for current molecule templates ------------------------------------------------------------------------- */ @@ -2677,6 +2708,7 @@ void FixBondReact::read(int myrxn) else if (strstr(line,"equivalences")) sscanf(line,"%d",&nequivalent); else if (strstr(line,"customIDs")) sscanf(line,"%d",&ncustom); else if (strstr(line,"deleteIDs")) sscanf(line,"%d",&ndelete); + else if (strstr(line,"constraints")) sscanf(line,"%d",&nconstraints[myrxn]); else break; } @@ -2708,6 +2740,8 @@ void FixBondReact::read(int myrxn) CustomEdges(line, myrxn); } else if (strcmp(keyword,"DeleteIDs") == 0) { DeleteAtoms(line, myrxn); + } else if (strcmp(keyword,"Constraints") == 0) { + Constraints(line, myrxn); } else error->one(FLERR,"Unknown section in superimpose file"); parse_keyword(1,line,keyword); @@ -2783,6 +2817,27 @@ void FixBondReact::DeleteAtoms(char *line, int myrxn) } } +void FixBondReact::Constraints(char *line, int myrxn) +{ + double tmp[MAXCONARGS]; + int n = strlen("distance") + 1; + char *constraint_type = new char[n]; + for (int i = 0; i < nconstraints[myrxn]; i++) { + readline(line); + sscanf(line,"%s",constraint_type); + if (strcmp(constraint_type,"distance") == 0) { + constraints[myrxn][0] = 0; // 0 = 'distance' ...maybe use another enum eventually + sscanf(line,"%*s %lg %lg %lg %lg",&tmp[0],&tmp[1],&tmp[2],&tmp[3]); + constraints[myrxn][1] = tmp[0]; + constraints[myrxn][2] = tmp[1]; + constraints[myrxn][3] = tmp[2]*tmp[2]; // using square of distance + constraints[myrxn][4] = tmp[3]*tmp[3]; + } else + error->one(FLERR,"Illegal constraint type in 'Constraints' section of map file"); + } + delete [] constraint_type; +} + void FixBondReact::open(char *file) { fp = fopen(file,"r"); diff --git a/src/USER-MISC/fix_bond_react.h b/src/USER-MISC/fix_bond_react.h index fffdfe4369..74d53b8f21 100644 --- a/src/USER-MISC/fix_bond_react.h +++ b/src/USER-MISC/fix_bond_react.h @@ -61,6 +61,8 @@ class FixBondReact : public Fix { int custom_exclude_flag; int *stabilize_steps_flag; int *update_edges_flag; + int *nconstraints; + double **constraints; int status; int *groupbits; @@ -141,6 +143,7 @@ class FixBondReact : public Fix { void Equivalences(char *,int); void CustomEdges(char *,int); void DeleteAtoms(char *,int); + void Constraints(char *,int); void make_a_guess (); void neighbor_loop(); @@ -148,6 +151,7 @@ class FixBondReact : public Fix { void crosscheck_the_neighbor(); void inner_crosscheck_loop(); void ring_check(); + int check_constraints(); void open(char *); void readline(char *);