diff --git a/doc/src/Errors_messages.rst b/doc/src/Errors_messages.rst index 71de83afa9..a1795a6776 100644 --- a/doc/src/Errors_messages.rst +++ b/doc/src/Errors_messages.rst @@ -518,6 +518,16 @@ Doc page with :doc:`WARNING messages ` *Bond/react: A deleted atom cannot remain bonded to an atom that is not deleted* Self-explanatory. +*Bond/react: First neighbors of chiral atoms must be of mutually different types* + Self-explanatory. + +*Bond/react: Chiral atoms must have exactly four first neighbors* + Self-explanatory. + +*Bond/react: Molecule template 'Coords' section required for chiralIDs keyword* + The coordinates of atoms in the pre-reacted template are used to determine + chirality. + *Bond/react special bond generation overflow* The number of special bonds per-atom created by a reaction exceeds the system setting. See the read\_data or create\_box command for how to diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst index ab63aa232c..b2866eb9c7 100644 --- a/doc/src/fix_bond_react.rst +++ b/doc/src/fix_bond_react.rst @@ -20,9 +20,9 @@ Syntax * the common keyword/values may be appended directly after 'bond/react' * this applies to all reaction specifications (below) * common\_keyword = *stabilization* - + .. parsed-literal:: - + *stabilization* values = *no* or *yes* *group-ID* *xmax* *no* = no reaction site stabilization *yes* = perform reaction site stabilization @@ -40,9 +40,9 @@ Syntax * map\_file = name of file specifying corresponding atom-IDs in the pre- and post-reacted templates * zero or more individual keyword/value pairs may be appended to each react argument * individual\_keyword = *prob* or *max\_rxn* or *stabilize\_steps* or *update\_edges* - + .. parsed-literal:: - + *prob* values = fraction seed fraction = initiate reaction with this probability if otherwise eligible seed = random number seed (positive integer) @@ -253,7 +253,7 @@ 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 four optional keywords. The +contains one mandatory keyword and five optional keywords. The mandatory keyword is 'equivalences': @@ -269,10 +269,11 @@ The optional keywords are 'edgeIDs', 'deleteIDs', 'customIDs' and N *edgeIDs* = # of edge atoms N in the pre-reacted molecule template N *deleteIDs* = # of atoms N that are specified for deletion + N *chiralIDs* = # of specified chiral centers N N *customIDs* = # of atoms N that are specified for a custom update N *constraints* = # of specified reaction constraints N -The body of the map file contains two mandatory sections and four +The body of the map file contains two mandatory sections and five 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 @@ -284,12 +285,14 @@ 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 'DeleteIDs' and lists the atom IDs of pre-reaction template atoms to -delete. The third optional section begins with the keyword 'Custom +delete. The third optional section begins with the keyword 'ChiralIDs' +lists the atom IDs of chiral atoms whose handedness should be +enforced. The fourth 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. The fourth optional section +discussion of the 'update\_edges' keyword. The fifth optional section begins with the keyword 'Constraints' and lists additional criteria that must be satisfied in order for the reaction to occur. Currently, there are three types of constraints available, as discussed below. @@ -332,6 +335,15 @@ A sample map file is given below: ---------- +The handedness of atoms that are chiral centers can be enforced by +listing their IDs in the ChiralIDs section. A chiral atom must be +bonded to four atoms with mutually different atom types. This feature +uses the coordinates and types of the involved atoms in the +pre-reaction template to determine handedness. Three atoms bonded to +the chiral center are arbitrarily chosen, to define an oriented plane, +and the relative position of the fourth bonded atom determines the +chiral center's handedness. + Any number of additional constraints may be specified in the Constraints section of the map file. The constraint of type 'distance' has syntax as follows: diff --git a/src/USER-MISC/fix_bond_react.cpp b/src/USER-MISC/fix_bond_react.cpp index 33aab5b4ad..2995a596a5 100644 --- a/src/USER-MISC/fix_bond_react.cpp +++ b/src/USER-MISC/fix_bond_react.cpp @@ -36,6 +36,7 @@ Contributing Author: Jacob Gissinger (jacob.gissinger@colorado.edu) #include "group.h" #include "citeme.h" #include "math_const.h" +#include "math_extra.h" #include "memory.h" #include "error.h" @@ -297,6 +298,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : 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"); + memory->create(chiral_atoms,max_natoms,6,nreacts,"bond/react:chiral_atoms"); for (int j = 0; j < nreacts; j++) for (int i = 0; i < max_natoms; i++) { @@ -304,6 +306,9 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : if (update_edges_flag[j] == 1) custom_edges[i][j] = 1; else custom_edges[i][j] = 0; delete_atoms[i][j] = 0; + for (int k = 0; k < 6; k++) { + chiral_atoms[i][k][j] = 0; + } } // read all map files afterward @@ -430,6 +435,7 @@ FixBondReact::~FixBondReact() memory->destroy(landlocked_atoms); memory->destroy(custom_edges); memory->destroy(delete_atoms); + memory->destroy(chiral_atoms); memory->destroy(nevery); memory->destroy(cutsq); @@ -1675,6 +1681,7 @@ int FixBondReact::check_constraints() delx1 = x[atom1][0] - x[atom2][0]; dely1 = x[atom1][1] - x[atom2][1]; delz1 = x[atom1][2] - x[atom2][2]; + domain->minimum_image(delx1,dely1,delz1); // ghost location fix rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; r1 = sqrt(rsq1); @@ -1682,6 +1689,7 @@ int FixBondReact::check_constraints() delx2 = x[atom3][0] - x[atom2][0]; dely2 = x[atom3][1] - x[atom2][1]; delz2 = x[atom3][2] - x[atom2][2]; + domain->minimum_image(delx2,dely2,delz2); // ghost location fix rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; r2 = sqrt(rsq2); @@ -1699,6 +1707,30 @@ int FixBondReact::check_constraints() } } } + + // let's also check chirality within 'check_constraint' + for (int i = 0; i < onemol->natoms; i++) { + if (chiral_atoms[i][0][rxnID] == 1) { + double my4coords[12]; + // already ensured, by transitive property, that chiral simulation atom has four neighs + for (int j = 0; j < 4; j++) { + atom1 = atom->map(glove[i][1]); + // loop over known types involved in chiral center + for (int jj = 0; jj < 4; jj++) { + if (atom->type[atom->map(xspecial[atom1][j])] == chiral_atoms[i][jj+2][rxnID]) { + atom2 = atom->map(xspecial[atom1][j]); + atom2 = domain->closest_image(atom1,atom2); + for (int k = 0; k < 3; k++) { + my4coords[3*jj+k] = x[atom2][k]; + } + break; + } + } + } + if (get_chirality(my4coords) != chiral_atoms[i][1][rxnID]) return 0; + } + } + return 1; } @@ -1739,6 +1771,33 @@ double FixBondReact::get_temperature() return t; } +/* ---------------------------------------------------------------------- +return handedness (1 or -1) of a chiral center, given ordered set of coordinates +------------------------------------------------------------------------- */ + +int FixBondReact::get_chirality(double four_coords[12]) +{ + // define oriented plane with first three coordinates + double vec1[3],vec2[3],vec3[3],vec4[3],mean3[3],dot; + + for (int i = 0; i < 3; i++) { + vec1[i] = four_coords[i]-four_coords[i+3]; + vec2[i] = four_coords[i+3]-four_coords[i+6]; + } + + MathExtra::cross3(vec1,vec2,vec3); + + for (int i = 0; i < 3; i++) { + mean3[i] = (four_coords[i] + four_coords[i+3] + + four_coords[i+6])/3; + vec4[i] = four_coords[i+9] - mean3[i]; + } + + dot = MathExtra::dot3(vec3,vec4); + dot = dot/fabs(dot); + return (int) dot; +} + /* ---------------------------------------------------------------------- Get xspecials for current molecule templates ------------------------------------------------------------------------- */ @@ -2867,6 +2926,7 @@ void FixBondReact::read(int myrxn) } else if (strstr(line,"customIDs")) sscanf(line,"%d",&ncustom); else if (strstr(line,"deleteIDs")) sscanf(line,"%d",&ndelete); + else if (strstr(line,"chiralIDs")) sscanf(line,"%d",&nchiral); else if (strstr(line,"constraints")) { sscanf(line,"%d",&nconstr); memory->grow(constraints,nconstraints+nconstr,MAXCONARGS,"bond/react:constraints"); @@ -2898,6 +2958,8 @@ void FixBondReact::read(int myrxn) CustomEdges(line, myrxn); } else if (strcmp(keyword,"DeleteIDs") == 0) { DeleteAtoms(line, myrxn); + } else if (strcmp(keyword,"ChiralIDs") == 0) { + ChiralCenters(line, myrxn); } else if (strcmp(keyword,"Constraints") == 0) { Constraints(line, myrxn); } else error->one(FLERR,"Bond/react: Unknown section in map file"); @@ -2975,6 +3037,37 @@ void FixBondReact::DeleteAtoms(char *line, int myrxn) } } +void FixBondReact::ChiralCenters(char *line, int myrxn) +{ + int tmp; + for (int i = 0; i < nchiral; i++) { + readline(line); + sscanf(line,"%d",&tmp); + chiral_atoms[tmp-1][0][myrxn] = 1; + if (onemol->xflag == 0) + error->one(FLERR,"Bond/react: Molecule template 'Coords' section required for chiralIDs keyword"); + if ((int) onemol_nxspecial[tmp-1][0] != 4) + error->one(FLERR,"Bond/react: Chiral atoms must have exactly four first neighbors"); + for (int j = 0; j < 4; j++) { + for (int k = j+1; k < 4; k++) { + if (onemol->type[onemol_xspecial[tmp-1][j]-1] == + onemol->type[onemol_xspecial[tmp-1][k]-1]) + error->one(FLERR,"Bond/react: First neighbors of chiral atoms must be of mutually different types"); + } + } + // record order of atom types, and coords + double my4coords[12]; + for (int j = 0; j < 4; j++) { + chiral_atoms[tmp-1][j+2][myrxn] = onemol->type[onemol_xspecial[tmp-1][j]-1]; + for (int k = 0; k < 3; k++) { + my4coords[3*j+k] = onemol->x[onemol_xspecial[tmp-1][j]-1][k]; + } + } + // get orientation + chiral_atoms[tmp-1][1][myrxn] = get_chirality(my4coords); + } +} + void FixBondReact::Constraints(char *line, int myrxn) { double tmp[MAXCONARGS]; diff --git a/src/USER-MISC/fix_bond_react.h b/src/USER-MISC/fix_bond_react.h index f59e051ed1..a5bf2bc587 100644 --- a/src/USER-MISC/fix_bond_react.h +++ b/src/USER-MISC/fix_bond_react.h @@ -110,7 +110,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,ndelete,nconstr; // # edge, equivalent, custom atoms in mapping file + int nedge,nequivalent,ncustom,ndelete,nchiral,nconstr; // # edge, equivalent, custom atoms in mapping file int attempted_rxn; // there was an attempt! int *local_rxn_count; int *ghostly_rxn_count; @@ -126,6 +126,7 @@ class FixBondReact : public Fix { 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 ***chiral_atoms; // pre-react chiral atoms. 1) flag 2) orientation 3-4) ordered atom types int **nxspecial,**onemol_nxspecial,**twomol_nxspecial; // full number of 1-4 neighbors tagint **xspecial,**onemol_xspecial,**twomol_xspecial; // full 1-4 neighbor list @@ -149,6 +150,7 @@ class FixBondReact : public Fix { void Equivalences(char *,int); void CustomEdges(char *,int); void DeleteAtoms(char *,int); + void ChiralCenters(char *,int); void Constraints(char *,int); void make_a_guess (); @@ -159,6 +161,7 @@ class FixBondReact : public Fix { void ring_check(); int check_constraints(); double get_temperature(); + int get_chirality(double[12]); // get handedness given an ordered set of coordinates void open(char *); void readline(char *); @@ -249,6 +252,18 @@ E: Bond/react: A deleted atom cannot remain bonded to an atom that is not delete Self-explanatory. +E: Bond/react: First neighbors of chiral atoms must be of mutually different types + +Self-explanatory. + +E: Bond/react: Chiral atoms must have exactly four first neighbors + +Self-explanatory. + +E: Bond/react: Molecule template 'Coords' section required for chiralIDs keyword + +The coordinates of atoms in the pre-reacted template are used to determine chirality. + E: Bond/react special bond generation overflow The number of special bonds per-atom created by a reaction exceeds the