Merge pull request #1697 from jrgissing/bond/react-angle-constraints

discussion of the 'update_edges' keyword. The fourth optional section

begins with the keyword 'Constraints' and lists additional criteria

that must be satisfied in order for the reaction to occur. Currently,

there is one type of constraint available, as discussed below.

there are two types of constraints available, as discussed below.

A sample map file is given below:

:line

Any number of additional constraints may be specified in the

Constraints section of the map file. Currently there is one type of

additional constraint, of type 'distance', whose syntax is as follows:

Constraints section of the map file. The constraint of type 'distance'

has syntax as follows:

distance {ID1} {ID2} {rmin} {rmax} :pre

where 'distance' is the required keyword, {ID1} and {ID2} are

pre-reaction atom IDs, and these two atoms must be separated by a

distance between {rmin} and {rmax} for the reaction to occur. This

distance between {rmin} and {rmax} for the reaction to occur.

The constraint of type 'angle' has the following syntax:

angle {ID1} {ID2} {ID3} {amin} {amax} :pre

where 'angle' is the required keyword, {ID1}, {ID2} and {ID3} are

pre-reaction atom IDs, and these three atoms must form an angle

between {amin} and {amax} for the reaction to occur (where {ID2} is

the central atom). Angles must be specified in degrees. This

constraint can be used to enforce a certain orientation between

reacting molecules.

#include "molecule.h"

#include "group.h"

#include "citeme.h"

#include "math_const.h"

#include "memory.h"

#include "error.h"

using namespace LAMMPS_NS;

using namespace FixConst;

using namespace MathConst;

static const char cite_fix_bond_react[] =

  "fix bond/react:\n\n"

#define BIG 1.0e20

#define DELTA 16

#define MAXGUESS 20 // max # of guesses allowed by superimpose algorithm

#define MAXCONARGS 5 // max # of arguments for any type of constraint

#define MAXCONARGS 7 // max # of arguments for any type of constraint + rxnID

// various statuses of superimpose algorithm:

// ACCEPT: site successfully matched to pre-reacted template

// RESTORE: restore mode, load most recent restore point

enum{ACCEPT,REJECT,PROCEED,CONTINUE,GUESSFAIL,RESTORE};

// types of available reaction constraints

enum{DISTANCE,ANGLE};

/* ---------------------------------------------------------------------- */

FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :

  global_freq = 1;

  extvector = 0;

  rxnID = 0;

  nconstraints = 0;

  status = PROCEED;

  nxspecial = NULL;

  memory->create(limit_duration,nreacts,"bond/react:limit_duration");

  memory->create(stabilize_steps_flag,nreacts,"bond/react:stabilize_steps_flag");

  memory->create(update_edges_flag,nreacts,"bond/react:update_edges_flag");

  memory->create(nconstraints,nreacts,"bond/react:nconstraints");

  memory->create(constraints,nreacts,MAXCONARGS,"bond/react:constraints");

  memory->create(constraints,1,MAXCONARGS,"bond/react:constraints");

  memory->create(iatomtype,nreacts,"bond/react:iatomtype");

  memory->create(jatomtype,nreacts,"bond/react:jatomtype");

  memory->create(ibonding,nreacts,"bond/react:ibonding");

    max_rxn[i] = INT_MAX;

    stabilize_steps_flag[i] = 0;

    update_edges_flag[i] = 0;

    nconstraints[i] = 0;

    // set default limit duration to 60 timesteps

    limit_duration[i] = 60;

    reaction_count[i] = 0;

      glove[myjbonding-1][1] = created[lcl_inst][1][rxnID];

      glove_counter++;

      // special case, only two atoms in reaction templates

      // then: bonding onemol_nxspecials guaranteed to be equal, and either 0 or 1

      if (glove_counter == onemol->natoms) {

        tagint local_atom1 = atom->map(glove[myibonding-1][1]);

        tagint local_atom2 = atom->map(glove[myjbonding-1][1]);

        if ( (nxspecial[local_atom1][0] == onemol_nxspecial[myibonding-1][0] &&

              nxspecial[local_atom2][0] == nxspecial[local_atom1][0]) &&

             (nxspecial[local_atom1][0] == 0 ||

              xspecial[local_atom1][0] == atom->tag[local_atom2]) &&

             check_constraints() ) {

          status = ACCEPT;

          glove_ghostcheck();

        } else

          status = REJECT;

      }

      avail_guesses = 0;

      for (int i = 0; i < max_natoms; i++)

int FixBondReact::check_constraints()

{

  tagint atom1,atom2;

  tagint atom1,atom2,atom3;

  double delx,dely,delz,rsq;

  double delx1,dely1,delz1,delx2,dely2,delz2;

  double rsq1,rsq2,r1,r2,c;

  double **x = atom->x;

  for (int i = 0; i < nconstraints[rxnID]; i++) {

    if (constraints[rxnID][0] == 0) { // 'distance' type

      atom1 = atom->map(glove[(int) constraints[rxnID][1]-1][1]);

      atom2 = atom->map(glove[(int) constraints[rxnID][2]-1][1]);

  for (int i = 0; i < nconstraints; i++) {

    if (constraints[i][0] == rxnID) {

      if (constraints[i][1] == DISTANCE) {

        atom1 = atom->map(glove[(int) constraints[i][2]-1][1]);

        atom2 = atom->map(glove[(int) constraints[i][3]-1][1]);

        delx = x[atom1][0] - x[atom2][0];

        dely = x[atom1][1] - x[atom2][1];

        delz = x[atom1][2] - x[atom2][2];

        domain->minimum_image(delx,dely,delz); // ghost location fix

        rsq = delx*delx + dely*dely + delz*delz;

      if (rsq < constraints[rxnID][3] || rsq > constraints[rxnID][4]) return 0;

        if (rsq < constraints[i][4] || rsq > constraints[i][5]) return 0;

      } else if (constraints[i][1] == ANGLE) {

        atom1 = atom->map(glove[(int) constraints[i][2]-1][1]);

        atom2 = atom->map(glove[(int) constraints[i][3]-1][1]);

        atom3 = atom->map(glove[(int) constraints[i][4]-1][1]);

        // 1st bond

        delx1 = x[atom1][0] - x[atom2][0];

        dely1 = x[atom1][1] - x[atom2][1];

        delz1 = x[atom1][2] - x[atom2][2];

        rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;

        r1 = sqrt(rsq1);

        // 2nd bond

        delx2 = x[atom3][0] - x[atom2][0];

        dely2 = x[atom3][1] - x[atom2][1];

        delz2 = x[atom3][2] - x[atom2][2];

        rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;

        r2 = sqrt(rsq2);

        // angle (cos and sin)

        c = delx1*delx2 + dely1*dely2 + delz1*delz2;

        c /= r1*r2;

        if (c > 1.0) c = 1.0;

        if (c < -1.0) c = -1.0;

        if (acos(c) < constraints[i][5] || acos(c) > constraints[i][6]) return 0;

      }

    }

  }

  return 1;

    if (strspn(line," \t\n\r") == strlen(line)) continue;

    if (strstr(line,"edgeIDs")) sscanf(line,"%d",&nedge);

    else if (strstr(line,"equivalences")) sscanf(line,"%d",&nequivalent);

    else if (strstr(line,"equivalences")) {

      sscanf(line,"%d",&nequivalent);

      if (nequivalent != onemol->natoms)

        error->one(FLERR,"Bond/react: Number of equivalences in map file must "

                                  "equal number of atoms in reaction templates");

    }

    else if (strstr(line,"customIDs")) sscanf(line,"%d",&ncustom);

    else if (strstr(line,"deleteIDs")) sscanf(line,"%d",&ndelete);

    else if (strstr(line,"constraints")) sscanf(line,"%d",&nconstraints[myrxn]);

    else if (strstr(line,"constraints")) sscanf(line,"%d",&nconstr);

    else break;

  }

  memory->grow(constraints,nconstraints+nconstr,MAXCONARGS,"bond/react:constraints");

  //count = NULL;

  // grab keyword and skip next line

  double tmp[MAXCONARGS];

  int n = strlen("distance") + 1;

  char *constraint_type = new char[n];

  for (int i = 0; i < nconstraints[myrxn]; i++) {

  for (int i = 0; i < nconstr; i++) {

    readline(line);

    sscanf(line,"%s",constraint_type);

    constraints[nconstraints][0] = myrxn;

    if (strcmp(constraint_type,"distance") == 0) {

      constraints[myrxn][0] = 0; // 0 = 'distance' ...maybe use another enum eventually

      constraints[nconstraints][1] = DISTANCE;

      sscanf(line,"%*s %lg %lg %lg %lg",&tmp[0],&tmp[1],&tmp[2],&tmp[3]);

      constraints[myrxn][1] = tmp[0];

      constraints[myrxn][2] = tmp[1];

      constraints[myrxn][3] = tmp[2]*tmp[2]; // using square of distance

      constraints[myrxn][4] = tmp[3]*tmp[3];

      constraints[nconstraints][2] = tmp[0];

      constraints[nconstraints][3] = tmp[1];

      constraints[nconstraints][4] = tmp[2]*tmp[2]; // using square of distance

      constraints[nconstraints][5] = tmp[3]*tmp[3];

    } else if (strcmp(constraint_type,"angle") == 0) {

      constraints[nconstraints][1] = ANGLE;

      sscanf(line,"%*s %lg %lg %lg %lg %lg",&tmp[0],&tmp[1],&tmp[2],&tmp[3],&tmp[4]);

      constraints[nconstraints][2] = tmp[0];

      constraints[nconstraints][3] = tmp[1];

      constraints[nconstraints][4] = tmp[2];

      constraints[nconstraints][5] = tmp[3]/180.0 * MY_PI;

      constraints[nconstraints][6] = tmp[4]/180.0 * MY_PI;

    } else

      error->one(FLERR,"Bond/react: Illegal constraint type in 'Constraints' section of map file");

    nconstraints++;

  }

  delete [] constraint_type;

}

  int custom_exclude_flag;

  int *stabilize_steps_flag;

  int *update_edges_flag;

  int *nconstraints;

  int nconstraints;

  double **constraints;

  int status;

  int *groupbits;

  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; // number of edge, equivalent, custom atoms in mapping file

  int nedge,nequivalent,ncustom,ndelete,nconstr; // # edge, equivalent, custom atoms in mapping file

  int attempted_rxn; // there was an attempt!

  int *local_rxn_count;

  int *ghostly_rxn_count;