From 916ed7176072427013204b25ee0ccefc2353cf15 Mon Sep 17 00:00:00 2001
From: jrgissing <jrgiss05@gmail.com>
Date: Sun, 6 May 2018 00:34:03 -0600
Subject: [PATCH] fix bond/react final touches
---
doc/src/fix_bond_react.txt | 75 ++--
.../nylon,6-6_melt/in.large_nylon_melt | 4 +-
.../misc/bond_react/tiny_nylon/in.tiny_nylon | 4 +-
src/USER-MISC/fix_bond_react.cpp | 423 +++++++++++-------
src/USER-MISC/fix_bond_react.h | 9 +-
5 files changed, 306 insertions(+), 209 deletions(-)
diff --git a/doc/src/fix_bond_react.txt b/doc/src/fix_bond_react.txt
index b5f9a73084..13c4f49b46 100644
--- a/doc/src/fix_bond_react.txt
+++ b/doc/src/fix_bond_react.txt
@@ -11,9 +11,9 @@ fix bond/react command :h3
[Syntax:]
fix ID group-ID bond/react common_keyword values ...
- react react-ID react-group-ID Nevery Rmin template-ID(pre-reacted) template-ID(post-reacted) map_file individual_keyword values ...
- react react-ID react-group-ID Nevery Rmin template-ID(pre-reacted) template-ID(post-reacted) map_file individual_keyword values ...
- react react-ID react-group-ID Nevery Rmin template-ID(pre-reacted) template-ID(post-reacted) map_file individual_keyword values ...
+ react react-ID react-group-ID Nevery Rmin Rmax template-ID(pre-reacted) template-ID(post-reacted) map_file individual_keyword values ...
+ react react-ID react-group-ID Nevery Rmin Rmax template-ID(pre-reacted) template-ID(post-reacted) map_file individual_keyword values ...
+ react react-ID react-group-ID Nevery Rmin Rmax template-ID(pre-reacted) template-ID(post-reacted) map_file individual_keyword values ...
... :pre
ID, group-ID are documented in "fix"_fix.html command. Group-ID is ignored. :ulb,l
@@ -29,7 +29,8 @@ react = mandatory argument indicating new reaction specification
react-ID = user-assigned name for the reaction
react-group-ID = only atoms in this group are available for the reaction
Nevery = attempt reaction every this many steps :l
- Rmin = bonding pair atoms separated by less than Rmin can initiate reaction (distance units) :l
+ Rmin = bonding pair atoms must be separated by more than Rmin to initiate reaction (distance units) :l
+ Rmax = bonding pair atoms must be separated by less than Rmax to initiate reaction (distance units) :l
template-ID(pre-reacted) = ID of a molecule template containing pre-reaction topology :l
template-ID(post-reacted) = ID of a molecule template containing post-reaction topology :l
map_file = name of file specifying corresponding atomIDs in the pre- and post-reacted templates :l
@@ -46,15 +47,15 @@ react = mandatory argument indicating new reaction specification
molecule mol1 pre_reacted_topology.txt
molecule mol2 post_reacted_topology.txt
-fix 5 all bond/react stabilization no react myrxn1 all 1 3.25 mol1 mol2 map_file.txt
+fix 5 all bond/react stabilization no react myrxn1 all 1 0 3.25 mol1 mol2 map_file.txt
molecule mol1 pre_reacted_rxn1.txt
molecule mol2 post_reacted_rxn1.txt
molecule mol3 pre_reacted_rxn2.txt
molecule mol4 post_reacted_rxn2.txt
fix 5 all bond/react stabilization yes nvt_grp .03 &
- react myrxn1 all 1 3.25 mol1 mol2 map_file_rxn1.txt prob 0.50 12345 &
- react myrxn2 all 1 2.75 mol3 mol4 map_file_rxn2.txt prob 0.25 12345
+ react myrxn1 all 1 0 3.25 mol1 mol2 map_file_rxn1.txt prob 0.50 12345 &
+ react myrxn2 all 1 0 2.75 mol3 mol4 map_file_rxn2.txt prob 0.25 12345
fix 6 nvt_grp nvt temp 300 300 100 # system-wide thermostat must be defined after bond/react :pre
[Description:]
@@ -101,9 +102,14 @@ The group-ID set using the {stabilization} keyword should be a
previously unused group-ID. The fix bond/react command creates a
"dynamic group"_group.html of this name that excludes reacting atoms.
This dynamic group-ID should then be used by a subsequent system-wide
-time integrator, as shown in the second example above. It is necessary
-to place the time integration command after the fix bond/react command
-due to the internal dynamic grouping performed by fix bond/react.
+time integrator, 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 following comments pertain to each 'react' argument:
@@ -118,21 +124,20 @@ modified to match the post-reaction template.
A bonding atom pair will be identified if several conditions are met.
First, a pair of atoms within the specified react-group-ID of type
-typei and typej must be within a distance Rmin of each other. The atom
-types typei and typej are specified in the pre- and post-reaction
-templates. The distance calculation uses the pair neighbor list,
-therefore bonded neighbor exclusions may prevent a reaction between
-1st, 2nd or 3rd bonded neighbor atoms. If multiple bonding atom pairs
-are identified for an atom, the closest bonding atom partner is set as
-its "nearest" bonding partner. Then, if both an atomi and atomj have
-each other as their nearest bonding partners, these two atoms are
-identified as the bonding atom pair of the reaction site. Once this
-unique bonding atom pair is identified for each reaction, there could
-two or more reactions that involve a given atom on the same timestep.
-If this is the case, only one such reaction is permitted to occur.
-This reaction is chosen randomly from all potential reactions. This
-capability allows e.g. for different reaction pathways to proceed from
-identical reaction sites with user-specified probabilities.
+typei and typej must separated by a distance between Rmin and Rmax. It
+is possible that multiple bonding atom pairs are identified: if the
+bonding atoms in the pre-reacted template are not 1-2, 1-3, or 1-4
+neighbors, the closest bonding atom partner is set as its bonding
+partner; otherwise, the farthest potential partner is chosen. Then, if
+both an atomi and atomj have each other as their nearest bonding
+partners, these two atoms are identified as the bonding atom pair of
+the reaction site. Once this unique bonding atom pair is identified
+for each reaction, there could two or more reactions that involve a
+given atom on the same timestep. If this is the case, only one such
+reaction is permitted to occur. This reaction is chosen randomly from
+all potential reactions. This capability allows e.g. for different
+reaction pathways to proceed from identical reaction sites with
+user-specified probabilities.
The pre-reacted molecule template is specified by a molecule command.
This molecule template file contains a sample reaction site and its
@@ -262,9 +267,11 @@ angles, dihedrals or impropers are supported.
A few capabilities to note: 1) You may specify as many 'react'
arguments as desired. For example, you could break down a complicated
reaction mechanism into several reaction steps, each defined by its
-own 'react' argument. 2) While typically a bond is formed between the
-bonding atom pairs specified in the pre-reacted molecule template,
-this is not required.
+own 'react' argument. 2) While typically a bond is formed or removed
+between the bonding atom pairs specified in the pre-reacted molecule
+template, this is not required. 3) By reversing the order of the pre-
+and post- reacted molecule templates in another 'react' argument, you
+can allow for the possibility of one or more reverse reactions.
The optional keywords deal with the probability of a given reaction
occurring as well as the stable equilibration of each reaction site as
@@ -300,14 +307,14 @@ reaction:
fix 1 bond_react_MASTER_group temp/rescale 1 300 300 10 1
NOTE: This command must be added after the fix bond/react command, and
-will apply to all reaction steps.
+will apply to all reactions.
Computationally, each timestep this fix operates, it loops over
-neighbor lists and computes distances between pairs of atoms in the
-list. It also communicates between neighboring processors to
-coordinate which bonds are created. All of these operations increase
-the cost of a timestep. Thus you should be cautious about invoking
-this fix too frequently.
+neighbor lists (for bond-forming reactions) and computes distances
+between pairs of atoms in the list. It also communicates between
+neighboring processors to coordinate which bonds are created and/or
+removed. All of these operations increase the cost of a timestep. Thus
+you should be cautious about invoking this fix too frequently.
You can dump out snapshots of the current bond topology via the dump
local command.
diff --git a/examples/USER/misc/bond_react/nylon,6-6_melt/in.large_nylon_melt b/examples/USER/misc/bond_react/nylon,6-6_melt/in.large_nylon_melt
index 072a8e3c45..f2dc506dde 100644
--- a/examples/USER/misc/bond_react/nylon,6-6_melt/in.large_nylon_melt
+++ b/examples/USER/misc/bond_react/nylon,6-6_melt/in.large_nylon_melt
@@ -32,8 +32,8 @@ thermo 50
# dump 1 all xyz 100 test_vis.xyz
fix myrxns all bond/react stabilization yes statted_grp .03 &
- react rxn1 all 1 2.9 mol1 mol2 rxn1_stp1_map &
- react rxn2 all 1 5 mol3 mol4 rxn1_stp2_map
+ react rxn1 all 1 0.0 2.9 mol1 mol2 rxn1_stp1_map &
+ react rxn2 all 1 0.0 5.0 mol3 mol4 rxn1_stp2_map
# stable at 800K
fix 1 statted_grp nvt temp 800 800 100
diff --git a/examples/USER/misc/bond_react/tiny_nylon/in.tiny_nylon b/examples/USER/misc/bond_react/tiny_nylon/in.tiny_nylon
index 3e21f69331..1f7e9c42b7 100644
--- a/examples/USER/misc/bond_react/tiny_nylon/in.tiny_nylon
+++ b/examples/USER/misc/bond_react/tiny_nylon/in.tiny_nylon
@@ -33,8 +33,8 @@ thermo 50
# dump 1 all xyz 1 test_vis.xyz
fix myrxns all bond/react stabilization yes statted_grp .03 &
- react rxn1 all 1 2.9 mol1 mol2 rxn1_stp1_map &
- react rxn2 all 1 5 mol3 mol4 rxn1_stp2_map
+ react rxn1 all 1 0.0 2.9 mol1 mol2 rxn1_stp1_map &
+ react rxn2 all 1 0.0 5.0 mol3 mol4 rxn1_stp2_map
fix 1 statted_grp nvt temp 300 300 100
diff --git a/src/USER-MISC/fix_bond_react.cpp b/src/USER-MISC/fix_bond_react.cpp
index b3c92eed67..c3d7866334 100644
--- a/src/USER-MISC/fix_bond_react.cpp
+++ b/src/USER-MISC/fix_bond_react.cpp
@@ -43,6 +43,17 @@ Contributing Author: Jacob Gissinger (jacob.gissinger@colorado.edu)
using namespace LAMMPS_NS;
using namespace FixConst;
+static const char cite_fix_bond_react[] =
+ "fix bond/react:\n\n"
+ "@Article{Gissinger17,\n"
+ " author = {J. R. Gissinger, B. D. Jensen, K. E. Wise},\n"
+ " title = {Modeling chemical reactions in classical molecular dynamics simulations},\n"
+ " journal = {Polymer},\n"
+ " year = 2017,\n"
+ " volume = 128,\n"
+ " pages = {211--217}\n"
+ "}\n\n";
+
#define BIG 1.0e20
#define DELTA 16
#define MAXLINE 256
@@ -62,6 +73,8 @@ enum{ACCEPT,REJECT,PROCEED,CONTINUE,GUESSFAIL,RESTORE};
FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
+ if (lmp->citeme) lmp->citeme->add(cite_fix_bond_react);
+
fix1 = NULL;
fix2 = NULL;
@@ -71,7 +84,6 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
MPI_Comm_size(world,&nprocs);
attempted_rxn = 0;
- ghostcheck_flag = 0;
force_reneighbor = 1;
next_reneighbor = -1;
vector_flag = 1;
@@ -127,7 +139,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
// this looks excessive
// the price of vectorization (all reactions in one command)?
memory->create(nevery,nreacts,"bond/react:nevery");
- memory->create(cutsq,nreacts,"bond/react:cutsq");
+ memory->create(cutsq,nreacts,2,"bond/react:cutsq");
memory->create(unreacted_mol,nreacts,"bond/react:unreacted_mol");
memory->create(reacted_mol,nreacts,"bond/react:reacted_mol");
memory->create(fraction,nreacts,"bond/react:fraction");
@@ -138,6 +150,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
memory->create(jatomtype,nreacts,"bond/react:jatomtype");
memory->create(ibonding,nreacts,"bond/react:ibonding");
memory->create(jbonding,nreacts,"bond/react:jbonding");
+ memory->create(closeneigh,nreacts,"bond/react:closeneigh");
memory->create(groupbits,nreacts,"bond/react:groupbits");
memory->create(reaction_count,nreacts,"bond/react:reaction_count");
memory->create(local_rxn_count,nreacts,"bond/react:local_rxn_count");
@@ -176,7 +189,11 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
double cutoff = force->numeric(FLERR,arg[iarg++]);
if (cutoff < 0.0) error->all(FLERR,"Illegal fix bond/react command 0.5");
- cutsq[rxn] = cutoff*cutoff;
+ cutsq[rxn][0] = cutoff*cutoff;
+
+ cutoff = force->numeric(FLERR,arg[iarg++]);
+ if (cutoff < 0.0) error->all(FLERR,"Illegal fix bond/react command 0.55");
+ cutsq[rxn][1] = cutoff*cutoff;
unreacted_mol[rxn] = atom->find_molecule(arg[iarg++]);
if (unreacted_mol[rxn] == -1) error->all(FLERR,"Unreacted molecule template ID for "
@@ -241,6 +258,23 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
if (atom->molecular != 1)
error->all(FLERR,"Cannot use fix bond/react with non-molecular systems");
+ // check if bonding atoms are 1-2, 1-3, or 1-4 bonded neighbors
+ // if so, we don't need non-bonded neighbor list
+ for (int myrxn = 0; myrxn < nreacts; myrxn++) {
+ closeneigh[myrxn] = -1; // indicates will search non-bonded neighbors
+ onemol = atom->molecules[unreacted_mol[myrxn]];
+ for (int k = 0; k < onemol->nspecial[ibonding[myrxn]-1][2]; k++) {
+ if (onemol->special[ibonding[myrxn]-1][k] == jbonding[myrxn]) {
+ closeneigh[myrxn] = 2; // index for 1-4 neighbor
+ if (k < onemol->nspecial[ibonding[myrxn]-1][1])
+ closeneigh[myrxn] = 1; // index for 1-3 neighbor
+ if (k < onemol->nspecial[ibonding[myrxn]-1][0])
+ closeneigh[myrxn] = 0; // index for 1-2 neighbor
+ break;
+ }
+ }
+ }
+
// initialize Marsaglia RNG with processor-unique seed
random = new class RanMars*[nreacts];
@@ -259,6 +293,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
nmax = 0;
partner = finalpartner = NULL;
distsq = NULL;
+ probability = NULL;
maxcreate = 0;
created = NULL;
local_ncreate = NULL;
@@ -302,6 +337,7 @@ FixBondReact::~FixBondReact()
memory->destroy(local_ncreate);
memory->destroy(ncreate);
memory->destroy(distsq);
+ memory->destroy(probability);
memory->destroy(created);
memory->destroy(edge);
memory->destroy(equivalences);
@@ -320,6 +356,7 @@ FixBondReact::~FixBondReact()
memory->destroy(jatomtype);
memory->destroy(ibonding);
memory->destroy(jbonding);
+ memory->destroy(closeneigh);
memory->destroy(groupbits);
memory->destroy(reaction_count);
memory->destroy(local_rxn_count);
@@ -372,7 +409,7 @@ it will have the name 'i_limit_tags' and will be intitialized to 0 (not in group
void FixBondReact::post_constructor()
{
- //let's add the limit_tags per-atom property fix
+ // let's add the limit_tags per-atom property fix
int len = strlen("per_atom_props") + 1;
id_fix2 = new char[len];
strcpy(id_fix2,"per_atom_props");
@@ -521,7 +558,7 @@ void FixBondReact::init()
// check cutoff for iatomtype,jatomtype
for (int i = 0; i < nreacts; i++) {
- if (force->pair == NULL || cutsq[i] > force->pair->cutsq[iatomtype[i]][jatomtype[i]])
+ if (force->pair == NULL || cutsq[i][1] > force->pair->cutsq[iatomtype[i]][jatomtype[i]])
error->all(FLERR,"Fix bond/react cutoff is longer than pairwise cutoff");
}
@@ -548,10 +585,6 @@ void FixBondReact::init_list(int id, NeighList *ptr)
void FixBondReact::post_integrate()
{
- int inum,jnum,itype,jtype,possible;
- double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
- int *ilist,*jlist,*numneigh,**firstneigh;
-
// check if any reactions could occur on this timestep
int nevery_check = 1;
for (int i = 0; i < nreacts; i++) {
@@ -578,7 +611,6 @@ void FixBondReact::post_integrate()
comm->forward_comm();
// resize bond partner list and initialize it
- // probability array overlays distsq array
// needs to be atom->nmax in length
if (atom->nmax > nmax) {
@@ -587,13 +619,14 @@ void FixBondReact::post_integrate()
memory->destroy(distsq);
memory->destroy(local_ncreate);
memory->destroy(ncreate);
+ memory->destroy(probability);
nmax = atom->nmax;
memory->create(partner,nmax,"bond/react:partner");
memory->create(finalpartner,nmax,"bond/react:finalpartner");
- memory->create(distsq,nmax,"bond/react:distsq");
+ memory->create(distsq,nmax,2,"bond/react:distsq");
memory->create(local_ncreate,nreacts,"bond/react:local_ncreate");
memory->create(ncreate,nreacts,"bond/react:ncreate");
- probability = distsq;
+ memory->create(probability,nmax,"bond/react:probability");
}
// reset create counts
@@ -605,129 +638,28 @@ void FixBondReact::post_integrate()
int nlocal = atom->nlocal;
int nall = atom->nlocal + atom->nghost;
- // NOTE: this might be faster if we remembered neighbor distances from
- // previous timestep and used those --JG
-
// loop over neighbors of my atoms
// each atom sets one closest eligible partner atom ID to bond with
- double **x = atom->x;
tagint *tag = atom->tag;
- int **nspecial = atom->nspecial;
- tagint **special = atom->special;
- int *mask = atom->mask;
int *type = atom->type;
neighbor->build_one(list,1);
- inum = list->inum;
- ilist = list->ilist;
- numneigh = list->numneigh;
- firstneigh = list->firstneigh;
-
- // per-atom property indicating if in bond/react master group
- int flag;
- int index1 = atom->find_custom("limit_tags",flag);
- int *i_limit_tags = atom->ivector[index1];
int i,j;
- for (int myrxn = 0; myrxn < nreacts; myrxn++) {
+ for (rxnID = 0; rxnID < nreacts; rxnID++) {
for (int ii = 0; ii < nall; ii++) {
partner[ii] = 0;
finalpartner[ii] = 0;
- distsq[ii] = BIG;
+ distsq[ii][0] = 0.0;
+ distsq[ii][1] = BIG;
}
- for (int ii = 0; ii < inum; ii++) { // inum vs nlocal
- i = ilist[ii];
- if (!(mask[i] & groupbits[myrxn])) continue;
- if (i_limit_tags[i] != 0) continue;
- itype = type[i];
- xtmp = x[i][0];
- ytmp = x[i][1];
- ztmp = x[i][2];
- jlist = firstneigh[i];
- jnum = numneigh[i];
-
- for (int jj = 0; jj < jnum; jj++) {
- j = jlist[jj];
- j &= NEIGHMASK;
-
- if (!(mask[j] & groupbits[myrxn])) {
- continue;
- }
-
- if (i_limit_tags[j] != 0) {
- continue;
- }
-
- jtype = type[j];
- possible = 0;
-
- if (itype == iatomtype[myrxn] && jtype == jatomtype[myrxn]) {
- possible = 1;
- } else if (itype == jatomtype[myrxn] && jtype == iatomtype[myrxn]) {
- possible = 1;
- }
-
- if (possible == 0) continue;
-
- for (int k = 0; k < nspecial[i][0]; k++)
- if (special[i][k] == tag[j]) possible = 0;
- if (!possible) continue;
-
- // NOTE(for below): certain neighbor list settings prevent 3-cycles anyway!
- // e.g., in my examples, must use kspace command to include 1-3 neighbors for consideration here
-
- // do not allow a three-membered ring to be created (by the new bond)
- // check 1-3 neighbors of atom I
- for (int k = nspecial[i][0]; k < nspecial[i][1]; k++)
- if (special[i][k] == tag[j]) possible = 0;
- if (possible == 0) {
- continue;
- }
-
- // do not allow a four-membered ring to be created (by the new bond)
- // check 1-4 neighbors of atom I -> probably make an option
-
- /*
- for (k = nspecial[i][1]; k < nspecial[i][2]; k++)
- if (special[i][k] == tag[j]) {
- possible = 0;
- }
- if (possible == 0) continue;
- */
-
- // do not allow 5 membered rings -> probably make this an option
-
- /*
- for (k = nspecial[i][1]; k < nspecial[i][2]; k++) {
- for (m = nspecial[atom->map(special[i][k])][1]; m < nspecial[atom->map(special[i][k])][2]; m++) {
- if (special[atom->map(special[i][k])][m] == tag[j]) possible = 0;
- / }
- }
- if (possible == 0) continue;
- */
-
- delx = xtmp - x[j][0];
- dely = ytmp - x[j][1];
- delz = ztmp - x[j][2];
- rsq = delx*delx + dely*dely + delz*delz;
-
- if (rsq >= cutsq[myrxn]) {
- continue;
- }
- if (rsq < distsq[i]) {
- partner[i] = tag[j];
- distsq[i] = rsq;
- }
- if (rsq < distsq[j]) {
- partner[j] = tag[i];
- distsq[j] = rsq;
- }
- }
- }
+ // fork between far and close_partner here
+ if (closeneigh[rxnID] < 0) far_partner();
+ else close_partner();
// reverse comm of distsq and partner
// not needed if newton_pair off since I,J pair was seen by both procs
@@ -739,9 +671,9 @@ void FixBondReact::post_integrate()
// for prob check, generate random value for each atom with a bond partner
// forward comm of partner and random value, so ghosts have it
- if (fraction[myrxn] < 1.0) {
+ if (fraction[rxnID] < 1.0) {
for (int i = 0; i < nlocal; i++)
- if (partner[i]) probability[i] = random[myrxn]->uniform();
+ if (partner[i]) probability[i] = random[rxnID]->uniform();
}
commflag = 2;
@@ -764,11 +696,11 @@ void FixBondReact::post_integrate()
// apply probability constraint using RN for atom with smallest ID
- if (fraction[myrxn] < 1.0) {
+ if (fraction[rxnID] < 1.0) {
if (tag[i] < tag[j]) {
- if (probability[i] >= fraction[myrxn]) continue;
+ if (probability[i] >= fraction[rxnID]) continue;
} else {
- if (probability[j] >= fraction[myrxn]) continue;
+ if (probability[j] >= fraction[rxnID]) continue;
}
}
@@ -780,31 +712,17 @@ void FixBondReact::post_integrate()
if (tag[i] < tag[j]) temp_ncreate++;
}
- local_ncreate[myrxn] = temp_ncreate;
+ local_ncreate[rxnID] = temp_ncreate;
// break loop if no even eligible bonding atoms were found (on any proc)
int some_chance;
MPI_Allreduce(&temp_ncreate,&some_chance,1,MPI_INT,MPI_SUM,world);
- if (!some_chance) {
- continue;
- }
+ if (!some_chance) continue;
// communicate final partner
commflag = 3;
comm->forward_comm_fix(this);
-//obsolete comment block
-/*
- // I think this also simplifies for bond/react.
- // but currently unsure how to adapt it - JG
- // create list of created bonds that influence my owned atoms
- // even if between owned-ghost or ghost-ghost atoms
- // finalpartner is now set for owned and ghost atoms so loop over nall
- // OK if duplicates in created list due to ghosts duplicating owned atoms
- // check J < 0 to insure a created bond to unknown atom is included
- // i.e. a bond partner outside of skin length
-*/
-
// add instance to 'created' only if this processor
// owns the atoms with smaller global ID
// NOTE: we no longer care about ghost-ghost instances as bond/create did
@@ -816,21 +734,21 @@ void FixBondReact::post_integrate()
j = atom->map(finalpartner[i]);
// if (j < 0 || tag[i] < tag[j]) {
if (tag[i] < tag[j]) { //atom->map(std::min(tag[i],tag[j])) <= nlocal &&
- if (ncreate[myrxn] == maxcreate) {
+ if (ncreate[rxnID] == maxcreate) {
maxcreate += DELTA;
// third column of 'created': bond/react integer ID
memory->grow(created,maxcreate,2,nreacts,"bond/react:created");
}
// to ensure types remain in same order
// unnecessary now taken from reaction map file
- if (iatomtype[myrxn] == type[i]) {
- created[ncreate[myrxn]][0][myrxn] = tag[i];
- created[ncreate[myrxn]][1][myrxn] = finalpartner[i];
+ if (iatomtype[rxnID] == type[i]) {
+ created[ncreate[rxnID]][0][rxnID] = tag[i];
+ created[ncreate[rxnID]][1][rxnID] = finalpartner[i];
} else {
- created[ncreate[myrxn]][0][myrxn] = finalpartner[i];
- created[ncreate[myrxn]][1][myrxn] = tag[i];
+ created[ncreate[rxnID]][0][rxnID] = finalpartner[i];
+ created[ncreate[rxnID]][1][rxnID] = tag[i];
}
- ncreate[myrxn]++;
+ ncreate[rxnID]++;
}
}
unlimit_bond(); //free atoms that have been relaxed
@@ -854,7 +772,152 @@ void FixBondReact::post_integrate()
superimpose_algorithm();
// free atoms that have been limited after reacting
unlimit_bond();
+}
+
+/* ----------------------------------------------------------------------
+ Search non-bonded neighbor lists if bonding atoms are not in special list
+------------------------------------------------------------------------- */
+
+void FixBondReact::far_partner()
+{
+ int inum,jnum,itype,jtype,possible;
+ double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+
+ // loop over neighbors of my atoms
+ // each atom sets one closest eligible partner atom ID to bond with
+
+ double **x = atom->x;
+ tagint *tag = atom->tag;
+ int **nspecial = atom->nspecial;
+ tagint **special = atom->special;
+ int *mask = atom->mask;
+ int *type = atom->type;
+
+ inum = list->inum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ // per-atom property indicating if in bond/react master group
+ int flag;
+ int index1 = atom->find_custom("limit_tags",flag);
+ int *i_limit_tags = atom->ivector[index1];
+
+ int i,j;
+
+ for (int ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ if (!(mask[i] & groupbits[rxnID])) continue;
+ if (i_limit_tags[i] != 0) continue;
+ itype = type[i];
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ for (int jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ j &= NEIGHMASK;
+
+ if (!(mask[j] & groupbits[rxnID])) {
+ continue;
+}
+
+ if (i_limit_tags[j] != 0) {
+ continue;
+ }
+
+ jtype = type[j];
+ possible = 0;
+
+ if (itype == iatomtype[rxnID] && jtype == jatomtype[rxnID]) {
+ possible = 1;
+ } else if (itype == jatomtype[rxnID] && jtype == iatomtype[rxnID]) {
+ possible = 1;
+ }
+
+ if (possible == 0) continue;
+
+ // do not allow bonding atoms within special list
+ for (int k = 0; k < nspecial[i][2]; k++)
+ if (special[i][k] == tag[j]) possible = 0;
+ if (!possible) continue;
+
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+
+ if (rsq >= cutsq[rxnID][1] || rsq <= cutsq[rxnID][0]) {
+ continue;
+ }
+ if (rsq < distsq[i][1]) {
+ partner[i] = tag[j];
+ distsq[i][1] = rsq;
+ }
+ if (rsq < distsq[j][1]) {
+ partner[j] = tag[i];
+ distsq[j][1] = rsq;
+ }
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ Slightly simpler to find bonding partner when a close neighbor
+------------------------------------------------------------------------- */
+
+void FixBondReact::close_partner()
+{
+ int n,i1,i2,itype,jtype;
+ double delx,dely,delz,rsq;
+
+ double **x = atom->x;
+ tagint *tag = atom->tag;
+ tagint *type = atom->type;
+ int *mask = atom->mask;
+ int **nspecial = atom->nspecial;
+ int **special = atom->special;
+
+ // per-atom property indicating if in bond/react master group
+ int flag;
+ int index1 = atom->find_custom("limit_tags",flag);
+ int *i_limit_tags = atom->ivector[index1];
+ // loop over special list
+ for (int ii = 0; ii < atom->nlocal; ii++) {
+ itype = type[ii];
+ n = 0;
+ if (closeneigh[rxnID] != 0)
+ n = nspecial[ii][closeneigh[rxnID]-1];
+ for (n; n < nspecial[ii][closeneigh[rxnID]]; n++) {
+ i1 = ii;
+ i2 = atom->map(special[ii][n]);
+ jtype = type[i2];
+ if (!(mask[i1] & groupbits[rxnID])) continue;
+ if (!(mask[i2] & groupbits[rxnID])) continue;
+ if (i_limit_tags[i1] != 0) continue;
+ if (i_limit_tags[i2] != 0) continue;
+ if (itype != iatomtype[rxnID] || jtype != jatomtype[rxnID]) continue;
+
+ delx = x[i1][0] - x[i2][0];
+ dely = x[i1][1] - x[i2][1];
+ delz = x[i1][2] - x[i2][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+ if (rsq >= cutsq[rxnID][1] || rsq <= cutsq[rxnID][0]) continue;
+
+ if (rsq > distsq[i1][0]) {
+ partner[i1] = tag[i2];
+ distsq[i1][0] = rsq;
+ }
+ if (rsq > distsq[i2][0]) {
+ partner[i2] = tag[i1];
+ distsq[i2][0] = rsq;
+ }
+ }
+ }
}
/* ----------------------------------------------------------------------
@@ -1076,8 +1139,8 @@ void FixBondReact::make_a_guess()
if (assigned_count == nfirst_neighs) status = GUESSFAIL;
// check if all neigh atom types are the same between simulation and unreacted mol
- int mol_ntypes[atom->ntypes];
- int lcl_ntypes[atom->ntypes];
+ int *mol_ntypes = new int[atom->ntypes];
+ int *lcl_ntypes = new int[atom->ntypes];
for (int i = 0; i < atom->ntypes; i++) {
mol_ntypes[i] = 0;
@@ -1096,6 +1159,9 @@ void FixBondReact::make_a_guess()
}
}
+ delete [] mol_ntypes;
+ delete [] lcl_ntypes;
+
// okay everything seems to be in order. let's assign some ID pairs!!!
neighbor_loop();
}
@@ -1384,7 +1450,7 @@ void FixBondReact::find_landlocked_atoms(int myrxn)
// if no edge atoms (small reacting molecule), all atoms are landlocked
// we can delete all current topology of landlocked atoms and replace
- // alwasy remove edge atoms from landlocked list
+ // always remove edge atoms from landlocked list
for (int i = 0; i < twomol->natoms; i++) {
if (edge[equivalences[i][1][myrxn]-1][myrxn] == 1) landlocked_atoms[i][myrxn] = 0;
else landlocked_atoms[i][myrxn] = 1;
@@ -1395,7 +1461,6 @@ void FixBondReact::find_landlocked_atoms(int myrxn)
if ((force->dihedral && twomol->dihedralflag) ||
(force->improper && twomol->improperflag)) nspecial_limit = 1;
-
if (nspecial_limit != -1) {
for (int i = 0; i < twomol->natoms; i++) {
for (int j = 0; j < twomol->nspecial[i][nspecial_limit]; j++) {
@@ -1410,11 +1475,21 @@ void FixBondReact::find_landlocked_atoms(int myrxn)
// bad molecule templates check
// if atoms change types, but aren't landlocked, that's bad
-
for (int i = 0; i < twomol->natoms; i++) {
if (twomol->type[i] != onemol->type[equivalences[i][1][myrxn]-1] && landlocked_atoms[i][myrxn] == 0)
error->one(FLERR,"Atom affected by reaction too close to template edge");
}
+
+ // also, if atoms change number of bonds, but aren't landlocked, that could be bad
+ if (me == 0)
+ for (int i = 0; i < twomol->natoms; i++) {
+ if (twomol->nspecial[i][0] != onemol->nspecial[equivalences[i][1][myrxn]-1][0] && landlocked_atoms[i][myrxn] == 0) {
+ char str[128];
+ sprintf(str,"An atom in 'react #%d' changes bond connectivity but not atom type",myrxn+1);
+ error->warning(FLERR,str);
+ break;
+ }
+ }
}
/* ----------------------------------------------------------------------
@@ -1438,7 +1513,8 @@ void FixBondReact::dedup_mega_gloves(int dedup_mode)
dedup_size = global_megasize;
}
- tagint dedup_glove[max_natoms+1][dedup_size];
+ tagint **dedup_glove;
+ memory->create(dedup_glove,max_natoms+1,dedup_size,"bond/react:dedup_glove");
if (dedup_mode == 0) {
for (int i = 0; i < dedup_size; i++) {
@@ -1456,8 +1532,8 @@ void FixBondReact::dedup_mega_gloves(int dedup_mode)
// dedup_mask is size dedup_size and filters reactions that have been deleted
// a value of 1 means this reaction instance has been deleted
- int dedup_mask[dedup_size];
- int dup_list[dedup_size];
+ int *dedup_mask = new int[dedup_size];
+ int *dup_list = new int[dedup_size];
for (int i = 0; i < dedup_size; i++) {
dedup_mask[i] = 0;
@@ -1467,14 +1543,12 @@ void FixBondReact::dedup_mega_gloves(int dedup_mode)
// let's randomly mix up our reaction instances first
// then we can feel okay about ignoring ones we've already deleted (or accepted)
// based off std::shuffle
- // will produce same 'random' sequences on different processors.
- // not really an issue, as reaction lists are independent
+ int *temp_rxn = new int[max_natoms+1];
for (int i = dedup_size-1; i > 0; --i) { //dedup_size
- //choose random entry to swap current one with
- int k = rand() % (i+1);
+ // choose random entry to swap current one with
+ int k = floor(random[0]->uniform()*(i+1));
// swap entries
- int temp_rxn[max_natoms+1];
for (int j = 0; j < max_natoms+1; j++)
temp_rxn[j] = dedup_glove[j][i];
@@ -1483,6 +1557,7 @@ void FixBondReact::dedup_mega_gloves(int dedup_mode)
dedup_glove[j][k] = temp_rxn[j];
}
}
+ delete [] temp_rxn;
for (int i = 0; i < dedup_size; i++) {
if (dedup_mask[i] == 0) {
@@ -1557,6 +1632,10 @@ void FixBondReact::dedup_mega_gloves(int dedup_mode)
}
global_megasize = new_global_megasize;
}
+
+ memory->destroy(dedup_glove);
+ delete [] dedup_mask;
+ delete [] dup_list;
}
/* ----------------------------------------------------------------------
@@ -1743,9 +1822,6 @@ void FixBondReact::ghost_glovecast()
for (int j = 0; j < global_megasize; j++)
global_mega_glove[i][j] = 0;
-
- ghostcheck_flag = 1;
-
if (ghostly_num_mega > 0) {
for (int i = 0; i < max_natoms+1; i++) {
for (int j = 0; j < ghostly_num_mega; j++) {
@@ -1806,7 +1882,8 @@ void FixBondReact::update_everything()
// add check for local atoms as well
int update_num_mega;
- tagint update_mega_glove[max_natoms+1][MAX(local_num_mega,global_megasize)];
+ tagint **update_mega_glove;
+ memory->create(update_mega_glove,max_natoms+1,MAX(local_num_mega,global_megasize),"bond/react:update_mega_glove");
for (int pass = 0; pass < 2; pass++) {
@@ -2217,6 +2294,9 @@ void FixBondReact::update_everything()
}
}
+
+ memory->destroy(update_mega_glove);
+
// 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
@@ -2522,7 +2602,10 @@ int FixBondReact::pack_reverse_comm(int n, int first, double *buf)
for (i = first; i < last; i++) {
buf[m++] = ubuf(partner[i]).d;
- buf[m++] = distsq[i];
+ if (closeneigh[rxnID] < 0)
+ buf[m++] = distsq[i][1];
+ else
+ buf[m++] = distsq[i][0];
}
return m;
}
@@ -2538,10 +2621,16 @@ void FixBondReact::unpack_reverse_comm(int n, int *list, double *buf)
if (commflag != 1) {
for (i = 0; i < n; i++) {
j = list[i];
- if (buf[m+1] < distsq[j]) {
+ if (closeneigh[rxnID] < 0)
+ if (buf[m+1] < distsq[j][1]) {
partner[j] = (tagint) ubuf(buf[m++]).i;
- distsq[j] = buf[m++];
+ distsq[j][1] = buf[m++];
} else m += 2;
+ else
+ if (buf[m+1] > distsq[j][0]) {
+ partner[j] = (tagint) ubuf(buf[m++]).i;
+ distsq[j][0] = buf[m++];
+ } else m += 2;
}
}
}
diff --git a/src/USER-MISC/fix_bond_react.h b/src/USER-MISC/fix_bond_react.h
index a0f52c0288..e60fff37fb 100644
--- a/src/USER-MISC/fix_bond_react.h
+++ b/src/USER-MISC/fix_bond_react.h
@@ -53,7 +53,7 @@ class FixBondReact : public Fix {
FILE *fp;
int *iatomtype,*jatomtype;
int *seed;
- double *cutsq,*fraction;
+ double **cutsq,*fraction;
tagint lastcheck;
int stabilization_flag;
int *stabilize_steps_flag;
@@ -66,7 +66,7 @@ class FixBondReact : public Fix {
int nmax; // max num local atoms
int max_natoms; // max natoms in a molecule template
tagint *partner,*finalpartner;
- double *distsq,*probability;
+ double **distsq,*probability;
int *ncreate;
int maxcreate;
int allncreate;
@@ -96,10 +96,9 @@ class FixBondReact : public Fix {
void superimpose_algorithm(); // main function of the superimpose algorithm
int *ibonding,*jbonding;
+ int *closeneigh; // indicates if bonding atoms of a rxn are 1-2, 1-3, or 1-4 neighbors
int nedge,nequivalent; // number of edge, equivalent atoms in mapping file
int attempted_rxn; // there was an attempt!
- int ghostcheck_flag; // idicates whether a reaction instances contains a nonlocal atom
- int this_rxn_count; // num of local reaction occurrences
int *local_rxn_count;
int *ghostly_rxn_count;
int avail_guesses; // num of restore points available
@@ -143,6 +142,8 @@ class FixBondReact : public Fix {
void skip_lines(int, char *);
int parse(char *, char **, int);
+ void far_partner();
+ void close_partner();
void find_landlocked_atoms(int);
void glove_ghostcheck();
void ghost_glovecast();
--
GitLab