Merge branch 'bond/react-delete_atoms' of https://github.com/jrgissing/lammps into collected-post-stable-patches

This commit is contained in:
Axel Kohlmeyer 2018-12-12 00:20:17 -05:00
commit 0ed4da0bf9
4 changed files with 127 additions and 11 deletions

View File

@ -69,8 +69,9 @@ Initiate complex covalent bonding (topology) changes. These topology
changes will be referred to as 'reactions' throughout this
documentation. Topology changes are defined in pre- and post-reaction
molecule templates and can include creation and deletion of bonds,
angles, dihedrals, impropers, bond-types, angle-types, dihedral-types,
atom-types, or atomic charges.
angles, dihedrals, impropers, bond types, angle types, dihedral types,
atom types, or atomic charges. In addition, reaction by-products or
other molecules can be identified and deleted.
Fix bond/react does not use quantum mechanical (eg. fix qmmm) or
pairwise bond-order potential (eg. Tersoff or AIREBO) methods to
@ -203,15 +204,16 @@ 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 two optional keywords. The
contains one mandatory keyword and three optional keywords. The
mandatory keyword is 'equivalences' and the optional keywords are
'edgeIDs' and 'customIDs':
'edgeIDs' and 'deleteIDs' and 'customIDs':
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
The body of the map file contains two mandatory sections and two
The body of the map file contains two mandatory sections and three
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
@ -222,10 +224,12 @@ second column is the corresponding atom ID of the post-reacted
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. 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
'DeleteIDs' and lists the atom IDs of pre-reaction template atoms to
delete. The third 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:
@ -309,7 +313,16 @@ edge are unaffected by this setting.
A few other considerations:
It may be beneficial to ensure reacting atoms are at a certain
Many reactions result in one or more atoms that are considered
unwanted by-products. Therefore, bond/react provides the option to
delete a user-specified set of atoms. These pre-reaction atoms are
identified in the map file. A deleted atom must still be included in
the post-reaction molecule template, in which it cannot be bonded to
an atom that is not deleted. In addition to deleting unwanted reaction
by-products, this feature can be used to remove specific topologies,
such as small rings, that may be otherwise indistinguishable.
Also, it may be beneficial to ensure reacting atoms are at a certain
temperature before being released to the overall thermostat. For this,
you can use the internally-created dynamic group named
"bond_react_MASTER_group." For example, adding the following command

View File

@ -506,6 +506,7 @@ deepskyblue
defgrad
deformable
del
deleteIDs
Dellago
delocalization
delocalized

View File

@ -275,12 +275,14 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
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");
memory->create(delete_atoms,max_natoms,nreacts,"bond/react:delete_atoms");
for (int j = 0; j < nreacts; j++)
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;
delete_atoms[i][j] = 0;
}
// read all map files afterward
@ -393,6 +395,7 @@ FixBondReact::~FixBondReact()
memory->destroy(equivalences);
memory->destroy(reverse_equiv);
memory->destroy(custom_edges);
memory->destroy(delete_atoms);
memory->destroy(nevery);
memory->destroy(cutsq);
@ -1650,6 +1653,18 @@ void FixBondReact::find_landlocked_atoms(int myrxn)
}
}
// additionally, if a deleted atom is bonded to an atom that is not deleted, bad
for (int i = 0; i < onemol->natoms; i++) {
if (delete_atoms[i][myrxn] == 1) {
int ii = reverse_equiv[i][1][myrxn] - 1;
for (int j = 0; j < twomol_nxspecial[ii][0]; j++) {
if (delete_atoms[equivalences[twomol_xspecial[ii][j]-1][1][myrxn]-1][myrxn] == 0) {
error->one(FLERR,"A deleted atom cannot remain bonded to an atom that is not deleted");
}
}
}
}
// also, if atoms change number of bonds, but aren't landlocked, that could be bad
if (me == 0)
for (int i = 0; i < twomol->natoms; i++) {
@ -2053,6 +2068,13 @@ void FixBondReact::update_everything()
tagint **bond_atom = atom->bond_atom;
int *num_bond = atom->num_bond;
// used when deleting atoms
int ndel,ndelone;
int *mark = new int[nlocal];
for (int i = 0; i < nlocal; i++) mark[i] = 0;
tagint *tag = atom->tag;
AtomVec *avec = atom->avec;
// update atom->nbonds, etc.
// TODO: correctly tally with 'newton off'
int delta_bonds = 0;
@ -2086,6 +2108,18 @@ void FixBondReact::update_everything()
}
}
// mark to-delete atoms
for (int i = 0; i < update_num_mega; i++) {
rxnID = update_mega_glove[0][i];
onemol = atom->molecules[unreacted_mol[rxnID]];
for (int j = 0; j < onemol->natoms; j++) {
int iatom = atom->map(update_mega_glove[j+1][i]);
if (delete_atoms[j][rxnID] == 1 && iatom >= 0 && iatom < nlocal) {
mark[iatom] = 1;
}
}
}
// update charges and types of landlocked atoms
for (int i = 0; i < update_num_mega; i++) {
rxnID = update_mega_glove[0][i];
@ -2486,6 +2520,59 @@ void FixBondReact::update_everything()
memory->destroy(update_mega_glove);
// delete atoms. taken from fix_evaporate. but don't think it needs to be in pre_exchange
// loop in reverse order to avoid copying marked atoms
ndel = ndelone = 0;
for (int i = atom->nlocal-1; i >= 0; i--) {
if (mark[i] == 1) {
avec->copy(atom->nlocal-1,i,1);
atom->nlocal--;
ndelone++;
if (atom->avec->bonds_allow) {
if (force->newton_bond) delta_bonds += atom->num_bond[i];
else {
for (int j = 0; j < atom->num_bond[i]; j++) {
if (tag[i] < atom->bond_atom[i][j]) delta_bonds++;
}
}
}
if (atom->avec->angles_allow) {
if (force->newton_bond) delta_angle += atom->num_angle[i];
else {
for (int j = 0; j < atom->num_angle[i]; j++) {
int m = atom->map(atom->angle_atom2[i][j]);
if (m >= 0 && m < nlocal) delta_angle++;
}
}
}
if (atom->avec->dihedrals_allow) {
if (force->newton_bond) delta_dihed += atom->num_dihedral[i];
else {
for (int j = 0; j < atom->num_dihedral[i]; j++) {
int m = atom->map(atom->dihedral_atom2[i][j]);
if (m >= 0 && m < nlocal) delta_dihed++;
}
}
}
if (atom->avec->impropers_allow) {
if (force->newton_bond) delta_imprp += atom->num_improper[i];
else {
for (int j = 0; j < atom->num_improper[i]; j++) {
int m = atom->map(atom->improper_atom2[i][j]);
if (m >= 0 && m < nlocal) delta_imprp;
}
}
}
}
}
delete [] mark;
MPI_Allreduce(&ndelone,&ndel,1,MPI_INT,MPI_SUM,world);
atom->natoms -= ndel;
// done deleting atoms
// something to think about: this could done much more concisely if
// all atom-level info (bond,angles, etc...) were kinda inherited from a common data struct --JG
@ -2536,6 +2623,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 if (strstr(line,"deleteIDs")) sscanf(line,"%d",&ndelete);
else break;
}
@ -2565,6 +2653,8 @@ void FixBondReact::read(int myrxn)
} else if (strcmp(keyword,"Custom Edges") == 0) {
customedgesflag = 1;
CustomEdges(line, myrxn);
} else if (strcmp(keyword,"DeleteIDs") == 0) {
DeleteAtoms(line, myrxn);
} else error->one(FLERR,"Unknown section in superimpose file");
parse_keyword(1,line,keyword);
@ -2630,6 +2720,16 @@ void FixBondReact::CustomEdges(char *line, int myrxn)
delete [] edgemode;
}
void FixBondReact::DeleteAtoms(char *line, int myrxn)
{
int tmp;
for (int i = 0; i < ndelete; i++) {
readline(line);
sscanf(line,"%d",&tmp);
delete_atoms[tmp-1][myrxn] = 1;
}
}
void FixBondReact::open(char *file)
{
fp = fopen(file,"r");

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,ncustom; // number of edge, equivalent, custom atoms in mapping file
int nedge,nequivalent,ncustom,ndelete; // number of edge, equivalent, custom atoms in mapping file
int attempted_rxn; // there was an attempt!
int *local_rxn_count;
int *ghostly_rxn_count;
@ -116,6 +116,7 @@ class FixBondReact : public Fix {
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 **delete_atoms; // atoms in pre-reacted templates to delete
int **nxspecial,**onemol_nxspecial,**twomol_nxspecial; // full number of 1-4 neighbors
tagint **xspecial,**onemol_xspecial,**twomol_xspecial; // full 1-4 neighbor list
@ -138,6 +139,7 @@ class FixBondReact : public Fix {
void EdgeIDs(char *,int);
void Equivalences(char *,int);
void CustomEdges(char *,int);
void DeleteAtoms(char *,int);
void make_a_guess ();
void neighbor_loop();