ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #1195 from jrgissing/bond/react-custom_edges

Browse Source 
Bond/react: charge update for custom edges
This commit is contained in:
Axel Kohlmeyer
2018-11-08 16:19:23 -05:00
committed by GitHub
parent a6e9b99295 84fcf01bed
commit 3681bc853d
3 changed files with 82 additions and 21 deletions
Download Patch File Download Diff File Expand all files Collapse all files

32
doc/src/fix_bond_react.txt
View File

@ -44,7 +44,8 @@ react = mandatory argument indicating new reaction specification :l
timesteps = number of timesteps to apply internally created nve/limit.html
{update_edges} value = {none} or {charges} :l
none = do not update topology near the edges of reaction templates
charges = update atomic charges of all atoms in reaction templates :pre
charges = update atomic charges of all atoms in reaction templates
custom = force the update of user-specified atomic charges :pre
:ule
[Examples:]
@ -201,23 +202,30 @@ A discussion of correctly handling this is also provided on the
The map file is a text document with the following format:
A map file has a header and a body. The header of map file the
contains one mandatory keyword and one optional keyword. The mandatory
keyword is 'equivalences' and the optional keyword is 'edgeIDs':
contains one mandatory keyword and two optional keywords. The mandatory
keyword is 'equivalences' and the optional keywords are 'edgeIDs' and
'customIDs':
N {equivalences} = # of atoms N in the reaction molecule templates
N {edgeIDs} = # of edge atoms N in the pre-reacted molecule template :pre
N {edgeIDs} = # of edge atoms N in the pre-reacted molecule template
N {customIDs} = # of atoms N that are specified for a custom update :pre
The body of the map file contains two mandatory sections and one
optional section. The first mandatory section begins with the keyword
The body of the map file contains two mandatory sections and two
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
with the keyword 'Equivalences' and lists a one-to-one correspondence
between atom IDs of the pre- and post-reacted templates. The first
column is an atom ID of the pre-reacted molecule template, and the
second column is the corresponding atom ID of the post-reacted
molecule template. The optional section begins with the keyword
molecule template. The first optional section begins with the keyword
'EdgeIDs' and lists the atom IDs of edge atoms in the pre-reacted
molecule template.
molecule template. The second optional section begins with the keyword
'Custom 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.
A sample map file is given below:
@ -282,7 +290,13 @@ The {update_edges} keyword can increase the number of atoms whose
atomic charges are updated, when the pre-reaction template contains
edge atoms. When the value is set to 'charges,' all atoms' atomic
charges are updated to those specified by the post-reaction template,
including atoms near the edge of reaction templates.
including atoms near the edge of reaction templates. When the value is
set to 'custom,' an additional section must be included in the map
file that specifies whether to update charges, on a per-atom basis.
The format of this section is detailed above. Listing a pre-reaction
atom ID with a value of 'charges' will force the update of the atom's
charge, even if it is near a template edge. Atoms not near a template
edge are unaffected by this setting.
In order to produce the most physical behavior, this 'reaction site
equilibration time' should be tuned to be as small as possible while

67
src/USER-MISC/fix_bond_react.cpp
View File

@ -255,7 +255,10 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
} else if (strcmp(arg[iarg],"update_edges") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/react command: "
"'update_edges' has too few arguments");
if (strcmp(arg[iarg+1],"charges") == 0) update_edges_flag[rxn] = 1;
if (strcmp(arg[iarg+1],"none") == 0) update_edges_flag[rxn] = 0;
else if (strcmp(arg[iarg+1],"charges") == 0) update_edges_flag[rxn] = 1;
else if (strcmp(arg[iarg+1],"custom") == 0) update_edges_flag[rxn] = 2;
else error->all(FLERR,"Illegal value for 'update_edges' keyword'");
iarg += 2;
} else error->all(FLERR,"Illegal fix bond/react command: unknown keyword");
}
@ -271,11 +274,16 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
memory->create(reverse_equiv,max_natoms,2,nreacts,"bond/react:reverse_equiv");
memory->create(edge,max_natoms,nreacts,"bond/react:edge");
memory->create(landlocked_atoms,max_natoms,nreacts,"bond/react:landlocked_atoms");
memory->create(custom_edges,max_natoms,nreacts,"bond/react:custom_edges");
for (int j = 0; j < nreacts; j++)
for (int i = 0; i < max_natoms; i++) edge[i][j] = 0;
for (int i = 0; i < max_natoms; i++) {
edge[i][j] = 0;
if (update_edges_flag[j] == 1) custom_edges[i][j] = 1;
else custom_edges[i][j] = 0;
}
// read all superimpose files afterward
// read all map files afterward
for (int i = 0; i < nreacts; i++) {
open(files[i]);
onemol = atom->molecules[unreacted_mol[i]];
@ -384,6 +392,7 @@ FixBondReact::~FixBondReact()
memory->destroy(edge);
memory->destroy(equivalences);
memory->destroy(reverse_equiv);
memory->destroy(custom_edges);
memory->destroy(nevery);
memory->destroy(cutsq);
@ -1854,8 +1863,11 @@ 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_id,flag);
int *i_statted_tags = atom->ivector[index2];
int *i_statted_tags;
if (stabilization_flag == 1) {
int index2 = atom->find_custom(statted_id,flag);
i_statted_tags = atom->ivector[index2];
}
int index3 = atom->find_custom("react_tags",flag);
int *i_react_tags = atom->ivector[index3];
@ -1863,7 +1875,7 @@ void FixBondReact::limit_bond(int limit_bond_mode)
for (int i = 0; i < temp_limit_num; i++) {
// update->ntimestep could be 0. so add 1 throughout
i_limit_tags[atom->map(temp_limit_glove[i])] = update->ntimestep + 1;
i_statted_tags[atom->map(temp_limit_glove[i])] = 0;
if (stabilization_flag == 1) i_statted_tags[atom->map(temp_limit_glove[i])] = 0;
i_react_tags[atom->map(temp_limit_glove[i])] = rxnID;
}
@ -1884,8 +1896,11 @@ 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_id,flag);
int *i_statted_tags = atom->ivector[index2];
int *i_statted_tags;
if (stabilization_flag == 1) {
int index2 = atom->find_custom(statted_id,flag);
i_statted_tags = atom->ivector[index2];
}
int index3 = atom->find_custom("react_tags",flag);
int *i_react_tags = atom->ivector[index3];
@ -1895,7 +1910,7 @@ void FixBondReact::unlimit_bond()
// first '1': indexing offset, second '1': for next step
if (i_limit_tags[i] != 0 && (update->ntimestep + 1 - i_limit_tags[i]) > limit_duration[i_react_tags[i]]) { // + 1
i_limit_tags[i] = 0;
i_statted_tags[i] = 1;
if (stabilization_flag == 1) i_statted_tags[i] = 1;
i_react_tags[i] = 0;
}
}
@ -2077,7 +2092,7 @@ void FixBondReact::update_everything()
twomol = atom->molecules[reacted_mol[rxnID]];
for (int j = 0; j < twomol->natoms; j++) {
int jj = equivalences[j][1][rxnID]-1;
if ((landlocked_atoms[j][rxnID] == 1 || update_edges_flag[rxnID] == 1) &&
if ((landlocked_atoms[j][rxnID] == 1 || custom_edges[jj][rxnID] == 1) &&
atom->map(update_mega_glove[jj+1][i]) >= 0 &&
atom->map(update_mega_glove[jj+1][i]) < nlocal) {
type[atom->map(update_mega_glove[jj+1][i])] = twomol->type[j];
@ -2520,6 +2535,7 @@ void FixBondReact::read(int myrxn)
if (strstr(line,"edgeIDs")) sscanf(line,"%d",&nedge);
else if (strstr(line,"equivalences")) sscanf(line,"%d",&nequivalent);
else if (strstr(line,"customIDs")) sscanf(line,"%d",&ncustom);
else break;
}
@ -2532,7 +2548,7 @@ void FixBondReact::read(int myrxn)
// loop over sections of superimpose file
int equivflag = 0, edgeflag = 0, bondflag = 0;
int equivflag = 0, edgeflag = 0, bondflag = 0, customedgesflag = 0;
while (strlen(keyword)) {
if (strcmp(keyword,"BondingIDs") == 0) {
bondflag = 1;
@ -2546,6 +2562,9 @@ void FixBondReact::read(int myrxn)
} else if (strcmp(keyword,"Equivalences") == 0) {
equivflag = 1;
Equivalences(line, myrxn);
} else if (strcmp(keyword,"Custom Edges") == 0) {
customedgesflag = 1;
CustomEdges(line, myrxn);
} else error->one(FLERR,"Unknown section in superimpose file");
parse_keyword(1,line,keyword);
@ -2555,6 +2574,12 @@ void FixBondReact::read(int myrxn)
// error check
if (bondflag == 0 || equivflag == 0)
error->all(FLERR,"Superimpose file missing BondingIDs or Equivalences section\n");
if (update_edges_flag[myrxn] == 2 && customedgesflag == 0)
error->all(FLERR,"Map file must have a Custom Edges section when using 'update_edges custom'\n");
if (update_edges_flag[myrxn] != 2 && customedgesflag == 1)
error->all(FLERR,"Specify 'update_edges custom' to include Custom Edges section in map file\n");
}
void FixBondReact::EdgeIDs(char *line, int myrxn)
@ -2585,6 +2610,26 @@ void FixBondReact::Equivalences(char *line, int myrxn)
}
}
void FixBondReact::CustomEdges(char *line, int myrxn)
{
// 0 for 'none', 1 for 'charges'
int tmp;
int n = MAX(strlen("none"),strlen("charges")) + 1;
char *edgemode = new char[n];
for (int i = 0; i < ncustom; i++) {
readline(line);
sscanf(line,"%d %s",&tmp,edgemode);
if (strcmp(edgemode,"none") == 0)
custom_edges[tmp-1][myrxn] = 0;
else if (strcmp(edgemode,"charges") == 0)
custom_edges[tmp-1][myrxn] = 1;
else
error->one(FLERR,"Illegal value in 'Custom Edges' section of map file");
}
delete [] edgemode;
}
void FixBondReact::open(char *file)
{
fp = fopen(file,"r");

4
src/USER-MISC/fix_bond_react.h
View File

@ -101,7 +101,7 @@ class FixBondReact : public Fix {
int *ibonding,*jbonding;
int *closeneigh; // indicates if bonding atoms of a rxn are 1-2, 1-3, or 1-4 neighbors
int nedge,nequivalent; // number of edge, equivalent atoms in mapping file
int nedge,nequivalent,ncustom; // number of edge, equivalent, custom atoms in mapping file
int attempted_rxn; // there was an attempt!
int *local_rxn_count;
int *ghostly_rxn_count;
@ -115,6 +115,7 @@ class FixBondReact : public Fix {
int ***equivalences; // relation between pre- and post-reacted templates
int ***reverse_equiv; // re-ordered equivalences
int **landlocked_atoms; // all atoms at least three bonds away from edge atoms
int **custom_edges; // atoms in molecule templates with incorrect valences
int **nxspecial,**onemol_nxspecial,**twomol_nxspecial; // full number of 1-4 neighbors
tagint **xspecial,**onemol_xspecial,**twomol_xspecial; // full 1-4 neighbor list
@ -136,6 +137,7 @@ class FixBondReact : public Fix {
void read(int);
void EdgeIDs(char *,int);
void Equivalences(char *,int);
void CustomEdges(char *,int);
void make_a_guess ();
void neighbor_loop();