From 53e66dcd154d17c3887ba224ae1f0f6abfa97356 Mon Sep 17 00:00:00 2001 From: jrgissing Date: Fri, 30 Nov 2018 22:35:10 -0700 Subject: [PATCH 1/4] bond/react: delete atoms allows deleting of a user-specified set of atoms, based on topology --- doc/src/fix_bond_react.txt | 33 ++++++++---- src/USER-MISC/fix_bond_react.cpp | 88 ++++++++++++++++++++++++++++++++ src/USER-MISC/fix_bond_react.h | 4 +- 3 files changed, 114 insertions(+), 11 deletions(-) diff --git a/doc/src/fix_bond_react.txt b/doc/src/fix_bond_react.txt index b62bf68a75..e55e06b7d2 100644 --- a/doc/src/fix_bond_react.txt +++ b/doc/src/fix_bond_react.txt @@ -70,7 +70,8 @@ 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. +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 @@ -79,7 +80,7 @@ probabilistic criteria to effect predetermined topology changes in simulations using standard force fields. This fix was created to facilitate the dynamic creation of polymeric, -amorphous or highly cross-linked systems. A suggested workflow for +amorphous or highly-crosslinked systems. A suggested workflow for using this fix is: 1) identify a reaction to be simulated 2) build a molecule template of the reaction site before the reaction has occurred 3) build a molecule template of the reaction site after the @@ -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 undeleted atom. 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 diff --git a/src/USER-MISC/fix_bond_react.cpp b/src/USER-MISC/fix_bond_react.cpp index 394be64460..b06f40c3f3 100644 --- a/src/USER-MISC/fix_bond_react.cpp +++ b/src/USER-MISC/fix_bond_react.cpp @@ -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); @@ -2053,6 +2056,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 +2096,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 +2508,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 +2611,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 +2641,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 +2708,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"); diff --git a/src/USER-MISC/fix_bond_react.h b/src/USER-MISC/fix_bond_react.h index d54ab7c385..d6e7b785e7 100644 --- a/src/USER-MISC/fix_bond_react.h +++ b/src/USER-MISC/fix_bond_react.h @@ -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(); From d14404254e8dcd1ef1f2884463c073cd0412a870 Mon Sep 17 00:00:00 2001 From: jrgissing Date: Sat, 1 Dec 2018 15:18:11 -0700 Subject: [PATCH 2/4] check for illegally deleted atoms --- src/USER-MISC/fix_bond_react.cpp | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/USER-MISC/fix_bond_react.cpp b/src/USER-MISC/fix_bond_react.cpp index b06f40c3f3..50efca8d09 100644 --- a/src/USER-MISC/fix_bond_react.cpp +++ b/src/USER-MISC/fix_bond_react.cpp @@ -1653,6 +1653,18 @@ void FixBondReact::find_landlocked_atoms(int myrxn) } } + // additionally, if a deleted atom is bonded to a undeleted atom, 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 be bonded to an undeleted atom"); + } + } + } + } + // 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++) { From cf3aee908349692d84fb2e2ba59677aed8cece7f Mon Sep 17 00:00:00 2001 From: jrgissing Date: Mon, 3 Dec 2018 20:12:38 -0700 Subject: [PATCH 3/4] bond/react doc tweaks --- doc/src/fix_bond_react.txt | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/src/fix_bond_react.txt b/doc/src/fix_bond_react.txt index e55e06b7d2..5d6349a426 100644 --- a/doc/src/fix_bond_react.txt +++ b/doc/src/fix_bond_react.txt @@ -69,8 +69,8 @@ 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. In addition, reaction by-products or +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 @@ -80,7 +80,7 @@ probabilistic criteria to effect predetermined topology changes in simulations using standard force fields. This fix was created to facilitate the dynamic creation of polymeric, -amorphous or highly-crosslinked systems. A suggested workflow for +amorphous or highly cross-linked systems. A suggested workflow for using this fix is: 1) identify a reaction to be simulated 2) build a molecule template of the reaction site before the reaction has occurred 3) build a molecule template of the reaction site after the From afaaf442d3d4d3debc5bbcbf82971f696acf3fba Mon Sep 17 00:00:00 2001 From: jrgissing Date: Thu, 6 Dec 2018 21:32:10 -0700 Subject: [PATCH 4/4] bond/react doc tweaks 2 --- doc/src/fix_bond_react.txt | 2 +- doc/utils/sphinx-config/false_positives.txt | 1 + src/USER-MISC/fix_bond_react.cpp | 4 ++-- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/doc/src/fix_bond_react.txt b/doc/src/fix_bond_react.txt index 5d6349a426..8f71cd14ec 100644 --- a/doc/src/fix_bond_react.txt +++ b/doc/src/fix_bond_react.txt @@ -318,7 +318,7 @@ 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 undeleted atom. In addition to deleting unwanted reaction +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. diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index de60206304..aab1dc12f8 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -507,6 +507,7 @@ deepskyblue defgrad deformable del +deleteIDs Dellago delocalization delocalized diff --git a/src/USER-MISC/fix_bond_react.cpp b/src/USER-MISC/fix_bond_react.cpp index 50efca8d09..600bb6a540 100644 --- a/src/USER-MISC/fix_bond_react.cpp +++ b/src/USER-MISC/fix_bond_react.cpp @@ -1653,13 +1653,13 @@ void FixBondReact::find_landlocked_atoms(int myrxn) } } - // additionally, if a deleted atom is bonded to a undeleted atom, bad + // 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 be bonded to an undeleted atom"); + error->one(FLERR,"A deleted atom cannot remain bonded to an atom that is not deleted"); } } }