Merge pull request #1795 from jrgissing/bond/react-chiral_centers

Bond/react chiral centers

doc/src/Errors_messages.rst

@@ -518,6 +518,16 @@ Doc page with :doc:`WARNING messages <Errors_warnings>`
 *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

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:

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];

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