Unverified Commit e2fd95d8 authored by Axel Kohlmeyer's avatar Axel Kohlmeyer Committed by GitHub
Browse files

Merge pull request #2240 from jrgissing/bond/react-reset_mol_ids 

support molecule ID resets in fix bond/react
parents dc241abb 10080079

doc/src/Errors_messages.rst

+2 −1
Original line number Diff line number Diff line
@@ -502,7 +502,8 @@ Doc page with :doc:`WARNING messages <Errors_warnings>` 
*Bond/react: Unknown section in map file*

   Please ensure reaction map files are properly formatted.



*Bond/react: Atom affected by reaction too close to template edge*

*Bond/react: Atom type affected by reaction too close to template edge*

*Bond/react: Bond type affected by reaction too close to template edge*

   This means an atom which changes type or connectivity during the

   reaction is too close to an 'edge' atom defined in the map

   file. This could cause incorrect assignment of bonds, angle, etc.

doc/src/fix_bond_react.rst

+21 −10
Original line number Diff line number Diff line
@@ -14,19 +14,22 @@ Syntax 
     react react-ID react-group-ID Nevery Rmin Rmax template-ID(pre-reacted) template-ID(post-reacted) map_file individual_keyword values ...

     ...



* ID, group-ID are documented in :doc:`fix <fix>` command. Group-ID is ignored.

* ID, group-ID are documented in :doc:`fix <fix>` command.

* bond/react = style name of this fix command

* the common keyword/values may be appended directly after 'bond/react'

* this applies to all reaction specifications (below)

* common_keyword = *stabilization*

* common_keyword = *stabilization* or *reset_mol_ids*



  .. parsed-literal::



       *stabilization* values = *no* or *yes* *group-ID* *xmax*

         *no* = no reaction site stabilization

         *no* = no reaction site stabilization (default)

         *yes* = perform reaction site stabilization

           *group-ID* = user-assigned prefix for the dynamic group of atoms not currently involved in a reaction

           *xmax* = xmax value that is used by an internally-created :doc:`nve/limit <fix_nve_limit>` integrator

       *reset_mol_ids* values = *yes* or *no*

         *yes* = update molecule IDs based on new global topology (default)

         *no* = do not update molecule IDs



* react = mandatory argument indicating new reaction specification

* react-ID = user-assigned name for the reaction
@@ -50,9 +53,9 @@ Syntax 
         *stabilize_steps* value = timesteps

           timesteps = number of timesteps to apply the internally-created :doc:`nve/limit <fix_nve_limit>` fix to reacting atoms

         *update_edges* value = *none* or *charges* or *custom*

           none = do not update topology near the edges of reaction templates

           charges = update atomic charges of all atoms in reaction templates

           custom = force the update of user-specified atomic charges

           *none* = do not update topology near the edges of reaction templates

           *charges* = update atomic charges of all atoms in reaction templates

           *custom* = force the update of user-specified atomic charges



Examples

""""""""
@@ -154,6 +157,13 @@ due to the internal dynamic grouping performed by fix bond/react. 
   If the group-ID is an existing static group, react-group-IDs

   should also be specified as this static group, or a subset.



The *reset_mol_ids* keyword invokes the :doc:`reset_mol_ids <reset_mol_ids>`

command after a reaction occurs, to ensure that molecule IDs are

consistent with the new bond topology. The group-ID used for

:doc:`reset_mol_ids <reset_mol_ids>` is the group-ID for this fix.

Resetting molecule IDs is necessarily a global operation, and so can

be slow for very large systems.



The following comments pertain to each *react* argument (in other

words, can be customized for each reaction, or reaction step):
@@ -203,9 +213,10 @@ surrounding topology. As described below, the bonding atom pairs of 
the pre-reacted template are specified by atom ID in the map file. The

pre-reacted molecule template should contain as few atoms as possible

while still completely describing the topology of all atoms affected

by the reaction. For example, if the force field contains dihedrals,

the pre-reacted template should contain any atom within three bonds of

reacting atoms.

by the reaction (which includes all atoms that change atom type or

connectivity, and all bonds that change bond type). For example, if

the force field contains dihedrals, the pre-reacted template should

contain any atom within three bonds of reacting atoms.



Some atoms in the pre-reacted template that are not reacting may have

missing topology with respect to the simulation. For example, the
@@ -554,7 +565,7 @@ Default 
"""""""



The option defaults are stabilization = no, prob = 1.0, stabilize_steps = 60,

update_edges = none

reset_mol_ids = yes, update_edges = none



----------

src/USER-REACTION/fix_bond_react.cpp

+32 −12
Original line number Diff line number Diff line
@@ -33,6 +33,7 @@ Contributing Author: Jacob Gissinger (jacob.gissinger@colorado.edu) 
#include "neigh_list.h"

#include "neigh_request.h"

#include "random_mars.h"

#include "reset_mol_ids.h"

#include "molecule.h"

#include "group.h"

#include "citeme.h"
@@ -92,6 +93,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : 
  fix1 = NULL;

  fix2 = NULL;

  fix3 = NULL;

  reset_mol_ids = NULL;



  if (narg < 8) error->all(FLERR,"Illegal fix bond/react command: "

                           "too few arguments");
@@ -144,7 +146,8 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : 


  int iarg = 3;

  stabilization_flag = 0;

  int num_common_keywords = 1;

  reset_mol_ids_flag = 1;

  int num_common_keywords = 2;

  for (int m = 0; m < num_common_keywords; m++) {

    if (strcmp(arg[iarg],"stabilization") == 0) {

      if (strcmp(arg[iarg+1],"no") == 0) {
@@ -162,11 +165,23 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : 
        nve_limit_xmax = arg[iarg+3];

        iarg += 4;

      }

    } else if (strcmp(arg[iarg],"reset_mol_ids") == 0) {

      if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/react command: "

                                    "'reset_mol_ids' keyword has too few arguments");

      if (strcmp(arg[iarg+1],"yes") == 0) ; // default

      if (strcmp(arg[iarg+1],"no") == 0) reset_mol_ids_flag = 0;

      iarg += 2;

    } else if (strcmp(arg[iarg],"react") == 0) {

      break;

    } else error->all(FLERR,"Illegal fix bond/react command: unknown keyword");

  }



  if (reset_mol_ids_flag) {

    delete reset_mol_ids;

    reset_mol_ids = new ResetMolIDs(lmp);

    reset_mol_ids->create_computes(id,group->names[igroup]);

  }



  // set up common variables as vectors of length 'nreacts'

  // nevery, cutoff, onemol, twomol, superimpose file
@@ -229,11 +244,11 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : 


    int n = strlen(arg[iarg]) + 1;

    if (n > MAXLINE) error->all(FLERR,"Reaction name (react-ID) is too long (limit: 256 characters)");

    strncpy(rxn_name[rxn],arg[iarg++],n);

    strcpy(rxn_name[rxn],arg[iarg++]);



    int igroup = group->find(arg[iarg++]);

    if (igroup == -1) error->all(FLERR,"Could not find fix group ID");

    groupbits[rxn] = group->bitmask[igroup];

    int groupid = group->find(arg[iarg++]);

    if (groupid == -1) error->all(FLERR,"Could not find fix group ID");

    groupbits[rxn] = group->bitmask[groupid];



    if (strncmp(arg[iarg],"v_",2) == 0) {

      n = strlen(&arg[iarg][2]) + 1;
@@ -499,6 +514,8 @@ FixBondReact::~FixBondReact() 
  }

  delete [] random;



  delete reset_mol_ids;



  memory->destroy(partner);

  memory->destroy(finalpartner);

  memory->destroy(ncreate);
@@ -623,9 +640,9 @@ void FixBondReact::post_constructor() 
  group->assign(cmd);



  if (stabilization_flag == 1) {

    int igroup = group->find(exclude_group);

    int groupid = group->find(exclude_group);

    // create exclude_group if not already existing, or use as parent group if static

    if (igroup == -1 || group->dynamic[igroup] == 0) {

    if (groupid == -1 || group->dynamic[groupid] == 0) {

      // create stabilization per-atom property

      cmd = std::string("bond_react_stabilization_internal");

      id_fix3 = new char[cmd.size()+1];
@@ -655,7 +672,7 @@ void FixBondReact::post_constructor() 
      strcat(exclude_group,"_REACT");



      group->find_or_create(exclude_group);

      if (igroup == -1)

      if (groupid == -1)

        cmd = fmt::format("{} dynamic all property statted_tags",exclude_group);

      else

        cmd = fmt::format("{} dynamic {} property statted_tags",exclude_group,exclude_PARENT_group);
@@ -2061,7 +2078,7 @@ void FixBondReact::find_landlocked_atoms(int myrxn) 
  for (int i = 0; i < twomol->natoms; i++) {

    if (twomol->type[i] != onemol->type[equivalences[i][1][myrxn]-1] && landlocked_atoms[i][myrxn] == 0) {

      char str[128];

      snprintf(str,128,"Bond/react: Atom affected by reaction %s too close to template edge",rxn_name[myrxn]);

      snprintf(str,128,"Bond/react: Atom type affected by reaction %s too close to template edge",rxn_name[myrxn]);

      error->all(FLERR,str);

    }

  }
@@ -2080,7 +2097,7 @@ void FixBondReact::find_landlocked_atoms(int myrxn) 
            if (onemol_batom == equivalences[twomol_atomj-1][1][myrxn]) {

              if (twomol->bond_type[i][j] != onemol->bond_type[onemol_atomi-1][m]) {

                char str[128];

                snprintf(str,128,"Bond/react: Atom affected by reaction %s too close to template edge",rxn_name[myrxn]);

                snprintf(str,128,"Bond/react: Bond type affected by reaction %s too close to template edge",rxn_name[myrxn]);

                error->all(FLERR,str);

              }

            }
@@ -2092,7 +2109,7 @@ void FixBondReact::find_landlocked_atoms(int myrxn) 
              if (onemol_batom == equivalences[i][1][myrxn]) {

                if (twomol->bond_type[i][j] != onemol->bond_type[onemol_atomj-1][m]) {

                  char str[128];

                  snprintf(str,128,"Bond/react: Atom affected by reaction %s too close to template edge",rxn_name[myrxn]);

                  snprintf(str,128,"Bond/react: Bond type affected by reaction %s too close to template edge",rxn_name[myrxn]);

                  error->all(FLERR,str);

                }

              }
@@ -2503,7 +2520,7 @@ void FixBondReact::ghost_glovecast() 
}



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

update charges, types, special lists and all topology

update molecule IDs, charges, types, special lists and all topology

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



void FixBondReact::update_everything()
@@ -3042,6 +3059,9 @@ void FixBondReact::update_everything() 
  atom->natoms -= ndel;

  // done deleting atoms



  // reset mol ids

  if (reset_mol_ids_flag) reset_mol_ids->reset();



  // 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

src/USER-REACTION/fix_bond_react.h

+5 −1
Original line number Diff line number Diff line
@@ -61,6 +61,7 @@ class FixBondReact : public Fix { 
  int *max_rxn,*nlocalskips,*nghostlyskips;

  tagint lastcheck;

  int stabilization_flag;

  int reset_mol_ids_flag;

  int custom_exclude_flag;

  int *stabilize_steps_flag;

  int *update_edges_flag;
@@ -93,6 +94,7 @@ class FixBondReact : public Fix { 
  class RanMars **random; // random number for 'prob' keyword

  class RanMars **rrhandom; // random number for Arrhenius constraint

  class NeighList *list;

  class ResetMolIDs *reset_mol_ids; // class for resetting mol IDs



  int *reacted_mol,*unreacted_mol;

  int *limit_duration; // indicates how long to relax
@@ -240,8 +242,10 @@ E: Bond/react: Invalid template atom ID in map file 
Atom IDs in molecule templates range from 1 to the number of atoms in the template.



E or W: Bond/react: Atom affected by reaction %s too close to template edge

        Bond/react: Atom type affected by reaction %s too close to template edge

        Bond/react: Bond type affected by reaction %s too close to template edge



This means an atom which changes type or connectivity during the

This means an atom (or bond) that changes type or connectivity during the

reaction is too close to an 'edge' atom defined in the map file. This

could cause incorrect assignment of bonds, angle, etc. Generally, this

means you must include more atoms in your templates, such that there

src/reset_mol_ids.cpp

+85 −46
Original line number Diff line number Diff line
@@ -17,7 +17,6 @@ 


#include "reset_mol_ids.h"

#include <mpi.h>

#include <string>

#include "atom.h"

#include "domain.h"

#include "comm.h"
@@ -33,10 +32,29 @@ using namespace LAMMPS_NS; 


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



ResetMolIDs::ResetMolIDs(LAMMPS *lmp) : Pointers(lmp) {}

ResetMolIDs::ResetMolIDs(LAMMPS *lmp) : Pointers(lmp) {

  cfa = NULL;

  cca = NULL;



  // default settings



  compressflag = 1;

  singleflag = 0;

  offset = -1;

}



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



ResetMolIDs::~ResetMolIDs()

{

  modify->delete_compute(idfrag);

  if (compressflag) modify->delete_compute(idchunk);

}



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

   called as reset_mol_ids command in input script

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



void ResetMolIDs::command(int narg, char **arg)

{

  if (domain->box_exist == 0)
@@ -49,13 +67,7 @@ void ResetMolIDs::command(int narg, char **arg) 
  // process args



  if (narg < 1) error->all(FLERR,"Illegal reset_mol_ids command");

  int igroup = group->find(arg[0]);

  if (igroup == -1) error->all(FLERR,"Could not find reset_mol_ids group ID");

  int groupbit = group->bitmask[igroup];



  int compressflag = 1;

  int singleflag = 0;

  tagint offset = -1;

  char *groupid = arg[0];



  int iarg = 1;

  while (iarg < narg) {
@@ -86,20 +98,6 @@ void ResetMolIDs::command(int narg, char **arg) 
  MPI_Barrier(world);

  double time1 = MPI_Wtime();



  // create instances of compute fragment/atom, compute reduce (if needed),

  // and compute chunk/atom.  all use the group-ID for this command



  const std::string idfrag = "reset_mol_ids_FRAGMENT_ATOM";

  if (singleflag)

    modify->add_compute(fmt::format("{} {} fragment/atom single yes",idfrag,arg[0]));

  else

    modify->add_compute(fmt::format("{} {} fragment/atom single no",idfrag,arg[0]));



  const std::string idchunk = "reset_mol_ids_CHUNK_ATOM";

  if (compressflag)

    modify->add_compute(fmt::format("{} {} chunk/atom molecule compress yes",

                                    idchunk,arg[0]));



  // initialize system since comm->borders() will be invoked



  lmp->init();
@@ -116,12 +114,73 @@ void ResetMolIDs::command(int narg, char **arg) 
  comm->borders();

  if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost);



  // create computes



  create_computes((char *) "COMMAND",groupid);



  // reset molecule IDs



  reset();



  // total time



  MPI_Barrier(world);



  if (comm->me == 0) {

    if (nchunk < 0)

      utils::logmesg(lmp,fmt::format("  number of new molecule IDs = unknown\n"));

    else

      utils::logmesg(lmp,fmt::format("  number of new molecule IDs = {}\n",nchunk));

    utils::logmesg(lmp,fmt::format("  reset_mol_ids CPU = {:.3f} seconds\n",

                                   MPI_Wtime()-time1));

  }

}



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

   create computes used by reset_mol_ids

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



void ResetMolIDs::create_computes(char *fixid, char *groupid)

{

  int igroup = group->find(groupid);

  if (igroup == -1) error->all(FLERR,"Could not find reset_mol_ids group ID");

  groupbit = group->bitmask[igroup];



  // create instances of compute fragment/atom, compute reduce (if needed),

  // and compute chunk/atom.  all use the group-ID for this command.

  // 'fixid' allows for creating independent instances of the computes



  idfrag = fmt::format("{}_reset_mol_ids_FRAGMENT_ATOM",fixid);

  if (singleflag)

    modify->add_compute(fmt::format("{} {} fragment/atom single yes",idfrag,groupid));

  else

    modify->add_compute(fmt::format("{} {} fragment/atom single no",idfrag,groupid));



  idchunk = fmt::format("{}_reset_mol_ids_CHUNK_ATOM",fixid);

  if (compressflag)

    modify->add_compute(fmt::format("{} {} chunk/atom molecule compress yes",

                                    idchunk,groupid));



  int icompute = modify->find_compute(idfrag);

  cfa = (ComputeFragmentAtom *) modify->compute[icompute];





  if (compressflag) {

    icompute = modify->find_compute(idchunk);

    cca = (ComputeChunkAtom *) modify->compute[icompute];

  }

}



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

   called from command() and directly from fixes that update molecules

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



void ResetMolIDs::reset()

{

  // invoke peratom method of compute fragment/atom

  // walks bond connectivity and assigns each atom a fragment ID

  // if singleflag = 0, atoms w/out bonds will be assigned fragID = 0



  int icompute = modify->find_compute(idfrag);

  ComputeFragmentAtom *cfa = (ComputeFragmentAtom *) modify->compute[icompute];

  cfa->compute_peratom();

  double *fragIDs = cfa->vector_atom;
@@ -138,15 +197,13 @@ void ResetMolIDs::command(int narg, char **arg) 
  // if compressflag = 0, done

  // set nchunk = -1 since cannot easily determine # of new molecule IDs



  int nchunk = -1;

  nchunk = -1;



  // if compressflag = 1, invoke peratom method of compute chunk/atom

  // will compress new molecule IDs to be contiguous 1 to Nmol

  // NOTE: use of compute chunk/atom limits Nmol to a 32-bit int



  if (compressflag) {

    icompute = modify->find_compute(idchunk);

    ComputeChunkAtom *cca = (ComputeChunkAtom *) modify->compute[icompute];

    cca->compute_peratom();

    double *chunkIDs = cca->vector_atom;

    nchunk = cca->nchunk;
@@ -193,22 +250,4 @@ void ResetMolIDs::command(int narg, char **arg) 
      }

    }

  }



  // clean up



  modify->delete_compute(idfrag);

  if (compressflag) modify->delete_compute(idchunk);



  // total time



  MPI_Barrier(world);



  if (comm->me == 0) {

    if (nchunk < 0)

      utils::logmesg(lmp,fmt::format("  number of new molecule IDs = unknown\n"));

    else

      utils::logmesg(lmp,fmt::format("  number of new molecule IDs = {}\n",nchunk));

    utils::logmesg(lmp,fmt::format("  reset_mol_ids CPU = {:.3f} seconds\n",

                                   MPI_Wtime()-time1));

  }

}

CRA Git | Maintained and supported by SUSTech CRA and CCSE