Commit 5821a5ec authored by jrgissing's avatar jrgissing
Browse files

bond/react: custom group with stabilization option

parent 20b9c7fd
Loading
Loading
Loading
Loading
+21 −13
Original line number Diff line number Diff line
@@ -102,19 +102,27 @@ involved in any new reactions. The {xmax} value keyword should
typically be set to the maximum distance that non-reacting atoms move
during the simulation.

The group-ID set using the {stabilization} keyword should be a
previously unused group-ID. It cannot be specified as 'all'. The fix
bond/react command creates a "dynamic group"_group.html of this name
that includes all non-reacting atoms. This dynamic group-ID should
then be used by a subsequent system-wide time integrator such as nvt,
npt, or nve, as shown in the second example above. It is currently
necessary to place the time integration command after the fix
bond/react command due to the internal dynamic grouping performed by
fix bond/react.

NOTE: The internally created group currently applies to all atoms in
the system, i.e. you should generally not have a separate thermostat
which acts on the 'all' group.
The group-ID set using the {stabilization} keyword can be an existing
static group or a previously-unused group-ID. It cannot be specified
as 'all'. If the group-ID is previously unused, fix bond/react command
creates a "dynamic group"_group.html of this name that is initialized
to include all atoms. If the group-ID is that of an existing static
group, the group is converted into a dynamic group, whose atoms are
limited to those belonging to the original static group. In either
case, this dynamic group-ID should then be used by a subsequent
system-wide time integrator such as nvt, npt, or nve, as shown in the
second example above. The time integration command should be placed
after the fix bond/react command due to the internal dynamic grouping
performed by fix bond/react. By specifying an existing group, you may
thermostat non-reacting parts of your system separately.

NOTE: If the group-ID is an existing static group, react-group-IDs
should also be specified as this group, or a subset.

NOTE: If the group-ID is previously unused, the internally created
group applies to all atoms in the system, i.e. you should generally
not have a separate thermostat which acts on the 'all' group, or any
other group.

The following comments pertain to each {react} argument:

+100 −59
Original line number Diff line number Diff line
@@ -354,6 +354,9 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :

  id_fix1 = NULL;
  id_fix2 = NULL;
  id_fix3 = NULL;
  statted_id = NULL;
  custom_exclude_flag = 0;
}

/* ---------------------------------------------------------------------- */
@@ -426,11 +429,17 @@ FixBondReact::~FixBondReact()
    // check nfix in case all fixes have already been deleted
    if (id_fix1 == NULL && modify->nfix) modify->delete_fix(id_fix1);
    delete [] id_fix1;

    if (custom_exclude_flag == 0) {
      if (id_fix3 == NULL && modify->nfix) modify->delete_fix(id_fix3);
      delete [] id_fix3;
    }
  }

  if (id_fix2 == NULL && modify->nfix) modify->delete_fix(id_fix2);
  delete [] id_fix2;

  delete [] statted_id;
  delete [] guess_branch;
  delete [] pioneer_count;
}
@@ -461,21 +470,20 @@ void FixBondReact::post_constructor()

  int ifix = modify->find_fix(id_fix2);
  if (ifix == -1) {
    char **newarg = new char*[8];
    char **newarg = new char*[7];
    newarg[0] = (char *) "bond_react_props_internal";
    newarg[1] = (char *) "all"; // group ID is ignored
    newarg[2] = (char *) "property/atom";
    newarg[3] = (char *) "i_limit_tags";
    newarg[4] = (char *) "i_statted_tags";
    newarg[5] = (char *) "i_react_tags";
    newarg[6] = (char *) "ghost";
    newarg[7] = (char *) "yes";
    modify->add_fix(8,newarg);
    fix2 = modify->fix[modify->nfix-1];
    newarg[4] = (char *) "i_react_tags";
    newarg[5] = (char *) "ghost";
    newarg[6] = (char *) "yes";
    modify->add_fix(7,newarg);
    delete [] newarg;
  }

  // create master_group if not already existing
  // NOTE: limit_tags and react_tags automaticaly intitialized to zero (unless read from restart)
  if (group->find(master_group) == -1) {
    group->find_or_create(master_group);
    char **newarg;
@@ -489,33 +497,100 @@ void FixBondReact::post_constructor()
    delete [] newarg;
  }

  // on to statted_tags (system-wide thermostat)
  // intialize per-atom statted_flags to 1
  // (only if not already initialized by restart)
  // NOTE: limit_tags and react_tags automaticaly intitialized to zero (unless read from restart)
  if (fix2->restart_reset != 1) {
    int flag;
    int index = atom->find_custom("statted_tags",flag);
    int *i_statted_tags = atom->ivector[index];
  if (stabilization_flag == 1) {
    int igroup = 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) {
      // create stabilization per-atom property
      len = strlen("bond_react_stabilization_internal") + 1;
      id_fix3 = new char[len];
      strcpy(id_fix3,"bond_react_stabilization_internal");

      ifix = modify->find_fix(id_fix3);
      if (ifix == -1) {
        char **newarg = new char*[6];
        newarg[0] = (char *) "bond_react_stabilization_internal";
        newarg[1] = (char *) "all"; // group ID is ignored
        newarg[2] = (char *) "property/atom";
        newarg[3] = (char *) "i_statted_tags";
        newarg[4] = (char *) "ghost";
        newarg[5] = (char *) "yes";
        modify->add_fix(6,newarg);
        fix2 = modify->fix[modify->nfix-1];
        delete [] newarg;
      }

    for (int i = 0; i < atom->nlocal; i++)
      i_statted_tags[i] = 1;
      len = strlen("statted_tags") + 1;
      statted_id = new char[len];
      strcpy(statted_id,"statted_tags");

      // if static group exists, duplicate it, use duplicate as parent group
      // original will be converted into dynamic per-atom property group
      if (igroup != -1) {
        char **newarg;
        newarg = new char*[3];
        newarg[0] = (char *) "exclude_PARENT_group";
        newarg[1] = (char *) "union";
        newarg[2] = exclude_group;
        group->assign(3,newarg);
        delete [] newarg;
      }

  if (stabilization_flag == 1) {
    // create exclude_group if not already existing
    if (group->find(exclude_group) == -1) {
      group->find_or_create(exclude_group);
      char **newarg;
      newarg = new char*[5];
      newarg[0] = exclude_group;
      newarg[1] = (char *) "dynamic";
      newarg[2] = (char *) "all";
      if (igroup == -1) newarg[2] = (char *) "all";
      else newarg[2] = (char *) "exclude_PARENT_group";
      newarg[3] = (char *) "property";
      newarg[4] = (char *) "statted_tags";
      group->assign(5,newarg);
      delete [] newarg;

      // on to statted_tags (system-wide thermostat)
      // intialize per-atom statted_flags to 1
      // (only if not already initialized by restart)
      if (fix2->restart_reset != 1) {
        int flag;
        int index = atom->find_custom("statted_tags",flag);
        int *i_statted_tags = atom->ivector[index];

        for (int i = 0; i < atom->nlocal; i++)
          i_statted_tags[i] = 1;
      }
    } else {
        // sleeping code, for future capabilities
        custom_exclude_flag = 1;
        // first we have to find correct fix group reference
        int n = strlen("GROUP_") + strlen(exclude_group) + 1;
        char *fix_group = new char[n];
        strcpy(fix_group,"GROUP_");
        strcat(fix_group,exclude_group);
        int ifix = modify->find_fix(fix_group);
        Fix *fix = modify->fix[ifix];
        delete [] fix_group;

        // this returns names of corresponding property
        int unused;
        char * idprop;
        idprop = (char *) fix->extract("property",unused);
        if (idprop == NULL)
          error->all(FLERR,"Exclude group must be a per-atom property group");

        len = strlen(idprop) + 1;
        statted_id = new char[len];
        strcpy(statted_id,idprop);

        // intialize per-atom statted_tags to 1
        // need to correct for smooth restarts
        //int flag;
        //int index = atom->find_custom(statted_id,flag);
        //int *i_statted_tags = atom->ivector[index];
        //for (int i = 0; i < atom->nlocal; i++)
        //  i_statted_tags[i] = 1;
      }


    // let's create a new nve/limit fix to limit newly reacted atoms
    len = strlen("bond_react_MASTER_nve_limit") + 1;
@@ -534,40 +609,6 @@ void FixBondReact::post_constructor()
      fix1 = modify->fix[modify->nfix-1];
      delete [] newarg;
    }

  }

  // currently must redefine dynamic groups so they are updated at proper time
  // -> should double check as to why

  int must_redefine_groups = 1;

  if (must_redefine_groups) {
    group->find_or_create(master_group);
    char **newarg;
    newarg = new char*[5];
    newarg[0] = master_group;
    newarg[1] = (char *) "dynamic";
    newarg[2] = (char *) "all";
    newarg[3] = (char *) "property";
    newarg[4] = (char *) "limit_tags";
    group->assign(5,newarg);
    delete [] newarg;
  }

  if (stabilization_flag == 1) {
    if (must_redefine_groups) {
      group->find_or_create(exclude_group);
      char **newarg;
      newarg = new char*[5];
      newarg[0] = exclude_group;
      newarg[1] = (char *) "dynamic";
      newarg[2] = (char *) "all";
      newarg[3] = (char *) "property";
      newarg[4] = (char *) "statted_tags";
      group->assign(5,newarg);
      delete [] newarg;
    }
  }
}

@@ -1812,7 +1853,7 @@ void FixBondReact::limit_bond(int limit_bond_mode)
  int index1 = atom->find_custom("limit_tags",flag);
  int *i_limit_tags = atom->ivector[index1];

  int index2 = atom->find_custom("statted_tags",flag);
  int index2 = atom->find_custom(statted_id,flag);
  int *i_statted_tags = atom->ivector[index2];

  int index3 = atom->find_custom("react_tags",flag);
@@ -1842,7 +1883,7 @@ void FixBondReact::unlimit_bond()
  int index1 = atom->find_custom("limit_tags",flag);
  int *i_limit_tags = atom->ivector[index1];

  int index2 = atom->find_custom("statted_tags",flag);
  int index2 = atom->find_custom(statted_id,flag);
  int *i_statted_tags = atom->ivector[index2];

  int index3 = atom->find_custom("react_tags",flag);
+3 −0
Original line number Diff line number Diff line
@@ -57,6 +57,7 @@ class FixBondReact : public Fix {
  double **cutsq,*fraction;
  tagint lastcheck;
  int stabilization_flag;
  int custom_exclude_flag;
  int *stabilize_steps_flag;
  int *update_edges_flag;
  int status;
@@ -88,6 +89,8 @@ class FixBondReact : public Fix {
  char *nve_limit_xmax; // indicates max distance allowed to move when relaxing
  char *id_fix1; // id of internally created fix nve/limit
  char *id_fix2; // id of internally created fix per-atom properties
  char *id_fix3; // id of internally created 'stabilization group' per-atom property fix
  char *statted_id; // name of 'stabilization group' per-atom property
  char *master_group; // group containing relaxing atoms from all fix rxns
  char *exclude_group; // group for system-wide thermostat