diff --git a/doc/src/fix_bond_react.txt b/doc/src/fix_bond_react.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a1d8f748a0dfa8f5bddba88b96357085379c887d
--- /dev/null
+++ b/doc/src/fix_bond_react.txt
@@ -0,0 +1,206 @@
+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Section_commands.html#comm)
+
+:line
+
+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 ...
+ ... :pre
+
+ID, group-ID are documented in "fix"_fix.html command. Group-ID is ignored. :ulb,l
+bond/react = style name of this fix command :l
+zero or more common keyword/value pairs may be appended directly after 'bond/react' :l
+these apply to all reaction specifications (below) :l
+common_keyword = {stabilization}
+ {stabilization} values = group-ID xmax
+ group-ID = user-assigned ID of an internally-created dynamic group that excludes reacting atoms, and can be used by a subsequent time integration fix such as nvt, npt, or nve (cannot be 'all')
+ {xmax} value = distance
+ distance = xmax value that is used by an internally created "nve/limit"_nve_limit.html integrator
+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
+ 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
+ zero or more individual keyword/value pairs may be appended to each react argument :l
+ individual_keyword = {prob} or {stabilize_steps} :l
+ {prob} values = fraction seed
+ fraction = initiate reaction with this probability if otherwise eligible
+ seed = random number seed (positive integer)
+ {stabilize_steps} value = timesteps
+ timesteps = number of timesteps to apply internally created nve/limit:pre
+:ule
+
+[Examples:]
+
+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
+
+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
+fix 6 nvt_grp nvt temp 300 300 100 # system-wide thermostat must be defined after bond/react :pre
+
+[Description:]
+
+Initiate complex covalent bonding (topology) changes. These topology changes will be referred to as "reactions" throughout this documentation. Topology changes are defined in pre- and post-reaction molecule templates and can include creation and deletion of bonds, angles, dihedrals, impropers, bond-types, angle-types, dihedral-types, atom-types, or atomic charges.
+
+Fix bond/react does not use quantum mechanical (eg. fix qmmm) or pairwise bond-order potential (eg. Tersoff or AIREBO) methods to determine bonding changes a priori. Rather, it uses a distance-based probabilistic criteria to effect predetermined topology changes in simulations using standard force fields.
+
+This fix was created to facilitate the dynamic creation of polymeric, amorphous or highly-crosslinked systems. A suggested workflow for using this fix is: 1) identify a reaction to be simulated 2) build a molecule template of the reaction site before the reaction has occurred 3) build a molecule template of the reaction site after the reaction has occurred 4) create a map that relates the template-atom-IDs of each atom between pre- and post-reaction molecule templates 5) fill a simulation box with molecules and run a simulation with fix/bond react.
+
+Only one 'fix bond/react' command can be used at a time. Multiple reactions can be simultaneously applied by specifying multiple 'react' arguments to a single 'fix bond/react' command. This syntax is necessary because the ‘common keywords’ are applied to all reactions.
+
+The {stabilization} keyword enables reaction site stabilization. Reaction site stabilization is performed by including reacting atoms in an internally created fix "nve/limit"_fix_nve_limit.html time integrator for a set number of timesteps given by the {stabilize_steps} keyword. While reacting atoms are being time integrated by the internal nve/limit, they are prevented from being 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. 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.
+
+The following comments pertain to each 'react' argument:
+
+A check for possible new reaction sites is performed every Nevery timesteps.
+
+Two conditions must be met for a reaction to occur. First a bonding atom pair must be identified. Second, the topology surrounding the bonding atom pair must match the topology of the pre-reaction template. If both these conditions are met, the reaction site is 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.
+
+The pre-reacted molecule template is specified by a molecule command. This molecule template file contains a sample reaction site and its 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.
+
+Some atoms in the pre-reacted template that are not reacting may have missing topology with respect to the simulation. For example, the pre-reacted template may contain an atom that would connect to the rest of a long polymer chain. These are referred to as edge atoms, and are also specified in the map file.
+
+Note that some care must be taken when a building a molecule template for a given simulation. All atom types in the pre-reacted template must be the same as those of a potential reaction site in the simulation. A detailed discussion of matching molecule template atom types with the simulation is provided on the "molecule"_molecule.html command page.
+
+The post-reacted molecule template contains a sample of the reaction site and its surrounding topology after the reaction has occurred. It must contain the same number of atoms as the pre-reacted template. A one-to-one correspondence between the atom IDs in the pre- and post-reacted templates is specified in the map file as described below. Note that during a reaction, an atom, bond, etc. type may change to one that was previously not present in the simulation. These new types must also be defined during the setup of a given simulation. A discussion of correctly handling this is also provided on the "molecule"_molecule.html command page.
+
+The map file is a text document with the following format:
+
+Format of the map file
+
+A map file has a header and a body. The header appears first. The first line of the header is always skipped; it typically contains a description of the file. Lines can have a trailing comment starting with '#' that is ignored. If the line is blank (only whitespace after comment is deleted), it is skipped. If the line contains a header keyword, the corresponding value(s) is read from the line. If it doesn’t contain a header keyword, the line begins the body of the file.
+
+The header contains one mandatory keyword and one optional keyword. The mandatory keyword is 'equivalences' and the optional keyword is 'edgeIDs.' These specify the number of atoms in the pre- and post-reacted templates and the number of edge atoms in pre-reacted template, respectively.
+
+The body contains two mandatory sections and one optional section. The first 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 with the keyword 'Equivalences' and lists a one-to-one correspondence between atom IDs of the pre- and post-reacted templates. The optional section begins with the keyword 'EdgeIDs' and list the atom IDs of edge atoms in the pre-reacted molecule template.
+
+Format of the header of the map file
+
+These are the recognized header keywords. Header lines can come in any order. The value(s) are read from the beginning of the line. Thus the keyword 'equivalences' should be in a line like "25 equivalences."
+
+equivalences = # of atoms in the pre- and post-reacted molecule templates
+edgeIDs = # of edge atoms in the pre-reacted molecule template
+
+The edgeIDs keyword is optional.
+
+Format of the body of the map file
+
+These are the section keywords for the body of the file.
+
+BondingIDs, EdgeIDs = list of atom IDs of bonding and edge atoms in the pre-reacted molecule template
+
+Equivalences = a two column list where the first column is an atom ID of the pre-reacted
+molecule template, and the second column is the corresponding atom ID of the post-reacted
+molecule template
+
+The bondingIDs section will always contain two atom IDs, corresponding to the bonding atom pairs of the pre-reacted map file. The Equivalences section will contain as many rows as there are atoms in the pre- and post-reacted molecule templates. The edgeIDs section is optional, but would contain an atom ID for each edge atom in the pre-reacted molecule template.
+
+A sample map file is given below:
+
+---------------
+this is a map file
+
+2 edgeIDs
+7 equivalences
+
+BondingIDs
+
+3
+5
+
+EdgeIDs
+
+1
+7
+
+Equivalences
+
+1 1
+2 2
+3 3
+4 4
+5 5
+6 6
+7 7
+---------------
+
+Once a reaction site has been successfully identified, data structures within LAMMPS that store bond topology are updated to reflect the post-reacted molecule template. All force fields with fixed bonds, 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 a fix bond/react. 2) While typically a bond is formed 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 template in another fix bond/react command, 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 it occurs.
+
+The {prob} keyword can affect whether an eligible reaction actually occurs. The fraction setting must be a value between 0.0 and 1.0. A uniform random number between 0.0 and 1.0 is generated and the eligible reaction only occurs if the random number is less than the fraction.
+
+The {stabilize_steps} keyword allows for the specification of how many timesteps a reaction site is stabilized before being returned to the overall system thermostat.
+
+In order to produce the most physical behavior, this 'reaction site equilibration time' should be tuned to be as small as possible while retaining stability for a given system or reaction step. After a limited number of case studies, this number has been set to a default of 60 timesteps. Ideally, it should be individually tuned for each fix bond/react. Note that in some situations, decreasing rather than increasing this parameter will result in an increase in stability.
+
+A few other considerations:
+
+It may be beneficial to ensure reacting atoms are at a certain temperature before being released to the overall thermostat. For this, you can use the internally-created dynamic group named "bond_react_MASTER_group." For example, adding the following command would thermostat the group of all atoms currently involved in a 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 specifications.
+
+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.
+
+You can dump out snapshots of the current bond topology via the dump local command.
+
+:line
+
+[Restart, fix_modify, output, run start/stop, minimize info:]
+
+No information about this fix is written to "binary restart files"_restart.html. None of the "fix_modify"_fix_modify.html options are relevant to this fix.
+
+This fix computes one statistic for each 'react' argument that it stores in a global vector, of length 'number of react arguments', that can be accessed by various "output commands"_Section_howto.html#howto_15. The vector values calculated by this fix are "intensive".
+
+These is 1 quantity for each react argument:
+
+(1) cumulative # of reactions occurred :ul
+
+No parameter of this fix can be used with the {start/stop} keywords of the "run"_run.html command. This fix is not invoked during "energy minimization"_minimize.html.
+
+[Restrictions:]
+
+This fix is part of the MC package. It is only enabled if LAMMPS was built with that package. See the "Making LAMMPS"_Section_start.html#start_3 section for more info.
+
+[Related commands:]
+
+"fix bond/create"_fix_bond_create.html,
+"fix bond/break"_fix_bond_break.html, "fix
+bond/swap"_fix_bond_swap.html, "dump local"_dump.html,
+"special_bonds"_special_bonds.html
+
+[Default:]
+
+The option defaults are stabilization = no, stabilize_steps = 60
+
+:line
+
+:link(Gissinger)
+[(Gissinger)] Gissinger, Jensen and Wise, Polymer, 128, 211 (2017).
diff --git a/examples/bond_react/tiny_nylon/infromdata.class2 b/examples/bond_react/tiny_nylon/infromdata.class2
new file mode 100644
index 0000000000000000000000000000000000000000..52f79f8057c3e1222e755e8503cf64037d237ebf
--- /dev/null
+++ b/examples/bond_react/tiny_nylon/infromdata.class2
@@ -0,0 +1,51 @@
+# Jake practice run
+
+units real
+
+boundary p p p
+
+atom_style full
+
+log testing_log
+
+kspace_style pppm 1.0e-2 # actually, this appears required to make LAMMPS to check 1-3 neighbors
+
+pair_style lj/class2/coul/long 8.5
+
+angle_style class2
+
+bond_style class2
+
+dihedral_style class2
+
+improper_style class2
+
+read_data tiny_nylon.data
+
+velocity all create 300.0 4928459 dist gaussian
+
+molecule mol1 rxn1_stp1_unreacted.data_template
+molecule mol2 rxn1_stp1_reacted.data_template
+molecule mol3 rxn1_stp2_unreacted.data_template
+molecule mol4 rxn1_stp2_reacted.data_template
+
+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
+
+fix 1 statted_grp nvt temp 300 300 100
+
+fix 4 bond_react_MASTER_group temp/rescale 1 300 300 10 1
+
+thermo_style custom step temp press density f_myrxns[1] f_myrxns[2]
+
+restart 100 restart1 restart2
+
+run 10000
+
+write_restart restart_longrun
+write_data restart_longrun.data
diff --git a/examples/bond_react/tiny_nylon/rxn1_stp1_map b/examples/bond_react/tiny_nylon/rxn1_stp1_map
new file mode 100644
index 0000000000000000000000000000000000000000..cfbea2e0701f7344bc10212c0997864eb161f682
--- /dev/null
+++ b/examples/bond_react/tiny_nylon/rxn1_stp1_map
@@ -0,0 +1,35 @@
+this is a marvellous superimpose file
+
+2 edgeIDs
+18 equivalences
+
+BondingIDs
+
+10
+1
+
+EdgeIDs
+
+16
+8
+
+Equivalences
+
+1 1
+2 2
+3 3
+4 4
+5 5
+6 6
+7 7
+8 8
+9 9
+10 10
+11 11
+12 12
+13 13
+14 14
+15 15
+16 16
+17 17
+18 18
diff --git a/examples/bond_react/tiny_nylon/rxn1_stp1_reacted.data_template b/examples/bond_react/tiny_nylon/rxn1_stp1_reacted.data_template
new file mode 100644
index 0000000000000000000000000000000000000000..a188cc611f78942e1228dbf90d5e8bcc5128900b
--- /dev/null
+++ b/examples/bond_react/tiny_nylon/rxn1_stp1_reacted.data_template
@@ -0,0 +1,189 @@
+LAMMPS data file. msi2lmp v3.9.6 / 11 Sep 2014 / CGCMM for rxn1_stp1_reacted
+
+18 atoms
+17 bonds
+31 angles
+39 dihedrals
+20 impropers
+
+Types
+
+1 9
+2 1
+3 1
+4 8
+5 8
+6 4
+7 4
+8 1
+9 1
+10 2
+11 6
+12 3
+13 4
+14 4
+15 5
+16 1
+17 4
+18 4
+
+Charges
+
+1 -0.300000
+2 0.000000
+3 0.000000
+4 0.000000
+5 0.000000
+6 0.000000
+7 0.000000
+8 0.000000
+9 0.000000
+10 0.300000
+11 0.000000
+12 0.000000
+13 0.000000
+14 0.000000
+15 0.000000
+16 0.000000
+17 0.000000
+18 0.000000
+
+Coords
+
+1 -5.522237 -0.752722 1.631158
+2 -5.170398 -0.545733 0.178130
+3 -6.469695 -0.553072 -0.648889
+4 -6.052076 -1.721152 1.744648
+5 -6.183059 0.071387 1.971497
+6 -4.489340 -1.389197 -0.173156
+7 -4.637591 0.453703 0.051252
+8 -5.618658 0.138919 4.386107
+9 -4.669492 -0.989819 3.943591
+10 -4.270194 -0.766405 2.474102
+11 -3.348470 -1.875393 2.024289
+12 -3.569794 0.564183 2.345995
+13 -5.201079 -1.993301 4.044219
+14 -3.736682 -0.984819 4.598305
+15 -4.255402 1.370923 2.679069
+16 -6.136394 -0.339866 -2.136775
+17 -6.996331 -1.555519 -0.517408
+18 -7.153308 0.284949 -0.289930
+
+Bonds
+
+1 11 1 2
+2 12 1 4
+3 12 1 5
+4 13 1 10
+5 2 2 3
+6 1 2 6
+7 1 2 7
+8 2 3 16
+9 1 3 17
+10 1 3 18
+11 2 8 9
+12 4 9 10
+13 1 9 13
+14 1 9 14
+15 5 10 11
+16 3 10 12
+17 6 12 15
+
+Angles
+
+1 17 2 1 4
+2 17 2 1 5
+3 18 2 1 10
+4 19 4 1 5
+5 20 4 1 10
+6 20 5 1 10
+7 21 1 2 3
+8 22 1 2 6
+9 22 1 2 7
+10 2 3 2 6
+11 2 3 2 7
+12 1 6 2 7
+13 3 2 3 16
+14 2 2 3 17
+15 2 2 3 18
+16 2 16 3 17
+17 2 16 3 18
+18 1 17 3 18
+19 8 8 9 10
+20 2 8 9 13
+21 2 8 9 14
+22 23 13 9 10
+23 23 14 9 10
+24 1 13 9 14
+25 6 9 10 11
+26 4 9 10 12
+27 24 1 10 9
+28 25 11 10 12
+29 26 1 10 11
+30 27 1 10 12
+31 7 10 12 15
+
+Dihedrals
+
+1 19 4 1 2 3
+2 20 4 1 2 6
+3 20 4 1 2 7
+4 19 5 1 2 3
+5 20 5 1 2 6
+6 20 5 1 2 7
+7 21 10 1 2 3
+8 22 10 1 2 6
+9 22 10 1 2 7
+10 23 2 1 10 9
+11 24 2 1 10 11
+12 25 2 1 10 12
+13 26 4 1 10 9
+14 27 4 1 10 11
+15 28 4 1 10 12
+16 26 5 1 10 9
+17 27 5 1 10 11
+18 28 5 1 10 12
+19 29 1 2 3 16
+20 30 1 2 3 17
+21 30 1 2 3 18
+22 4 16 3 2 6
+23 2 6 2 3 17
+24 2 6 2 3 18
+25 4 16 3 2 7
+26 2 7 2 3 17
+27 2 7 2 3 18
+28 10 8 9 10 11
+29 8 8 9 10 12
+30 31 8 9 10 1
+31 11 13 9 10 11
+32 9 13 9 10 12
+33 32 13 9 10 1
+34 11 14 9 10 11
+35 9 14 9 10 12
+36 32 14 9 10 1
+37 6 9 10 12 15
+38 7 11 10 12 15
+39 33 1 10 12 15
+
+Impropers
+
+1 1 2 1 4 5
+2 1 2 1 4 10
+3 1 2 1 5 10
+4 1 4 1 5 10
+5 1 1 2 3 6
+6 1 1 2 3 7
+7 1 1 2 6 7
+8 1 3 2 6 7
+9 1 2 3 16 17
+10 1 2 3 16 18
+11 1 2 3 17 18
+12 1 16 3 17 18
+13 1 8 9 13 10
+14 1 8 9 14 10
+15 1 8 9 13 14
+16 1 13 9 14 10
+17 1 9 10 11 12
+18 1 1 10 9 11
+19 1 1 10 9 12
+20 1 1 10 11 12
diff --git a/examples/bond_react/tiny_nylon/rxn1_stp1_unreacted.data_template b/examples/bond_react/tiny_nylon/rxn1_stp1_unreacted.data_template
new file mode 100644
index 0000000000000000000000000000000000000000..d97d75a81fa45418c66450909c58580d5b6ca1a8
--- /dev/null
+++ b/examples/bond_react/tiny_nylon/rxn1_stp1_unreacted.data_template
@@ -0,0 +1,160 @@
+LAMMPS data file. msi2lmp v3.9.6 / 11 Sep 2014 / CGCMM for rxn1_stp1_unreacted
+
+18 atoms
+16 bonds
+25 angles
+23 dihedrals
+14 impropers
+
+Types
+
+1 7
+2 1
+3 1
+4 8
+5 8
+6 4
+7 4
+8 1
+9 1
+10 2
+11 6
+12 3
+13 4
+14 4
+15 5
+16 1
+17 4
+18 4
+
+Charges
+
+1 -0.300000
+2 0.000000
+3 0.000000
+4 0.000000
+5 0.000000
+6 0.000000
+7 0.000000
+8 0.000000
+9 0.000000
+10 0.300000
+11 0.000000
+12 0.000000
+13 0.000000
+14 0.000000
+15 0.000000
+16 0.000000
+17 0.000000
+18 0.000000
+
+Coords
+
+1 -4.922858 -0.946982 1.146055
+2 -5.047195 -0.935267 -0.358173
+3 -6.526281 -0.755366 -0.743523
+4 -5.282604 0.020447 1.552710
+5 -3.860697 -1.095850 1.428305
+6 -4.662382 -1.920900 -0.781524
+7 -4.433977 -0.072765 -0.784071
+8 -5.506279 0.202610 4.825816
+9 -4.449177 -0.844592 4.423366
+10 -4.103916 -0.749629 2.925195
+11 -3.376249 -1.886171 2.245643
+12 -4.493235 0.477214 2.137199
+13 -4.849053 -1.888877 4.663994
+14 -3.491823 -0.662913 5.018510
+15 -5.020777 1.189745 2.805427
+16 -3.964987 2.900602 -1.551341
+17 -4.460694 2.836102 0.668882
+18 -4.828494 3.219656 -0.122111
+
+Bonds
+
+1 14 1 2
+2 10 1 4
+3 10 1 5
+4 2 2 3
+5 1 2 6
+6 1 2 7
+7 2 3 16
+8 1 3 17
+9 1 3 18
+10 2 8 9
+11 4 9 10
+12 1 9 13
+13 1 9 14
+14 5 10 11
+15 3 10 12
+16 6 12 15
+
+Angles
+
+1 15 2 1 4
+2 15 2 1 5
+3 16 4 1 5
+4 28 1 2 3
+5 14 1 2 6
+6 14 1 2 7
+7 2 3 2 6
+8 2 3 2 7
+9 1 6 2 7
+10 3 2 3 16
+11 2 2 3 17
+12 2 2 3 18
+13 2 16 3 17
+14 2 16 3 18
+15 1 17 3 18
+16 8 8 9 10
+17 2 8 9 13
+18 2 8 9 14
+19 23 13 9 10
+20 23 14 9 10
+21 1 13 9 14
+22 6 9 10 11
+23 4 9 10 12
+24 25 11 10 12
+25 7 10 12 15
+
+Dihedrals
+
+1 34 4 1 2 3
+2 35 4 1 2 6
+3 35 4 1 2 7
+4 34 5 1 2 3
+5 35 5 1 2 6
+6 35 5 1 2 7
+7 36 1 2 3 16
+8 12 1 2 3 17
+9 12 1 2 3 18
+10 4 16 3 2 6
+11 2 6 2 3 17
+12 2 6 2 3 18
+13 4 16 3 2 7
+14 2 7 2 3 17
+15 2 7 2 3 18
+16 10 8 9 10 11
+17 8 8 9 10 12
+18 11 13 9 10 11
+19 9 13 9 10 12
+20 11 14 9 10 11
+21 9 14 9 10 12
+22 6 9 10 12 15
+23 7 11 10 12 15
+
+Impropers
+
+1 6 2 1 4 5
+2 11 9 10 11 12
+3 1 1 2 3 6
+4 1 1 2 3 7
+5 1 1 2 6 7
+6 1 3 2 6 7
+7 1 2 3 16 17
+8 1 2 3 16 18
+9 1 2 3 17 18
+10 1 16 3 17 18
+11 1 8 9 13 10
+12 1 8 9 14 10
+13 1 8 9 13 14
+14 1 13 9 14 10
diff --git a/examples/bond_react/tiny_nylon/rxn1_stp2_map b/examples/bond_react/tiny_nylon/rxn1_stp2_map
new file mode 100644
index 0000000000000000000000000000000000000000..64589a9f38ce52a97433932c5ebe846f98e54a37
--- /dev/null
+++ b/examples/bond_react/tiny_nylon/rxn1_stp2_map
@@ -0,0 +1,32 @@
+this is a marvellous superimpose file
+
+2 edgeIDs
+15 equivalences
+
+BondingIDs
+
+4
+12
+
+EdgeIDs
+
+8
+3
+
+Equivalences
+
+1 1
+2 2
+3 3
+4 4
+5 5
+6 6
+7 7
+8 8
+9 9
+10 10
+11 11
+12 12
+13 13
+14 14
+15 15
diff --git a/examples/bond_react/tiny_nylon/rxn1_stp2_reacted.data_template b/examples/bond_react/tiny_nylon/rxn1_stp2_reacted.data_template
new file mode 100644
index 0000000000000000000000000000000000000000..5f990d5581a22bc04e2d368417616b6dda985601
--- /dev/null
+++ b/examples/bond_react/tiny_nylon/rxn1_stp2_reacted.data_template
@@ -0,0 +1,131 @@
+LAMMPS data file. msi2lmp v3.9.6 / 11 Sep 2014 / CGCMM for rxn1_stp2_reacted
+
+15 atoms
+13 bonds
+19 angles
+16 dihedrals
+10 impropers
+
+Types
+
+1 9
+2 1
+3 1
+4 10
+5 8
+6 4
+7 4
+8 1
+9 1
+10 2
+11 6
+12 11
+13 4
+14 4
+15 10
+
+Charges
+
+1 -0.300000
+2 0.000000
+3 0.000000
+4 0.410000
+5 0.000000
+6 0.000000
+7 0.000000
+8 0.000000
+9 0.000000
+10 0.300000
+11 0.000000
+12 -0.820000
+13 0.000000
+14 0.000000
+15 0.410000
+
+Coords
+
+1 -4.856280 -1.050468 1.432625
+2 -5.047195 -0.935267 -0.358173
+3 -6.526281 -0.755366 -0.743523
+4 -5.282604 0.020447 1.552710
+5 -3.860697 -1.095850 1.428305
+6 -4.662382 -1.920900 -0.781524
+7 -4.433977 -0.072765 -0.784071
+8 -5.506279 0.202610 4.825816
+9 -4.449177 -0.844592 4.423366
+10 -4.103916 -0.749629 2.925195
+11 -3.376249 -1.886171 2.245643
+12 -4.493235 0.477214 2.137199
+13 -4.849053 -1.888877 4.663994
+14 -3.491823 -0.662913 5.018510
+15 -5.020777 1.189745 2.805427
+
+Bonds
+
+1 11 1 2
+2 12 1 5
+3 13 1 10
+4 2 2 3
+5 1 2 6
+6 1 2 7
+7 15 4 12
+8 2 8 9
+9 4 9 10
+10 1 9 13
+11 1 9 14
+12 5 10 11
+13 15 15 12
+
+Angles
+
+1 17 2 1 5
+2 18 2 1 10
+3 20 5 1 10
+4 21 1 2 3
+5 22 1 2 6
+6 22 1 2 7
+7 2 3 2 6
+8 2 3 2 7
+9 1 6 2 7
+10 8 8 9 10
+11 2 8 9 13
+12 2 8 9 14
+13 23 13 9 10
+14 23 14 9 10
+15 1 13 9 14
+16 6 9 10 11
+17 24 1 10 9
+18 26 1 10 11
+19 29 15 12 4
+
+Dihedrals
+
+1 19 5 1 2 3
+2 20 5 1 2 6
+3 20 5 1 2 7
+4 21 10 1 2 3
+5 22 10 1 2 6
+6 22 10 1 2 7
+7 23 2 1 10 9
+8 24 2 1 10 11
+9 26 5 1 10 9
+10 27 5 1 10 11
+11 10 8 9 10 11
+12 31 8 9 10 1
+13 11 13 9 10 11
+14 32 13 9 10 1
+15 11 14 9 10 11
+16 32 14 9 10 1
+
+Impropers
+
+1 12 2 1 5 10
+2 13 1 10 9 11
+3 1 1 2 3 6
+4 1 1 2 3 7
+5 1 1 2 6 7
+6 1 3 2 6 7
+7 1 8 9 13 10
+8 1 8 9 14 10
+9 1 8 9 13 14
+10 1 13 9 14 10
diff --git a/examples/bond_react/tiny_nylon/rxn1_stp2_unreacted.data_template b/examples/bond_react/tiny_nylon/rxn1_stp2_unreacted.data_template
new file mode 100644
index 0000000000000000000000000000000000000000..59870cc58381ae61070dbaaa1d3c82c8660b499a
--- /dev/null
+++ b/examples/bond_react/tiny_nylon/rxn1_stp2_unreacted.data_template
@@ -0,0 +1,158 @@
+LAMMPS data file. msi2lmp v3.9.6 / 11 Sep 2014 / CGCMM for rxn1_stp2_unreacted
+
+15 atoms
+14 bonds
+25 angles
+30 dihedrals
+16 impropers
+
+Types
+
+1 9
+2 1
+3 1
+4 8
+5 8
+6 4
+7 4
+8 1
+9 1
+10 2
+11 6
+12 3
+13 4
+14 4
+15 5
+
+Charges
+
+1 -0.300000
+2 0.000000
+3 0.000000
+4 0.000000
+5 0.000000
+6 0.000000
+7 0.000000
+8 0.000000
+9 0.000000
+10 0.300000
+11 0.000000
+12 0.000000
+13 0.000000
+14 0.000000
+15 0.000000
+
+Coords
+
+1 -4.922858 -0.946982 1.146055
+2 -5.047195 -0.935267 -0.358173
+3 -6.526281 -0.755366 -0.743523
+4 -5.282604 0.020447 1.552710
+5 -3.860697 -1.095850 1.428305
+6 -4.662382 -1.920900 -0.781524
+7 -4.433977 -0.072765 -0.784071
+8 -5.506279 0.202610 4.825816
+9 -4.449177 -0.844592 4.423366
+10 -4.103916 -0.749629 2.925195
+11 -3.376249 -1.886171 2.245643
+12 -4.493235 0.477214 2.137199
+13 -4.849053 -1.888877 4.663994
+14 -3.491823 -0.662913 5.018510
+15 -5.020777 1.189745 2.805427
+
+Bonds
+
+1 11 1 2
+2 12 1 4
+3 12 1 5
+4 13 1 10
+5 2 2 3
+6 1 2 6
+7 1 2 7
+8 2 8 9
+9 4 9 10
+10 1 9 13
+11 1 9 14
+12 5 10 11
+13 3 10 12
+14 6 12 15
+
+Angles
+
+1 17 2 1 4
+2 17 2 1 5
+3 18 2 1 10
+4 19 4 1 5
+5 20 4 1 10
+6 20 5 1 10
+7 21 1 2 3
+8 22 1 2 6
+9 22 1 2 7
+10 2 3 2 6
+11 2 3 2 7
+12 1 6 2 7
+13 8 8 9 10
+14 2 8 9 13
+15 2 8 9 14
+16 23 13 9 10
+17 23 14 9 10
+18 1 13 9 14
+19 6 9 10 11
+20 4 9 10 12
+21 24 1 10 9
+22 25 11 10 12
+23 26 1 10 11
+24 27 1 10 12
+25 7 10 12 15
+
+Dihedrals
+
+1 19 4 1 2 3
+2 20 4 1 2 6
+3 20 4 1 2 7
+4 19 5 1 2 3
+5 20 5 1 2 6
+6 20 5 1 2 7
+7 21 10 1 2 3
+8 22 10 1 2 6
+9 22 10 1 2 7
+10 23 2 1 10 9
+11 24 2 1 10 11
+12 25 2 1 10 12
+13 26 4 1 10 9
+14 27 4 1 10 11
+15 28 4 1 10 12
+16 26 5 1 10 9
+17 27 5 1 10 11
+18 28 5 1 10 12
+19 10 8 9 10 11
+20 8 8 9 10 12
+21 31 8 9 10 1
+22 11 13 9 10 11
+23 9 13 9 10 12
+24 32 13 9 10 1
+25 11 14 9 10 11
+26 9 14 9 10 12
+27 32 14 9 10 1
+28 6 9 10 12 15
+29 7 11 10 12 15
+30 33 1 10 12 15
+
+Impropers
+
+1 1 2 1 4 5
+2 1 2 1 4 10
+3 1 2 1 5 10
+4 1 4 1 5 10
+5 1 1 2 3 6
+6 1 1 2 3 7
+7 1 1 2 6 7
+8 1 3 2 6 7
+9 1 8 9 13 10
+10 1 8 9 14 10
+11 1 8 9 13 14
+12 1 13 9 14 10
+13 1 9 10 11 12
+14 1 1 10 9 11
+15 1 1 10 9 12
+16 1 1 10 11 12
diff --git a/examples/bond_react/tiny_nylon/tiny_nylon.data b/examples/bond_react/tiny_nylon/tiny_nylon.data
new file mode 100644
index 0000000000000000000000000000000000000000..23c65b3be574bc8d4b1c6e74ac8a273cde13cc1a
--- /dev/null
+++ b/examples/bond_react/tiny_nylon/tiny_nylon.data
@@ -0,0 +1,795 @@
+LAMMPS data file via write_data, version 23 Jun 2017, timestep = 10013
+
+44 atoms
+11 atom types
+42 bonds
+15 bond types
+74 angles
+29 angle types
+100 dihedrals
+36 dihedral types
+44 impropers
+13 improper types
+5 extra bond per atom
+15 extra angle per atom
+15 extra dihedral per atom
+25 extra improper per atom
+25 extra special per atom
+
+-25 25 xlo xhi
+-25 25 ylo yhi
+-25 25 zlo zhi
+
+Masses
+
+1 12.0112
+2 12.0112
+3 15.9994
+4 1.00797
+5 1.00797
+6 15.9994
+7 14.0067
+8 1.00797
+9 14.0067
+10 1.00797
+11 15.9994
+
+Pair Coeffs # lj/class2/coul/cut
+
+1 0.054 4.01
+2 0.12 3.81
+3 0.24 3.535
+4 0.02 2.7
+5 0.013 1.098
+6 0.267 3.3
+7 0.065 4.07
+8 0.013 1.098
+9 0.106 4.07
+10 0.013 1.098
+11 0.26 3.61
+
+Bond Coeffs # class2
+
+1 1.101 345 -691.89 844.6
+2 1.53 299.67 -501.77 679.81
+3 1.3649 368.731 -832.478 1274.02
+4 1.5202 253.707 -423.037 396.9
+5 1.202 851.14 -1918.49 2160.77
+6 0.965 532.506 -1282.9 2004.77
+7 1.53 299.67 -501.77 679.81
+8 1.101 345 -691.89 844.6
+9 1.457 365.805 -699.637 998.484
+10 1.006 466.74 -1073.6 1251.11
+11 1.452 327.166 -547.899 526.5
+12 1.01 462.75 -1053.63 1545.76
+13 1.416 359.159 -558.473 1146.38
+14 1.457 365.805 -699.637 998.484
+15 0.97 563.28 -1428.22 1902.12
+
+Angle Coeffs # class2
+
+1 107.66 39.641 -12.921 -2.4318
+2 110.77 41.453 -10.604 5.129
+3 112.67 39.516 -7.443 -9.5583
+4 123.145 55.5431 -17.2123 0.1348
+5 118.986 98.6813 -22.2485 10.3673
+6 123.145 55.5431 -17.2123 0.1348
+7 111.254 53.5303 -11.8454 -11.5405
+8 108.53 51.9747 -9.4851 -10.9985
+9 107.734 40.6099 -28.8121 0
+10 110.77 41.453 -10.604 5.129
+11 112.67 39.516 -7.443 -9.5583
+12 107.66 39.641 -12.921 -2.4318
+13 111.91 60.7147 -13.3366 -13.0785
+14 110.62 51.3137 -6.7198 -2.6003
+15 110.954 50.8652 -4.4522 -10.0298
+16 107.067 45.252 -7.5558 -9.512
+17 113.868 45.9271 -20.0824 0
+18 111.037 31.8958 -6.6942 -6.837
+19 116.94 37.5749 -8.6676 0
+20 117.961 37.4964 -8.1837 0
+21 114.302 42.6589 -10.5464 -9.3243
+22 108.937 57.401 2.9374 0
+23 107.734 40.6099 -28.8121 0
+24 116.926 39.4193 -10.9945 -8.7733
+25 118.986 98.6813 -22.2485 10.3673
+26 125.542 92.572 -34.48 -11.1871
+27 0 0 0 0
+28 111.91 60.7147 -13.3366 -13.0785
+29 103.7 49.84 -11.6 -8
+
+BondBond Coeffs
+
+1 5.3316 1.101 1.101
+2 3.3872 1.53 1.101
+3 0 1.53 1.53
+4 0 1.5202 1.3649
+5 0 1.3649 1.202
+6 46.0685 1.5202 1.202
+7 0 1.3649 0.965
+8 5.4199 1.53 1.5202
+9 0.7115 1.5202 1.101
+10 3.3872 1.53 1.101
+11 0 1.53 1.53
+12 5.3316 1.101 1.101
+13 4.6217 1.53 1.457
+14 12.426 1.457 1.101
+15 -6.4168 1.457 1.006
+16 -1.8749 1.006 1.006
+17 -3.471 1.452 1.01
+18 12.1186 1.452 1.416
+19 -0.5655 1.01 1.01
+20 -4.3126 1.01 1.416
+21 3.5446 1.452 1.53
+22 15.2994 1.452 1.101
+23 0.7115 1.101 1.5202
+24 0 1.416 1.5202
+25 0 1.202 1.3649
+26 138.495 1.416 1.202
+27 0 1.416 1.3649
+28 4.6217 1.457 1.53
+29 -9.5 0.97 0.97
+
+BondAngle Coeffs
+
+1 18.103 18.103 1.101 1.101
+2 20.754 11.421 1.53 1.101
+3 8.016 8.016 1.53 1.53
+4 0 0 1.5202 1.3649
+5 0 0 1.3649 1.202
+6 34.9982 37.1298 1.5202 1.202
+7 0 0 1.3649 0.965
+8 18.1678 15.8758 1.53 1.5202
+9 12.4632 9.1765 1.5202 1.101
+10 20.754 11.421 1.53 1.101
+11 8.016 8.016 1.53 1.53
+12 18.103 18.103 1.101 1.101
+13 6.0876 16.5702 1.53 1.457
+14 42.4332 13.4582 1.457 1.101
+15 31.8096 20.5799 1.457 1.006
+16 28.0322 28.0322 1.006 1.006
+17 11.8828 5.9339 1.452 1.01
+18 3.7812 14.8633 1.452 1.416
+19 19.8125 19.8125 1.01 1.01
+20 10.8422 29.5743 1.01 1.416
+21 4.6031 -5.479 1.452 1.53
+22 34.8907 10.6917 1.452 1.101
+23 9.1765 12.4632 1.101 1.5202
+24 0 0 1.416 1.5202
+25 0 0 1.202 1.3649
+26 62.7124 52.4045 1.416 1.202
+27 0 0 1.416 1.3649
+28 16.5702 6.0876 1.457 1.53
+29 22.35 22.35 0.97 0.97
+
+Dihedral Coeffs # class2
+
+1 -0.0228 0 0.028 0 -0.1863 0
+2 -0.1432 0 0.0617 0 -0.1083 0
+3 0.0972 0 0.0722 0 -0.2581 0
+4 0 0 0.0316 0 -0.1681 0
+5 0 0 0.0514 0 -0.143 0
+6 0 0 0 0 0 0
+7 -2.7332 0 2.9646 0 -0.0155 0
+8 0 0 0 0 0 0
+9 0 0 0 0 0 0
+10 0.0442 0 0.0292 0 0.0562 0
+11 -0.1804 0 0.0012 0 0.0371 0
+12 -0.2428 0 0.4065 0 -0.3079 0
+13 -0.1432 0 0.0617 0 -0.1083 0
+14 0.1764 0 0.1766 0 -0.5206 0
+15 0 0 0.0316 0 -0.1681 0
+16 0 0 0.0514 0 -0.143 0
+17 -1.1506 0 -0.6344 0 -0.1845 0
+18 -0.5187 0 -0.4837 0 -0.1692 0
+19 -0.0483 0 -0.0077 0 -0.0014 0
+20 -0.0148 0 -0.0791 0 -0.0148 0
+21 0.0143 0 -0.0132 0 0.0091 0
+22 0.0219 0 -0.026 0 0.0714 0
+23 -0.7532 0 2.7392 0 0.0907 0
+24 0.8297 0 3.7234 0 -0.0495 0
+25 0 0 0 0 0 0
+26 0 0 0 0 0 0
+27 -1.6938 0 2.7386 0 -0.336 0
+28 0 0 0 0 0 0
+29 0.0972 0 0.0722 0 -0.2581 0
+30 -0.0228 0 0.028 0 -0.1863 0
+31 0.1693 0 -0.009 0 -0.0687 0
+32 0.1693 0 -0.009 0 -0.0687 0
+33 0 0 0 0 0 0
+34 -1.1506 0 -0.6344 0 -0.1845 0
+35 -0.5187 0 -0.4837 0 -0.1692 0
+36 0.1764 0 0.1766 0 -0.5206 0
+
+AngleAngleTorsion Coeffs
+
+1 -5.3624 108.53 110.77
+2 -12.564 110.77 110.77
+3 -0.3801 112.67 108.53
+4 -16.164 112.67 110.77
+5 -22.045 112.67 112.67
+6 0 0 111.254
+7 0 118.985 111.254
+8 0 108.53 0
+9 0 107.734 0
+10 -8.019 108.53 123.145
+11 -15.3496 107.734 123.145
+12 -15.7572 111.91 110.77
+13 -12.564 110.77 110.77
+14 -27.3953 112.67 111.91
+15 -16.164 112.67 110.77
+16 -22.045 112.67 112.67
+17 -7.5499 111.91 110.954
+18 -10.4258 110.62 110.954
+19 -4.6337 113.868 114.302
+20 -6.659 113.868 108.937
+21 -7.4314 111.037 114.302
+22 -8.1335 111.037 108.937
+23 -6.5335 111.037 116.926
+24 -15.5547 111.037 125.542
+25 0 111.037 0
+26 -1.3234 117.961 116.926
+27 -7.3186 117.961 125.542
+28 0 117.961 0
+29 -1.0631 114.302 112.67
+30 -12.7974 114.302 110.77
+31 -5.4514 108.53 116.926
+32 -12.2417 107.734 116.926
+33 0 0 111.254
+34 -7.5499 110.954 111.91
+35 -10.4258 110.954 110.62
+36 -27.3953 111.91 112.67
+
+EndBondTorsion Coeffs
+
+1 -0.0204 0.3628 -0.4426 -0.0097 -0.0315 -0.0755 1.5202 1.101
+2 0.213 0.312 0.0777 0.213 0.312 0.0777 1.101 1.101
+3 0.0062 -0.0002 0.0036 0.0055 0.006 -0.0009 1.53 1.5202
+4 0.2486 0.2422 -0.0925 0.0814 0.0591 0.2219 1.53 1.101
+5 -0.0732 0 0 -0.0732 0 0 1.53 1.53
+6 0 0 0 0 0 0 1.5202 0.965
+7 0 0 0 0 0 0 1.202 0.965
+8 0 0 0 0 0 0 1.53 1.3649
+9 0 0 0 0 0 0 1.101 1.3649
+10 0.2654 0.0503 0.1046 -0.281 0.0816 -0.1522 1.53 1.202
+11 1.2143 0.2831 0.3916 -0.2298 0.0354 0.3853 1.101 1.202
+12 0.1022 0.209 0.6433 0.196 0.7056 0.112 1.457 1.101
+13 0.213 0.312 0.0777 0.213 0.312 0.0777 1.101 1.101
+14 0.1032 0.5896 -0.4836 0.0579 -0.0043 -0.1906 1.53 1.457
+15 0.2486 0.2422 -0.0925 0.0814 0.0591 0.2219 1.53 1.101
+16 -0.0732 0 0 -0.0732 0 0 1.53 1.53
+17 -0.9466 0.9356 -0.5542 0.057 0.0625 0.4112 1.53 1.006
+18 -1.1685 0.9266 -0.0993 0.085 0.3061 0.2104 1.101 1.006
+19 -0.0992 -0.0727 -0.4139 0.132 0.0015 0.1324 1.01 1.53
+20 -0.4894 0.1644 0.3105 -0.8983 0.2826 0.0881 1.01 1.101
+21 -0.1245 -0.9369 0.7781 -0.2033 0.0035 0.056 1.416 1.53
+22 0.2292 1.1732 -0.058 -0.3667 0.8197 0.1335 1.416 1.101
+23 0.2299 -0.1141 -0.1424 0.0933 -0.4631 0.2883 1.452 1.5202
+24 0.1598 0.7253 -0.1007 0.1226 -2.1326 0.5581 1.452 1.202
+25 0 0 0 0 0 0 1.452 1.3649
+26 0.6413 0.1676 0.144 -0.6979 0.5619 0.4212 1.01 1.5202
+27 0.1214 0.1936 0.0816 -0.7604 -2.6431 1.2467 1.01 1.202
+28 0 0 0 0 0 0 1.01 1.3649
+29 -0.0797 -0.0406 0.0255 0.0742 0.0105 0.0518 1.452 1.53
+30 0.3022 0.2513 0.4641 -0.0601 -0.3763 -0.1876 1.452 1.101
+31 -0.2631 -0.0076 -0.1145 -0.2751 -0.3058 -0.1767 1.53 1.416
+32 -0.0268 0.7836 0.0035 0.3552 -0.2685 0.5834 1.101 1.416
+33 0 0 0 0 0 0 1.416 0.965
+34 0.057 0.0625 0.4112 -0.9466 0.9356 -0.5542 1.006 1.53
+35 0.085 0.3061 0.2104 -1.1685 0.9266 -0.0993 1.006 1.101
+36 0.0579 -0.0043 -0.1906 0.1032 0.5896 -0.4836 1.457 1.53
+
+MiddleBondTorsion Coeffs
+
+1 -3.5039 1.2458 -0.761 1.53
+2 -14.261 -0.5322 -0.4864 1.53
+3 -1.5945 0.2267 -0.6911 1.53
+4 -14.879 -3.6581 -0.3138 1.53
+5 -17.787 -7.1877 0 1.53
+6 0 0 0 1.3649
+7 0 0 0 1.3649
+8 0 0 0 1.5202
+9 0 0 0 1.5202
+10 0.3388 -0.1096 0.1219 1.5202
+11 0.2359 0.9139 0.9594 1.5202
+12 -10.4959 -0.7647 -0.0545 1.53
+13 -14.261 -0.5322 -0.4864 1.53
+14 -15.4174 -7.3055 -1.0749 1.53
+15 -14.879 -3.6581 -0.3138 1.53
+16 -17.787 -7.1877 0 1.53
+17 -2.2208 0.5479 -0.3527 1.457
+18 -3.4611 1.6996 -0.6007 1.457
+19 -3.5406 -3.3866 0.0352 1.452
+20 -1.1752 2.8058 0.8083 1.452
+21 -3.9501 -0.4002 -0.6798 1.452
+22 -0.6899 -2.2646 1.1579 1.452
+23 0 0 0 1.416
+24 -8.8301 14.3079 -1.7716 1.416
+25 0 0 0 1.416
+26 0 0 0 1.416
+27 -0.9084 6.1447 -0.4852 1.416
+28 0 0 0 1.416
+29 -4.2324 -3.3023 -1.3244 1.53
+30 -4.1028 -0.5941 -0.047 1.53
+31 0 0 0 1.5202
+32 0 0 0 1.5202
+33 0 0 0 1.3649
+34 -2.2208 0.5479 -0.3527 1.457
+35 -3.4611 1.6996 -0.6007 1.457
+36 -15.4174 -7.3055 -1.0749 1.53
+
+BondBond13 Coeffs
+
+1 0 1.5202 1.101
+2 0 1.101 1.101
+3 0 1.53 1.5202
+4 0 1.53 1.101
+5 0 1.53 1.53
+6 0 1.5202 0.965
+7 0 1.202 0.965
+8 0 1.53 1.3649
+9 0 1.101 1.3649
+10 0 1.53 1.202
+11 0 1.101 1.202
+12 0 1.457 1.101
+13 0 1.101 1.101
+14 0 1.53 1.457
+15 0 1.53 1.101
+16 0 1.53 1.53
+17 0 1.53 1.006
+18 0 1.101 1.006
+19 0 1.01 1.53
+20 0 1.01 1.101
+21 0 1.416 1.53
+22 0 1.416 1.101
+23 0 1.452 1.5202
+24 0 1.452 1.202
+25 0 1.452 1.3649
+26 0 1.01 1.5202
+27 0 1.01 1.202
+28 0 1.01 1.3649
+29 0 1.452 1.53
+30 0 1.452 1.101
+31 0 1.53 1.416
+32 0 1.101 1.416
+33 0 1.416 0.965
+34 0 1.006 1.53
+35 0 1.006 1.101
+36 0 1.457 1.53
+
+AngleTorsion Coeffs
+
+1 -0.7466 -0.9448 -0.6321 0.0162 1.4211 -1.4092 108.53 110.77
+2 -0.8085 0.5569 -0.2466 -0.8085 0.5569 -0.2466 110.77 110.77
+3 -0.2607 0.3203 -0.2283 0.0515 -0.0674 -0.0474 112.67 108.53
+4 -0.2454 0 -0.1136 0.3113 0.4516 -0.1988 112.67 110.77
+5 0.3886 -0.3139 0.1389 0.3886 -0.3139 0.1389 112.67 112.67
+6 0 0 0 0 0 0 0 111.254
+7 0 0 0 0 0 0 118.985 111.254
+8 0 0 0 0 0 0 108.53 0
+9 0 0 0 0 0 0 107.734 0
+10 0.0885 -1.3703 -0.5452 0.675 0.5965 0.6725 108.53 123.145
+11 9.1299 -0.4847 0.3582 -1.4946 0.7308 -0.2083 107.734 123.145
+12 -1.1075 0.282 0.8318 0.5111 1.6328 -1.0155 111.91 110.77
+13 -0.8085 0.5569 -0.2466 -0.8085 0.5569 -0.2466 110.77 110.77
+14 -1.9225 -1.345 0.221 2.0125 0.944 -2.7612 112.67 111.91
+15 -0.2454 0 -0.1136 0.3113 0.4516 -0.1988 112.67 110.77
+16 0.3886 -0.3139 0.1389 0.3886 -0.3139 0.1389 112.67 112.67
+17 -3.343 4.4558 -0.0346 0.2873 -0.8072 -0.096 111.91 110.954
+18 -3.9582 2.0063 0.3213 -0.4294 -0.4442 -0.6141 110.62 110.954
+19 -0.5807 0.2041 -0.1384 -2.8967 2.7084 -0.0375 113.868 114.302
+20 -0.3868 0.2041 0.0445 -3.7022 1.3876 0.2393 113.868 108.937
+21 -1.523 1.1296 0.7167 -0.7555 0.0564 1.2177 111.037 114.302
+22 0.0372 -0.3418 -0.0775 -1.5157 2.0781 0.5364 111.037 108.937
+23 5.916 1.7856 0.4052 4.2133 2.9302 3.2903 111.037 116.926
+24 7.4427 2.1505 -0.2206 4.4466 4.0317 1.7129 111.037 125.542
+25 0 0 0 0 0 0 111.037 0
+26 1.9306 0.2105 0.0557 -2.2134 1.2909 0.9726 117.961 116.926
+27 2.3848 0.703 0.1399 -2.6238 0.3606 0.5474 117.961 125.542
+28 0 0 0 0 0 0 117.961 0
+29 0.2039 0.1602 -0.7946 -0.5501 -1.6982 0.2485 114.302 112.67
+30 -1.982 0.2325 -0.3928 -1.2469 1.6933 -1.2081 114.302 110.77
+31 2.1802 -0.0335 -1.3816 2.1221 0.5032 -0.0767 108.53 116.926
+32 7.095 0.0075 0.691 2.0013 0.5068 0.8406 107.734 116.926
+33 0 0 0 0 0 0 0 111.254
+34 0.2873 -0.8072 -0.096 -3.343 4.4558 -0.0346 110.954 111.91
+35 -0.4294 -0.4442 -0.6141 -3.9582 2.0063 0.3213 110.954 110.62
+36 2.0125 0.944 -2.7612 -1.9225 -1.345 0.221 111.91 112.67
+
+Improper Coeffs # class2
+
+1 0 0
+2 0 0
+3 0 0
+4 0 0
+5 0 0
+6 0 0
+7 0 0
+8 0 0
+9 0 0
+10 0 0
+11 0 0
+12 0 0
+13 24.3329 0
+
+AngleAngle Coeffs
+
+1 0 0 0 0 118.985 123.145
+2 0.2738 -0.4825 0.2738 110.77 107.66 110.77
+3 -1.3199 -1.3199 0.1184 112.67 110.77 110.77
+4 2.0403 -1.8202 1.0827 108.53 107.734 110.77
+5 -3.3867 -3.4976 -3.3867 107.734 107.66 107.734
+6 0 0 0 110.954 107.067 110.954
+7 0.2738 -0.4825 0.2738 110.77 107.66 110.77
+8 -1.3199 -1.3199 0.1184 112.67 110.77 110.77
+9 -2.5301 0.5381 2.4286 111.91 110.62 110.77
+10 2.4321 -3.5496 2.4321 110.62 107.66 110.62
+11 0 0 0 123.145 118.985 0
+12 0 0 0 113.868 117.961 111.037
+13 0 0 0 116.926 123.145 125.542
+
+Atoms # full
+
+1 1 1 0.0000000000000000e+00 12.288168 0.738732 4.374280 0 0 0
+2 1 2 2.9999999999999999e-01 13.959928 -0.883144 5.090597 0 0 0
+3 1 3 0.0000000000000000e+00 14.411288 -1.994419 5.682160 0 0 0
+4 1 4 0.0000000000000000e+00 12.881083 0.872503 3.506176 0 0 0
+5 1 4 0.0000000000000000e+00 11.232775 0.801641 3.998777 0 0 0
+6 1 5 0.0000000000000000e+00 13.704366 -2.470396 6.130105 0 0 0
+7 1 1 0.0000000000000000e+00 12.489752 -0.793693 4.710639 0 0 0
+8 1 1 0.0000000000000000e+00 12.455071 1.866388 5.385870 0 0 0
+9 1 1 0.0000000000000000e+00 11.248961 1.901849 6.347664 0 0 0
+10 1 2 2.9999999999999999e-01 10.005971 2.466710 5.772840 -1 1 0
+11 1 6 0.0000000000000000e+00 14.795360 -0.034436 4.807367 0 0 0
+12 1 6 0.0000000000000000e+00 9.115239 1.654547 5.617002 -1 0 0
+13 1 3 0.0000000000000000e+00 9.745096 3.807654 5.573585 -1 1 0
+14 1 4 0.0000000000000000e+00 12.248215 -1.371492 3.808598 0 0 0
+15 1 4 0.0000000000000000e+00 11.715755 -1.036825 5.500449 0 0 0
+16 1 4 0.0000000000000000e+00 12.559724 2.807687 4.858452 0 1 0
+17 1 4 0.0000000000000000e+00 13.299968 1.616570 6.123781 0 0 0
+18 1 4 0.0000000000000000e+00 11.650505 2.330454 7.282410 0 1 0
+19 1 4 0.0000000000000000e+00 10.888420 0.913219 6.637162 -1 0 0
+20 1 5 0.0000000000000000e+00 10.550073 4.294209 5.758192 -1 1 0
+21 2 1 0.0000000000000000e+00 5.851425 1.929552 6.038335 0 0 0
+22 2 1 0.0000000000000000e+00 6.741509 3.160751 6.233074 0 0 0
+23 2 7 -2.9999999999999999e-01 7.957761 3.121780 5.252257 1 0 0
+24 2 7 -2.9999999999999999e-01 2.599653 -2.258940 5.985863 0 -1 0
+25 2 1 0.0000000000000000e+00 3.834337 -1.907078 5.441528 0 -1 0
+26 2 1 0.0000000000000000e+00 4.810793 -1.083699 6.310184 0 -1 0
+27 2 4 0.0000000000000000e+00 6.505912 1.182799 5.449104 0 0 0
+28 2 4 0.0000000000000000e+00 5.156429 2.256468 5.348423 0 0 0
+29 2 4 0.0000000000000000e+00 7.232782 3.178785 7.181911 0 0 0
+30 2 4 0.0000000000000000e+00 6.251671 4.103621 6.222913 0 0 0
+31 2 8 0.0000000000000000e+00 8.249909 4.070668 4.881297 1 0 0
+32 2 8 0.0000000000000000e+00 7.813025 2.623184 4.400744 1 0 0
+33 2 8 0.0000000000000000e+00 2.626695 -2.857547 6.817247 0 -1 0
+34 2 8 0.0000000000000000e+00 1.955281 -2.684319 5.328460 0 -1 0
+35 2 4 0.0000000000000000e+00 3.637708 -1.322842 4.469265 0 -1 0
+36 2 4 0.0000000000000000e+00 4.415570 -2.739689 4.997336 0 -1 0
+37 2 4 0.0000000000000000e+00 5.710714 -1.010014 5.642798 0 -1 0
+38 2 4 0.0000000000000000e+00 5.103831 -1.696423 7.160345 0 -1 0
+39 2 1 0.0000000000000000e+00 5.270763 1.286629 7.308822 0 0 0
+40 2 4 0.0000000000000000e+00 4.834381 2.168531 7.931687 0 0 1
+41 2 4 0.0000000000000000e+00 6.118354 0.786724 7.794709 0 0 1
+42 2 1 0.0000000000000000e+00 4.273849 0.167695 6.957862 0 -1 0
+43 2 4 0.0000000000000000e+00 3.792544 -0.081782 7.904418 0 -1 1
+44 2 4 0.0000000000000000e+00 3.527495 0.674238 6.348869 0 0 0
+
+Velocities
+
+1 -2.4626989626218821e-03 -1.5920230003311222e-03 -3.0621927786115238e-03
+2 9.5082416704385837e-03 -6.9903166167507250e-03 1.3702671335945608e-02
+3 2.3431518493187576e-03 -2.9261683108242173e-03 1.4269399726982105e-03
+4 -1.8184451408256214e-02 3.1103803691687960e-02 -1.3358827768357973e-02
+5 2.6084132471017967e-02 -1.0819576493517332e-02 3.0403384454794881e-02
+6 -4.7312115958218744e-03 -1.9111462399478338e-02 -3.6793354156497558e-02
+7 -7.5068797595949869e-03 6.5661422055962489e-03 1.3226575122695422e-03
+8 3.3807881380161281e-03 3.0458732663557089e-03 2.2368826795446284e-03
+9 -3.1113905793879316e-03 8.2908867720754773e-03 -1.7561238039496530e-03
+10 2.4685206571693056e-03 1.3194776209841030e-03 -2.8041877032800441e-03
+11 -3.4945605770565296e-03 3.2323777135621814e-03 1.6223017668450866e-03
+12 -6.1153483612847778e-03 -5.1534857074262185e-03 1.7735747357354274e-03
+13 2.1384296781859011e-04 -4.5398902942729667e-03 6.1649769894413760e-03
+14 2.5004619864373401e-03 -1.5709184283264888e-03 2.0837548254667757e-02
+15 6.0547939205643532e-03 -1.2650704436910937e-02 -5.4430753266962190e-03
+16 -1.0374605775698001e-02 9.1408658463889240e-03 -1.1306875858287088e-02
+17 -1.2736499128987409e-02 -9.1726811852506501e-03 5.1136502685461254e-03
+18 7.6741778607048112e-03 1.8629856635459279e-02 -1.1300096447670932e-02
+19 -1.8616138775281121e-02 1.0848388547730185e-03 -5.7118433687798576e-03
+20 5.4137572241479059e-03 -1.4564578166395727e-02 -1.2618420441909540e-02
+21 5.8473521452312256e-03 -4.0595286000332086e-03 -6.2517801580146415e-03
+22 3.6402033824753104e-03 -1.4629540504663154e-03 -4.0030712318898046e-03
+23 9.0266305019107689e-03 -2.7511425384659687e-03 4.5576402565437142e-03
+24 -1.3102302415548614e-02 -4.7286703965305791e-03 -1.8966887841189517e-03
+25 7.8621682621103171e-03 -4.2046313540949568e-03 9.6887957374751301e-04
+26 -4.7380176438337968e-03 9.6090441940775827e-03 -8.7592431387039336e-03
+27 5.4311658811632517e-03 2.0032224663495989e-02 -9.4952076489808503e-03
+28 -2.9056381493904374e-03 3.3317109723156875e-03 1.6650350064426677e-02
+29 -6.4569944033489122e-03 2.8423983541959541e-03 -2.6066912906505167e-02
+30 -2.2173867823429387e-02 1.4628839880961319e-02 -2.3330833961402380e-02
+31 9.1925713381983114e-03 -2.5697556639281928e-03 -1.2822203161488303e-02
+32 -8.3206975051927905e-03 -2.2538429924858707e-03 7.7620244118580314e-03
+33 1.9920685674825727e-02 5.0317764848494097e-03 -2.1106672824976403e-02
+34 1.4118463330250982e-02 1.7455545466840316e-02 -1.2482101375598437e-02
+35 -6.1116505640437966e-03 1.3353021777303568e-02 -2.5492434283827668e-02
+36 9.1001521565859649e-03 5.5737774505222404e-03 1.4573768978939985e-02
+37 1.6523593470528035e-03 -2.2107518020000917e-02 2.0311423445130115e-02
+38 -1.0346275393471860e-02 1.6055856586351790e-02 5.5489127019262424e-03
+39 -3.2054811383248638e-03 1.6779208962376315e-03 2.9390509537535661e-03
+40 1.9649219364916443e-02 4.0815776523222859e-03 -9.8422441166041274e-03
+41 5.6961697588160361e-04 7.1361132234741477e-04 4.6335764220256257e-03
+42 2.2221300208006252e-03 3.6217319632558197e-03 -6.3299398503455151e-03
+43 2.5710172734841170e-03 8.0029179814482924e-03 1.9992986928468189e-02
+44 -6.0827581822674656e-03 -1.1834273655641976e-02 2.0526923045885208e-02
+
+Bonds
+
+1 1 1 5
+2 1 1 4
+3 2 1 7
+4 2 1 8
+5 3 2 3
+6 5 2 11
+7 6 3 6
+8 4 7 2
+9 1 7 14
+10 1 7 15
+11 2 8 9
+12 1 8 16
+13 1 8 17
+14 4 9 10
+15 1 9 18
+16 1 9 19
+17 5 10 12
+18 3 10 13
+19 6 13 20
+20 7 21 22
+21 8 21 27
+22 8 21 28
+23 7 21 39
+24 9 22 23
+25 8 22 29
+26 8 22 30
+27 10 23 31
+28 10 23 32
+29 10 24 33
+30 10 24 34
+31 9 25 24
+32 7 25 26
+33 8 25 35
+34 8 25 36
+35 8 26 37
+36 8 26 38
+37 7 26 42
+38 8 39 40
+39 8 39 41
+40 7 39 42
+41 8 42 43
+42 8 42 44
+
+Angles
+
+1 1 5 1 4
+2 2 7 1 5
+3 2 8 1 5
+4 2 7 1 4
+5 2 8 1 4
+6 3 7 1 8
+7 4 7 2 3
+8 5 3 2 11
+9 6 7 2 11
+10 7 2 3 6
+11 8 1 7 2
+12 2 1 7 14
+13 2 1 7 15
+14 9 2 7 14
+15 9 2 7 15
+16 1 14 7 15
+17 3 1 8 9
+18 2 1 8 16
+19 2 1 8 17
+20 2 9 8 16
+21 2 9 8 17
+22 1 16 8 17
+23 8 8 9 10
+24 2 8 9 18
+25 2 8 9 19
+26 9 10 9 18
+27 9 10 9 19
+28 1 18 9 19
+29 6 9 10 12
+30 4 9 10 13
+31 5 13 10 12
+32 7 10 13 20
+33 10 22 21 27
+34 10 22 21 28
+35 11 22 21 39
+36 12 27 21 28
+37 10 39 21 27
+38 10 39 21 28
+39 13 21 22 23
+40 10 21 22 29
+41 10 21 22 30
+42 14 23 22 29
+43 14 23 22 30
+44 12 29 22 30
+45 15 22 23 31
+46 15 22 23 32
+47 16 31 23 32
+48 15 25 24 33
+49 15 25 24 34
+50 16 33 24 34
+51 13 26 25 24
+52 14 24 25 35
+53 14 24 25 36
+54 10 26 25 35
+55 10 26 25 36
+56 12 35 25 36
+57 10 25 26 37
+58 10 25 26 38
+59 11 25 26 42
+60 12 37 26 38
+61 10 42 26 37
+62 10 42 26 38
+63 10 21 39 40
+64 10 21 39 41
+65 11 21 39 42
+66 12 40 39 41
+67 10 42 39 40
+68 10 42 39 41
+69 11 26 42 39
+70 10 26 42 43
+71 10 26 42 44
+72 10 39 42 43
+73 10 39 42 44
+74 12 43 42 44
+
+Dihedrals
+
+1 2 5 1 7 14
+2 2 5 1 7 15
+3 2 4 1 7 14
+4 2 4 1 7 15
+5 3 8 1 7 2
+6 4 8 1 7 14
+7 4 8 1 7 15
+8 2 5 1 8 16
+9 2 5 1 8 17
+10 2 4 1 8 16
+11 2 4 1 8 17
+12 5 7 1 8 9
+13 4 7 1 8 16
+14 4 7 1 8 17
+15 6 7 2 3 6
+16 7 11 2 3 6
+17 1 2 7 1 5
+18 1 2 7 1 4
+19 8 1 7 2 3
+20 9 14 7 2 3
+21 9 15 7 2 3
+22 10 1 7 2 11
+23 11 14 7 2 11
+24 11 15 7 2 11
+25 4 9 8 1 5
+26 4 9 8 1 4
+27 3 1 8 9 10
+28 4 1 8 9 18
+29 4 1 8 9 19
+30 2 16 8 9 18
+31 2 16 8 9 19
+32 2 17 8 9 18
+33 2 17 8 9 19
+34 1 10 9 8 16
+35 1 10 9 8 17
+36 10 8 9 10 12
+37 8 8 9 10 13
+38 11 18 9 10 12
+39 9 18 9 10 13
+40 11 19 9 10 12
+41 9 19 9 10 13
+42 6 9 10 13 20
+43 7 12 10 13 20
+44 13 27 21 22 29
+45 13 27 21 22 30
+46 13 28 21 22 29
+47 13 28 21 22 30
+48 14 39 21 22 23
+49 15 39 21 22 29
+50 15 39 21 22 30
+51 15 22 21 39 40
+52 15 22 21 39 41
+53 16 22 21 39 42
+54 13 27 21 39 40
+55 13 27 21 39 41
+56 13 28 21 39 40
+57 13 28 21 39 41
+58 12 23 22 21 27
+59 12 23 22 21 28
+60 17 21 22 23 31
+61 17 21 22 23 32
+62 18 29 22 23 31
+63 18 29 22 23 32
+64 18 30 22 23 31
+65 18 30 22 23 32
+66 17 26 25 24 33
+67 18 35 25 24 33
+68 18 36 25 24 33
+69 17 26 25 24 34
+70 18 35 25 24 34
+71 18 36 25 24 34
+72 12 24 25 26 37
+73 12 24 25 26 38
+74 13 35 25 26 37
+75 13 35 25 26 38
+76 13 36 25 26 37
+77 13 36 25 26 38
+78 14 42 26 25 24
+79 15 42 26 25 35
+80 15 42 26 25 36
+81 16 25 26 42 39
+82 15 25 26 42 43
+83 15 25 26 42 44
+84 13 37 26 42 43
+85 13 37 26 42 44
+86 13 38 26 42 43
+87 13 38 26 42 44
+88 15 42 39 21 27
+89 15 42 39 21 28
+90 16 21 39 42 26
+91 15 21 39 42 43
+92 15 21 39 42 44
+93 13 40 39 42 43
+94 13 40 39 42 44
+95 13 41 39 42 43
+96 13 41 39 42 44
+97 15 39 42 26 37
+98 15 39 42 26 38
+99 15 26 42 39 40
+100 15 26 42 39 41
+
+Impropers
+
+1 2 7 1 4 5
+2 2 8 1 4 5
+3 3 7 1 8 5
+4 3 7 1 8 4
+5 1 7 2 3 11
+6 4 1 7 2 14
+7 4 1 7 2 15
+8 2 1 7 14 15
+9 5 2 7 14 15
+10 3 1 8 9 16
+11 3 1 8 9 17
+12 2 1 8 16 17
+13 2 9 8 16 17
+14 4 8 9 10 18
+15 4 8 9 10 19
+16 2 8 9 18 19
+17 5 10 9 18 19
+18 1 9 10 13 12
+19 7 22 21 27 28
+20 8 22 21 39 27
+21 8 22 21 39 28
+22 7 39 21 28 27
+23 9 21 22 23 29
+24 9 21 22 23 30
+25 7 21 22 29 30
+26 10 23 22 29 30
+27 6 22 23 31 32
+28 6 25 24 33 34
+29 9 26 25 24 35
+30 9 26 25 24 36
+31 10 24 25 35 36
+32 7 26 25 35 36
+33 7 25 26 37 38
+34 8 25 26 42 37
+35 8 25 26 42 38
+36 7 42 26 38 37
+37 7 21 39 40 41
+38 8 21 39 42 40
+39 8 21 39 42 41
+40 7 42 39 41 40
+41 8 26 42 39 43
+42 8 26 42 39 44
+43 7 26 42 43 44
+44 7 39 42 43 44
diff --git a/src/fix_bond_react.cpp b/src/fix_bond_react.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ea05cc025097d444b0da7d11e171321afecf0602
--- /dev/null
+++ b/src/fix_bond_react.cpp
@@ -0,0 +1,2696 @@
+/* ----------------------------------------------------------------------
+LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+http://lammps.sandia.gov, Sandia National Laboratories
+Steve Plimpton, sjplimp@sandia.gov
+
+Copyright (2003) Sandia Corporation. Under the terms of Contract
+DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+certain rights in this software. This software is distributed under
+the GNU General Public License.
+
+See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+Contributing Author: Jacob Gissinger (jacob.gissinger@colorado.edu)
+------------------------------------------------------------------------- */
+
+#include <math.h>
+#include <string.h>
+#include "math.h"
+#include "mpi.h"
+#include "string.h"
+#include "stdlib.h"
+#include "fix_bond_react.h"
+#include "update.h"
+#include "modify.h"
+#include "respa.h"
+#include "atom.h"
+#include "atom_vec.h"
+#include "force.h"
+#include "pair.h"
+#include "comm.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "random_mars.h"
+#include "memory.h"
+#include "citeme.h"
+#include "error.h"
+#include "molecule.h"
+#include "group.h"
+#include "fix_nve_limit.h"
+#include <algorithm>
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+#define BIG 1.0e20
+#define DELTA 16
+#define MAXLINE 256
+#define MAXGUESS 20
+
+// various statuses of superimpose algorithm:
+// ACCEPT: site successfully matched to pre-reacted template
+// REJECT: site does not match pre-reacted template
+// PROCEED: normal execution (non-guessing mode)
+// CONTINUE: a neighbor has been assigned, skip to next neighbor
+// GUESSFAIL: a guess has failed (if no more restore points, status = 'REJECT')
+// RESTORE: restore mode, load most recent restore point
+enum{ACCEPT,REJECT,PROCEED,CONTINUE,GUESSFAIL,RESTORE};
+
+/* ---------------------------------------------------------------------- */
+
+FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
+Fix(lmp, narg, arg)
+{
+ setbuf(stdout, NULL);
+
+ fix1 = NULL;
+ fix2 = NULL;
+ fix3 = NULL;
+ fix4 = NULL;
+
+ if (narg < 8) error->all(FLERR,"Illegal fix bond/react command 0.0");
+
+ MPI_Comm_rank(world,&me);
+ MPI_Comm_size(world,&nprocs);
+
+ attempted_rxn = 0;
+ ghostcheck_flag = 0;
+ force_reneighbor = 1;
+ next_reneighbor = -1;
+ vector_flag = 1;
+ global_freq = 1;
+ extvector = 0;
+ rxnID = 0;
+ status = PROCEED;
+
+ // these group names are reserved for use exclusively by bond/react
+ master_group = (char *) "bond_react_MASTER_group";
+
+ // let's find number of reactions specified
+ nreacts = 0;
+ for (int i = 3; i < narg; i++) {
+ if (strcmp(arg[i],"react") == 0) {
+ nreacts++;
+ i = i + 6; // skip past mandatory arguments
+ if (i > narg) error->all(FLERR,"Illegal fix bond/react command 0.1");
+ }
+ }
+
+ if (nreacts == 0) error->all(FLERR,"Illegal fix bond/react command: missing mandatory 'react' argument");
+
+ size_vector = nreacts;
+
+ int iarg = 3;
+ stabilization_flag = 0;
+ while (strcmp(arg[iarg],"react") != 0) {
+ if (strcmp(arg[iarg],"stabilization") == 0) {
+ if (strcmp(arg[iarg+1],"no") == 0) {
+ if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/react command 0.2");
+ iarg += 2;
+ }
+ if (strcmp(arg[iarg+1],"yes") == 0) {
+ if (iarg+4 > narg) error->all(FLERR,"Illegal fix bond/react command 0.21");
+ int n = strlen(arg[iarg+2]) + 1;
+ exclude_group = new char[n];
+ strcpy(exclude_group,arg[iarg+2]);
+ stabilization_flag = 1;
+ nve_limit_xmax = arg[iarg+3];
+ iarg += 4;
+ }
+ }
+ }
+
+ // set up common variables as vectors of length 'nreacts'
+ // nevery, cutoff, onemol, twomol, superimpose file
+
+ // 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(unreacted_mol,nreacts,"bond/react:unreacted_mol");
+ memory->create(reacted_mol,nreacts,"bond/react:reacted_mol");
+ memory->create(fraction,nreacts,"bond/react:fraction");
+ memory->create(seed,nreacts,"bond/react:seed");
+ memory->create(limit_duration,nreacts,"bond/react:limit_duration");
+ memory->create(stabilize_steps_flag,nreacts,"bond/react:stabilize_steps_flag");
+ memory->create(iatomtype,nreacts,"bond/react:iatomtype");
+ memory->create(jatomtype,nreacts,"bond/react:jatomtype");
+ memory->create(ibonding,nreacts,"bond/react:ibonding");
+ memory->create(jbonding,nreacts,"bond/react:jbonding");
+ 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");
+ memory->create(ghostly_rxn_count,nreacts,"bond/react:ghostly_rxn_count");
+ memory->create(reaction_count_total,nreacts,"bond/react:reaction_count_total");
+
+ for (int i = 0; i < nreacts; i++) {
+ fraction[i] = 1;
+ seed[i] = 12345;
+ stabilize_steps_flag[i] = 0;
+ // set default limit duration to 60 timesteps
+ limit_duration[i] = 60;
+ reaction_count[i] = 0;
+ local_rxn_count[i] = 0;
+ ghostly_rxn_count[i] = 0;
+ reaction_count_total[i] = 0;
+ }
+
+ char **files;
+ files = new char*[nreacts];
+
+ int rxn = 0;
+ while (iarg < narg && strcmp(arg[iarg],"react") == 0) {
+
+ iarg++;
+
+ iarg++; // read in reaction name here
+ //for example, rxn_name[rxn] = ...
+
+ int igroup = group->find(arg[iarg++]);
+ if (igroup == -1) error->all(FLERR,"Could not find fix group ID");
+ groupbits[rxn] = group->bitmask[igroup];
+
+ nevery[rxn] = force->inumeric(FLERR,arg[iarg++]);
+ if (nevery[rxn] <= 0) error->all(FLERR,"Illegal fix bond/react command 0.4");
+
+ 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;
+
+ unreacted_mol[rxn] = atom->find_molecule(arg[iarg++]);
+ if (unreacted_mol[rxn] == -1) error->all(FLERR,"Unreacted molecule template ID for "
+ "fix bond/react does not exist");
+ reacted_mol[rxn] = atom->find_molecule(arg[iarg++]);
+ if (reacted_mol[rxn] == -1) error->all(FLERR,"Reacted molecule template ID for "
+ "fix bond/react does not exist");
+
+ //read superimpose file
+ files[rxn] = new char[strlen(arg[iarg])+1];
+ strcpy(files[rxn],arg[iarg]);
+ iarg++;
+
+ while (iarg < narg && strcmp(arg[iarg],"react") != 0 ) {
+ if (strcmp(arg[iarg],"prob") == 0) {
+ if (iarg+3 > narg) error->all(FLERR,"Illegal fix bond/react command 0.6");
+ fraction[rxn] = force->numeric(FLERR,arg[iarg+1]);
+ seed[rxn] = force->inumeric(FLERR,arg[iarg+2]);
+ if (fraction[rxn] < 0.0 || fraction[rxn] > 1.0)
+ error->all(FLERR,"Illegal fix bond/react command");
+ if (seed[rxn] <= 0) error->all(FLERR,"Illegal fix bond/react command 0.7");
+ iarg += 3;
+ } else if (strcmp(arg[iarg],"stabilize_steps") == 0) {
+ if (stabilization_flag == 0) error->all(FLERR,"Stabilize_steps keyword used without stabilization keyword");
+ if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/react command 0.8");
+ limit_duration[rxn] = force->numeric(FLERR,arg[iarg+1]);
+ stabilize_steps_flag[rxn] = 1;
+ iarg += 2;
+ } else error->all(FLERR,"Illegal fix bond/react command 0.9");
+ }
+ rxn++;
+ }
+
+ max_natoms = 0; // the number of atoms in largest molecule template
+ for (int myrxn = 0; myrxn < nreacts; myrxn++) {
+ twomol = atom->molecules[reacted_mol[myrxn]];
+ max_natoms = MAX(max_natoms,twomol->natoms);
+ }
+
+ memory->create(equivalences,max_natoms,2,nreacts,"bond/react:equivalences");
+ memory->create(reverse_equiv,max_natoms,2,nreacts,"bond/react:reverse_equiv");
+ memory->create(edge,max_natoms,nreacts,"bond/react:edge");
+ memory->create(landlocked_atoms,max_natoms,nreacts,"bond/react:landlocked_atoms");
+
+ // read all superimpose files afterward
+ for (int i = 0; i < nreacts; i++) {
+ open(files[i]);
+ onemol = atom->molecules[unreacted_mol[i]];
+ twomol = atom->molecules[reacted_mol[i]];
+ read(i);
+ fclose(fp);
+ iatomtype[i] = onemol->type[ibonding[i]-1];
+ jatomtype[i] = onemol->type[jbonding[i]-1];
+ find_landlocked_atoms(i);
+ }
+
+ for (int i = 0; i < nreacts; i++) {
+ delete [] files[i];
+ }
+ delete [] files;
+
+ if (atom->molecular != 1)
+ error->all(FLERR,"Cannot use fix bond/react with non-molecular systems");
+
+ // initialize Marsaglia RNG with processor-unique seed
+
+ random = new class RanMars*[nreacts];
+ for (int i = 0; i < nreacts; i++) {
+ random[i] = new RanMars(lmp,seed[i] + me);
+ }
+
+ // set comm sizes needed by this fix
+ // forward is big due to comm of broken bonds and 1-2 neighbors
+
+ comm_forward = MAX(2,2+atom->maxspecial);
+ comm_reverse = 2;
+
+ // allocate arrays local to this fix
+
+ nmax = 0;
+ partner = finalpartner = NULL;
+ distsq = NULL;
+ maxcreate = 0;
+ created = NULL;
+ local_ncreate = NULL;
+ ncreate = NULL;
+ allncreate = 0;
+ local_num_mega = 0;
+ ghostly_num_mega = 0;
+ restore = NULL;
+
+ // zero out stats
+ global_megasize = 0;
+ avail_guesses = 0;
+ glove_counter = 0;
+ guess_branch = new int[MAXGUESS]();
+ pioneer_count = new int[max_natoms];
+ local_mega_glove = NULL;
+ ghostly_mega_glove = NULL;
+ global_mega_glove = NULL;
+
+ // these are merely loop indices that became important
+ pion = neigh = trace = 0;
+
+ id_fix1 = NULL;
+ id_fix2 = NULL;
+ id_fix3 = NULL;
+ id_fix4 = NULL;
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixBondReact::~FixBondReact()
+{
+
+ // unregister callbacks to this fix from Atom class
+ atom->delete_callback(id,0);
+
+ for (int i = 0; i < nreacts; i++) {
+ delete random[i];
+ }
+ delete [] random;
+
+ memory->destroy(partner);
+ memory->destroy(finalpartner);
+ memory->destroy(local_ncreate);
+ memory->destroy(ncreate);
+ memory->destroy(distsq);
+ memory->destroy(created);
+ memory->destroy(edge);
+ memory->destroy(equivalences);
+ memory->destroy(reverse_equiv);
+
+ memory->destroy(nevery);
+ memory->destroy(cutsq);
+ memory->destroy(unreacted_mol);
+ memory->destroy(reacted_mol);
+ memory->destroy(fraction);
+ memory->destroy(seed);
+ memory->destroy(limit_duration);
+ memory->destroy(stabilize_steps_flag);
+
+ memory->destroy(iatomtype);
+ memory->destroy(jatomtype);
+ memory->destroy(ibonding);
+ memory->destroy(jbonding);
+ memory->destroy(groupbits);
+ memory->destroy(reaction_count);
+ memory->destroy(local_rxn_count);
+ memory->destroy(ghostly_rxn_count);
+ memory->destroy(reaction_count_total);
+
+ if (attempted_rxn == 1) {
+ memory->destroy(restore_pt);
+ memory->destroy(restore);
+ memory->destroy(glove);
+ memory->destroy(pioneers);
+ memory->destroy(landlocked_atoms);
+ memory->destroy(local_mega_glove);
+ memory->destroy(ghostly_mega_glove);
+ }
+
+ if (ghostcheck_flag == 1) memory->destroy(global_mega_glove);
+
+ if (stabilization_flag == 1) {
+ delete [] exclude_group;
+
+ // 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 (id_fix2 == NULL && modify->nfix) modify->delete_fix(id_fix2);
+ delete [] id_fix2;
+
+ if (id_fix3 == NULL && modify->nfix) modify->delete_fix(id_fix3);
+ delete [] id_fix3;
+
+ if (id_fix4 == NULL && modify->nfix) modify->delete_fix(id_fix4);
+ delete [] id_fix4;
+
+ delete [] guess_branch;
+ delete [] pioneer_count;
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixBondReact::setmask()
+{
+ int mask = 0;
+ mask |= POST_INTEGRATE;
+ mask |= POST_INTEGRATE_RESPA;
+ return mask;
+}
+
+/* ----------------------------------------------------------------------
+let's add an internal nve/limit fix for relaxation of reaction sites
+also let's add our per-atom property fix here!
+this per-atom property will state the timestep an atom was 'limited'
+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
+ //these are recently reacted atoms being relaxed
+ int len = strlen("limit_tags") + 1;
+ id_fix2 = new char[len];
+ strcpy(id_fix2,"limit_tags"); //note: this is purposefully same as property 'name'
+
+ int ifix = modify->find_fix(id_fix2);
+ if (ifix == -1) {
+ char **newarg = new char*[6];
+ newarg[0] = (char *) "limit_tags";
+ newarg[1] = (char *) "all"; // group ID is ignored
+ newarg[2] = (char *) "property/atom";
+ newarg[3] = (char *) "i_limit_tags";
+ newarg[4] = (char *) "ghost";
+ newarg[5] = (char *) "yes";
+ modify->add_fix(6,newarg);
+ fix2 = (FixPropertyAtom *) modify->fix[modify->nfix-1];
+ delete [] newarg;
+ }
+
+ // per-atom properties already initialized to zero (not in group)
+ // let's do it anyway for clarity
+ int flag;
+ int index = atom->find_custom(id_fix2,flag);
+ int *i_limit_tags = atom->ivector[index];
+
+ for (int i = 0; i < atom->nlocal; i++)
+ i_limit_tags[i] = 0;
+
+ // create master_group if not already existing
+ if (group->find(master_group) == -1) {
+ 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;
+ }
+
+ //let's add the statted_tags per-atom property fix
+ //this are atoms subjec to the system-wide thermostat
+ len = strlen("statted_tags") + 1;
+ id_fix3 = new char[len];
+ strcpy(id_fix3,"statted_tags"); //note: this is purposefully same as property 'name'
+
+ ifix = modify->find_fix(id_fix3);
+ if (ifix == -1) {
+ char **newarg = new char*[6];
+ newarg[0] = (char *) "statted_tags";
+ 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);
+ fix3 = (FixPropertyAtom *) modify->fix[modify->nfix-1];
+ delete [] newarg;
+ }
+
+ //intialize per-atom statted_flags to 1
+ index = atom->find_custom(id_fix3,flag);
+ int *i_statted_tags = atom->ivector[index];
+
+ for (int i = 0; i < atom->nlocal; i++)
+ i_statted_tags[i] = 1;
+
+ 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";
+ newarg[3] = (char *) "property";
+ newarg[4] = (char *) "statted_tags";
+ group->assign(5,newarg);
+ delete [] newarg;
+ }
+
+ // let's create a new nve/limit fix to limit newly reacted atoms
+ len = strlen("bond_react_MASTER_nve_limit") + 1;
+ id_fix1 = new char[len];
+ strcpy(id_fix1,"bond_react_MASTER_nve_limit");
+
+ ifix = modify->find_fix(id_fix1);
+
+ if (ifix == -1) {
+ char **newarg = new char*[4];
+ newarg[0] = id_fix1;
+ newarg[1] = master_group;
+ newarg[2] = (char *) "nve/limit";
+ newarg[3] = nve_limit_xmax;
+ modify->add_fix(4,newarg);
+ fix1 = (FixNVELimit *) modify->fix[modify->nfix-1];
+ delete [] newarg;
+ }
+
+ }
+
+ //let's add the react_tags per-atom property fix
+ //this per-atom property is the ID of the 'react' argument which recently caused atom to react
+ //so that atoms which wander between processors may be released to global thermostat at the proper time
+ len = strlen("react_tags") + 1;
+ id_fix4 = new char[len];
+ strcpy(id_fix4,"react_tags"); //note: this is purposefully same as property 'name'
+
+ ifix = modify->find_fix(id_fix4);
+ if (ifix == -1) { //ifix
+ char **newarg = new char*[6];
+ newarg[0] = (char *) "react_tags";
+ newarg[1] = (char *) "all"; // group ID is ignored
+ newarg[2] = (char *) "property/atom";
+ newarg[3] = (char *) "i_react_tags";
+ newarg[4] = (char *) "ghost";
+ newarg[5] = (char *) "yes";
+ modify->add_fix(6,newarg);
+ fix4 = (FixPropertyAtom *) modify->fix[modify->nfix-1];
+ delete [] newarg;
+ }
+
+ //per-atom values initalized to 0
+ index = atom->find_custom(id_fix4,flag);
+ int *i_react_tags = atom->ivector[index];
+
+ for (int i = 0; i < atom->nlocal; i++)
+ i_react_tags[i] = 0;
+
+ // 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;
+ }
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixBondReact::init()
+{
+
+ // warn if more than one bond/react fix
+
+ int count = 0;
+ for (int i = 0; i < modify->nfix; i++)
+ if (strcmp(modify->fix[i]->style,"bond/react") == 0) count++;
+ if (count > 1 && comm->me == 0) error->warning(FLERR,"More than one fix bond/react");
+
+ if (strstr(update->integrate_style,"respa"))
+ nlevels_respa = ((Respa *) update->integrate)->nlevels;
+
+ // 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]])
+ error->all(FLERR,"Fix bond/react cutoff is longer than pairwise cutoff");
+ }
+
+ // need a half neighbor list, built every Nevery steps
+ int irequest = neighbor->request(this,instance_me);
+ neighbor->requests[irequest]->pair = 0;
+ neighbor->requests[irequest]->fix = 1;
+ neighbor->requests[irequest]->occasional = 1;
+
+ lastcheck = -1;
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixBondReact::init_list(int id, NeighList *ptr)
+{
+ list = ptr;
+}
+
+/* ----------------------------------------------------------------------
+ Identify all pairs of potentially reactive atoms for this time step.
+ This function is modified from LAMMPS’ fix bond/create.
+---------------------------------------------------------------------- */
+
+void FixBondReact::post_integrate()
+{
+
+ int inum,jnum,itype,jtype,n1,n2,n3,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++) {
+ if (!(update->ntimestep % nevery[i])) {
+ nevery_check = 0;
+ break;
+ }
+ }
+
+ for (int i = 0; i < nreacts; i++) {
+ reaction_count[i] = 0;
+ local_rxn_count[i] = 0;
+ ghostly_rxn_count[i] = 0;
+ }
+
+ if (nevery_check) {
+ unlimit_bond();
+ return;
+ }
+
+ // acquire updated ghost atom positions
+ // necessary b/c are calling this after integrate, but before Verlet comm
+
+ 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) {
+ memory->destroy(partner);
+ memory->destroy(finalpartner);
+ memory->destroy(distsq);
+ memory->destroy(local_ncreate);
+ memory->destroy(ncreate);
+ 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(local_ncreate,nreacts,"bond/react:local_ncreate");
+ memory->create(ncreate,nreacts,"bond/react:ncreate");
+ probability = distsq;
+ }
+
+ // reset create counts
+ for (int i = 0; i < nreacts; i++) {
+ local_ncreate[i] = 0;
+ ncreate[i] = 0;
+ }
+
+ 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;
+ tagint **bond_atom = atom->bond_atom;
+ int *num_bond = atom->num_bond;
+ 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(id_fix2,flag);
+ int *i_limit_tags = atom->ivector[index1];
+
+ int i,j;
+
+ for (int myrxn = 0; myrxn < nreacts; myrxn++) {
+
+ for (int ii = 0; ii < nall; ii++) {
+ partner[ii] = 0;
+ finalpartner[ii] = 0;
+ distsq[ii] = 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;
+ }
+ }
+ }
+
+ // reverse comm of distsq and partner
+ // not needed if newton_pair off since I,J pair was seen by both procs
+
+ commflag = 2;
+ if (force->newton_pair) comm->reverse_comm_fix(this);
+
+ // each atom now knows its winning partner
+ // 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) {
+ for (int i = 0; i < nlocal; i++)
+ if (partner[i]) probability[i] = random[myrxn]->uniform();
+ }
+
+ commflag = 2;
+ comm->forward_comm_fix(this,2);
+
+ // consider for reaction:
+ // only if both atoms list each other as winning bond partner
+ // and probability constraint is satisfied
+ // if other atom is owned by another proc, it should do same thing
+
+ int **bond_type = atom->bond_type;
+ int newton_bond = force->newton_bond;
+
+ int temp_ncreate = 0;
+ for (int i = 0; i < nlocal; i++) {
+ if (partner[i] == 0) {
+ continue;
+ }
+ j = atom->map(partner[i]);
+ if (partner[j] != tag[i]) {
+ continue;
+ }
+
+ // apply probability constraint using RN for atom with smallest ID
+
+ if (fraction[myrxn] < 1.0) {
+ if (tag[i] < tag[j]) {
+ if (probability[i] >= fraction[myrxn]) continue;
+ } else {
+ if (probability[j] >= fraction[myrxn]) continue;
+ }
+ }
+
+ // store final created bond partners and count the rxn possibility once
+
+ finalpartner[i] = tag[j];
+ finalpartner[j] = tag[i];
+
+ if (tag[i] < tag[j]) temp_ncreate++;
+ }
+
+ local_ncreate[myrxn] = 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;
+ }
+
+ // 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
+ // this is because we take care of updating topology later (and differently)
+ for (int i = 0; i < nlocal; i++) {
+
+ if (finalpartner[i] == 0) continue;
+
+ 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) {
+ 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];
+ } else {
+ created[ncreate[myrxn]][0][myrxn] = finalpartner[i];
+ created[ncreate[myrxn]][1][myrxn] = tag[i];
+ }
+ ncreate[myrxn]++;
+ }
+ }
+ unlimit_bond(); //free atoms that have been relaxed
+ }
+
+ // break loop if no even eligible bonding atoms were found (on any proc)
+ int some_chance;
+
+ allncreate = 0;
+ for (int i = 0; i < nreacts; i++)
+ allncreate += ncreate[i];
+
+ MPI_Allreduce(&allncreate,&some_chance,1,MPI_INT,MPI_SUM,world);
+ if (!some_chance) {
+ unlimit_bond();
+ return;
+ }
+
+ // run through the superimpose algorithm
+ // this checks if simulation topology matches unreacted mol template
+ superimpose_algorithm();
+ // free atoms that have been limited after reacting
+ unlimit_bond();
+
+}
+
+/* ----------------------------------------------------------------------
+ Set up global variables. Loop through all pairs; loop through Pioneers
+ until Superimpose Algorithm is completed for each pair.
+------------------------------------------------------------------------- */
+
+void FixBondReact::superimpose_algorithm()
+{
+ local_num_mega = 0;
+ ghostly_num_mega = 0;
+
+ // quick description of important global indices you'll see floating about:
+ // 'pion' is the pioneer loop index
+ // 'neigh' in the first neighbor index
+ // 'trace' retraces the first nieghbors
+ // trace: once you choose a first neighbor, you then check for other nieghbors of same type
+
+ int **nspecial = atom->nspecial;
+ tagint **special = atom->special;
+ int *type = atom->type;
+
+ if (attempted_rxn == 1) {
+ memory->destroy(restore_pt);
+ memory->destroy(restore);
+ memory->destroy(glove);
+ memory->destroy(pioneers);
+ memory->destroy(local_mega_glove);
+ memory->destroy(ghostly_mega_glove);
+ }
+
+ memory->create(glove,max_natoms,2,"bond/react:glove");
+ memory->create(restore_pt,MAXGUESS,4,"bond/react:restore_pt");
+ memory->create(pioneers,max_natoms,"bond/react:pioneers");
+ memory->create(restore,max_natoms,MAXGUESS,"bond/react:restore");
+ memory->create(local_mega_glove,max_natoms+1,allncreate,"bond/react:local_mega_glove");
+ memory->create(ghostly_mega_glove,max_natoms+1,allncreate,"bond/react:ghostly_mega_glove");
+
+ attempted_rxn = 1;
+
+ for (int i = 0; i < max_natoms+1; i++) {
+ for (int j = 0; j < allncreate; j++) {
+ local_mega_glove[i][j] = 0;
+ ghostly_mega_glove[i][j] = 0;
+ }
+ }
+
+ //let's finally begin the superimpose loop
+ for (rxnID = 0; rxnID < nreacts; rxnID++) {
+ for (lcl_inst = 0; lcl_inst < ncreate[rxnID]; lcl_inst++) {
+
+ onemol = atom->molecules[unreacted_mol[rxnID]];
+ twomol = atom->molecules[reacted_mol[rxnID]];
+
+ status = PROCEED;
+
+ glove_counter = 0;
+ for (int i = 0; i < max_natoms; i++) {
+ for (int j = 0; j < 2; j++) {
+ glove[i][j] = 0;
+ }
+ }
+
+ for (int i = 0; i < MAXGUESS; i++) {
+ guess_branch[i] = 0;
+ }
+
+ int myibonding = ibonding[rxnID];
+ int myjbonding = jbonding[rxnID];
+
+ glove[myibonding-1][0] = myibonding;
+ glove[myibonding-1][1] = created[lcl_inst][0][rxnID];
+ glove_counter++;
+ glove[myjbonding-1][0] = myjbonding;
+ glove[myjbonding-1][1] = created[lcl_inst][1][rxnID];
+ glove_counter++;
+
+ avail_guesses = 0;
+
+ for (int i = 0; i < max_natoms; i++)
+ pioneer_count[i] = 0;
+
+ for (int i = 0; i < onemol->nspecial[myibonding-1][0]; i++)
+ pioneer_count[onemol->special[myibonding-1][i]-1]++;
+
+ for (int i = 0; i < onemol->nspecial[myjbonding-1][0]; i++)
+ pioneer_count[onemol->special[myjbonding-1][i]-1]++;
+
+
+ int hang_catch = 0;
+ while (!(status == ACCEPT || status == REJECT)) {
+
+ for (int i = 0; i < max_natoms; i++) {
+ pioneers[i] = 0;
+ }
+
+ for (int i = 0; i < onemol->natoms; i++) {
+ if (glove[i][0] !=0 && pioneer_count[i] < onemol->nspecial[i][0] && edge[i][rxnID] == 0) {
+ pioneers[i] = 1;
+ }
+ }
+
+ // run through the pioneers
+ // due to use of restore points, 'pion' index can change in loop
+ for (pion = 0; pion < onemol->natoms; pion++) {
+ if (pioneers[pion] || status == GUESSFAIL) {
+ make_a_guess();
+ if (status == ACCEPT || status == REJECT) break;
+ }
+ }
+
+ if (status == ACCEPT) { // reaction site found successfully!
+ glove_ghostcheck();
+ }
+ hang_catch++;
+ // let's go ahead and catch the simplest of hangs
+ //if (hang_catch > onemol->natoms*4)
+ if (hang_catch > atom->nlocal*3) {
+ error->all(FLERR,"Excessive iteration of superimpose algorithm");
+ }
+ }
+ }
+ }
+
+ global_megasize = 0;
+
+ ghost_glovecast(); // consolidate all mega_gloves to all processors
+ dedup_mega_gloves(0); // make sure atoms aren't added to more than one reaction
+
+ MPI_Allreduce(&local_rxn_count[0],&reaction_count[0],nreacts,MPI_INT,MPI_SUM,world);
+
+ for (int i = 0; i < nreacts; i++)
+ reaction_count_total[i] += reaction_count[i];
+
+ // this assumes compute_vector is called from process 0
+ // ...so doesn't bother to bcast ghostly_rxn_count
+ if (me == 0)
+ for (int i = 0; i < nreacts; i++)
+ reaction_count_total[i] += ghostly_rxn_count[i];
+
+ // this updates topology next step
+ next_reneighbor = update->ntimestep;
+
+ // call limit_bond in 'global_mega_glove mode.' oh, and local mode
+ limit_bond(0); // add reacting atoms to nve/limit
+ limit_bond(1);
+ update_everything(); // change topology
+
+}
+
+/* ----------------------------------------------------------------------
+ Screen for obvious algorithm fails. This is the return point when a guess
+ has failed: check for available restore points.
+------------------------------------------------------------------------- */
+
+void FixBondReact::make_a_guess()
+{
+
+ int **nspecial = atom->nspecial;
+ tagint **special = atom->special;
+ int *type = atom->type;
+ int *mask = atom->mask;
+ int nfirst_neighs = onemol->nspecial[pion][0];
+
+ // per-atom property indicating if in bond/react master group
+ int flag;
+ int index1 = atom->find_custom(id_fix2,flag);
+ int *i_limit_tags = atom->ivector[index1];
+
+ if (status == GUESSFAIL && avail_guesses == 0) {
+ status = REJECT;
+ return;
+ }
+
+ if (status == GUESSFAIL && avail_guesses > 0) {
+ // load restore point
+ for (int i = 0; i < onemol->natoms; i++) {
+ glove[i][0] = restore[i][(avail_guesses*4)-4];
+ glove[i][1] = restore[i][(avail_guesses*4)-3];
+ pioneer_count[i] = restore[i][(avail_guesses*4)-2];
+ pioneers[i] = restore[i][(avail_guesses*4)-1];
+ }
+ pion = restore_pt[avail_guesses-1][0];
+ neigh = restore_pt[avail_guesses-1][1];
+ trace = restore_pt[avail_guesses-1][2];
+ glove_counter = restore_pt[avail_guesses-1][3];
+ status = RESTORE;
+ neighbor_loop();
+ if (status != PROCEED) return;
+ }
+
+ nfirst_neighs = onemol->nspecial[pion][0];
+
+ // check if any of first neighbors are in bond_react_MASTER_group
+ // if so, this constitutes a fail
+ // because still undergoing a previous reaction!
+ // could technically fail unnecessarily during a wrong guess if near edge atoms
+ // we accept this temporary and infrequent decrease in reaction occurences
+
+ for (int i = 0; i < nspecial[atom->map(glove[pion][1])][0]; i++) {
+ if (atom->map(special[atom->map(glove[pion][1])][i]) < 0) {
+ error->all(FLERR,"Fix bond/react needs ghost atoms from further away1"); // parallel issues.
+ }
+ if (i_limit_tags[(int)atom->map(special[atom->map(glove[pion][1])][i])] != 0) {
+ status = GUESSFAIL;
+ return;
+ }
+ }
+
+ // check for same number of neighbors between unreacted mol and simulation
+ if (nfirst_neighs != nspecial[atom->map(glove[pion][1])][0]) {
+ status = GUESSFAIL;
+ return;
+ }
+
+ // make sure all neighbors aren't already assigned
+ // an issue discovered for coarse-grained example
+ int assigned_count = 0;
+ for (int i = 0; i < nfirst_neighs; i++)
+ for (int j = 0; j < onemol->natoms; j++)
+ if (special[atom->map(glove[pion][1])][i] == glove[j][1]) {
+ assigned_count++;
+ break;
+ }
+
+ 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];
+
+ for (int i = 0; i < atom->ntypes; i++) {
+ mol_ntypes[i] = 0;
+ lcl_ntypes[i] = 0;
+ }
+
+ for (int i = 0; i < nfirst_neighs; i++) {
+ mol_ntypes[(int)onemol->type[(int)onemol->special[pion][i]-1]-1]++;
+ lcl_ntypes[(int)type[(int)atom->map(special[atom->map(glove[pion][1])][i])]-1]++; //added -1
+ }
+
+ for (int i = 0; i < atom->ntypes; i++) {
+ if (mol_ntypes[i] != lcl_ntypes[i]) {
+ status = GUESSFAIL;
+ return;
+ }
+ }
+
+ // okay everything seems to be in order. let's assign some ID pairs!!!
+ neighbor_loop();
+
+}
+
+/* ----------------------------------------------------------------------
+ Loop through all First Bonded Neighbors of the current Pioneer.
+ Prepare appropriately if we are in Restore Mode.
+------------------------------------------------------------------------- */
+
+void FixBondReact::neighbor_loop()
+{
+
+ int **nspecial = atom->nspecial;
+ tagint **special = atom->special;
+ int *type = atom->type;
+
+ int nfirst_neighs = onemol->nspecial[pion][0];
+
+ if (status == RESTORE) {
+ check_a_neighbor();
+ return;
+ } else neigh = 0;
+
+ for (neigh; neigh < nfirst_neighs; neigh++) {
+ if (glove[(int)onemol->special[pion][neigh]-1][0] == 0) {
+ check_a_neighbor();
+ }
+ }
+
+ // status should still = PROCEED
+
+}
+
+/* ----------------------------------------------------------------------
+ Check if we can assign this First Neighbor to pre-reacted template
+ without guessing. If so, do it! If not, call crosscheck_the_nieghbor().
+------------------------------------------------------------------------- */
+
+void FixBondReact::check_a_neighbor()
+{
+
+ int **nspecial = atom->nspecial;
+ tagint **special = atom->special;
+ int *type = atom->type;
+ int nfirst_neighs = onemol->nspecial[pion][0];
+
+ if (status != RESTORE) {
+ // special consideration for hydrogen atoms (and all first neighbors bonded to no other atoms) (and aren't edge atoms)
+ if (onemol->nspecial[(int)onemol->special[pion][neigh]-1][0] == 1 && edge[(int)onemol->special[pion][neigh]-1][rxnID] == 0) {
+
+ for (int i = 0; i < nfirst_neighs; i++) {
+
+ if (type[(int)atom->map(special[(int)atom->map(glove[pion][1])][i])] == onemol->type[(int)onemol->special[pion][neigh]-1] &&
+ nspecial[(int)atom->map(special[(int)atom->map(glove[pion][1])][i])][0] == 1) {
+
+ int already_assigned = 0;
+ for (int j = 0; j < onemol->natoms; j++) {
+ if (glove[j][1] == special[atom->map(glove[pion][1])][i]) {
+ already_assigned = 1;
+ break;
+ }
+ }
+
+ if (already_assigned == 0) {
+ glove[(int)onemol->special[pion][neigh]-1][0] = onemol->special[pion][neigh];
+ glove[(int)onemol->special[pion][neigh]-1][1] = special[(int)atom->map(glove[pion][1])][i];
+
+ //another check for ghost atoms. perhaps remove the one in make_a_guess
+ if (atom->map(glove[(int)onemol->special[pion][neigh]-1][1]) < 0) {
+ error->all(FLERR,"Fix bond/react needs ghost atoms from further away2");
+ }
+
+ for (int j = 0; j < onemol->nspecial[onemol->special[pion][neigh]-1][0]; j++) {
+ pioneer_count[onemol->special[onemol->special[pion][neigh]-1][j]-1]++;
+ }
+
+ glove_counter++;
+ if (glove_counter == onemol->natoms) {
+ status = ACCEPT;
+ ring_check();
+ return;
+ }
+ // status should still == PROCEED
+ return;
+ }
+ }
+ }
+ // we are here if no matching atom found
+ status = GUESSFAIL;
+ return;
+ }
+ }
+
+ crosscheck_the_neighbor();
+ if (status != PROCEED) {
+ if (status == CONTINUE)
+ status = PROCEED;
+ return;
+ }
+
+ // finally ready to match non-duplicate, non-edge atom IDs!!
+
+ for (int i = 0; i < nfirst_neighs; i++) {
+
+ if (type[atom->map((int)special[(int)atom->map(glove[pion][1])][i])] == onemol->type[(int)onemol->special[pion][neigh]-1]) {
+ int already_assigned = 0;
+
+ //check if a first neighbor of the pioneer is already assigned to pre-reacted template
+ for (int j = 0; j < onemol->natoms; j++) {
+ if (glove[j][1] == special[atom->map(glove[pion][1])][i]) {
+ already_assigned = 1;
+ break;
+ }
+ }
+
+ if (already_assigned == 0) {
+ glove[(int)onemol->special[pion][neigh]-1][0] = onemol->special[pion][neigh];
+ glove[(int)onemol->special[pion][neigh]-1][1] = special[(int)atom->map(glove[pion][1])][i];
+
+ //another check for ghost atoms. perhaps remove the one in make_a_guess
+ if (atom->map(glove[(int)onemol->special[pion][neigh]-1][1]) < 0) {
+ error->all(FLERR,"Fix bond/react needs ghost atoms from further away3");
+ }
+
+ for (int ii = 0; ii < onemol->nspecial[onemol->special[pion][neigh]-1][0]; ii++) {
+ pioneer_count[onemol->special[onemol->special[pion][neigh]-1][ii]-1]++;
+ }
+
+ glove_counter++;
+ if (glove_counter == onemol->natoms) {
+ status = ACCEPT;
+ ring_check();
+ return;
+ // will never complete here when there are edge atoms
+ // ...actually that could be wrong if people get creative...shouldn't affect anything
+ }
+ // status should still = PROCEED
+ return;
+ }
+ }
+ }
+
+ // status is still 'PROCEED' if we are here!
+
+}
+
+/* ----------------------------------------------------------------------
+ Check if there a viable guess to be made. If so, prepare to make a
+ guess by recording a restore point.
+------------------------------------------------------------------------- */
+
+void FixBondReact::crosscheck_the_neighbor()
+{
+ int **nspecial = atom->nspecial;
+ tagint **special = atom->special;
+ int *type = atom->type;
+
+ int nfirst_neighs = onemol->nspecial[pion][0];
+
+ if (status == RESTORE) {
+ inner_crosscheck_loop();
+ return;
+ } else trace = 0;
+
+ for (trace; trace < nfirst_neighs; trace++) {
+ if (neigh!=trace && onemol->type[(int)onemol->special[pion][neigh]-1] == onemol->type[(int)onemol->special[pion][trace]-1] &&
+ glove[onemol->special[pion][trace]-1][0] == 0) {
+
+ if (avail_guesses == MAXGUESS) {
+ error->warning(FLERR,"Fix bond/react failed because MAXGUESS set too small. ask developer for info");
+ status = GUESSFAIL;
+ return;
+ }
+ avail_guesses++;
+ for (int i = 0; i < onemol->natoms; i++) {
+ restore[i][(avail_guesses*4)-4] = glove[i][0];
+ restore[i][(avail_guesses*4)-3] = glove[i][1];
+ restore[i][(avail_guesses*4)-2] = pioneer_count[i];
+ restore[i][(avail_guesses*4)-1] = pioneers[i];
+ restore_pt[avail_guesses-1][0] = pion;
+ restore_pt[avail_guesses-1][1] = neigh;
+ restore_pt[avail_guesses-1][2] = trace;
+ restore_pt[avail_guesses-1][3] = glove_counter;
+ }
+
+ inner_crosscheck_loop();
+ return;
+ }
+ }
+
+ // status is still 'PROCEED' if we are here!
+
+}
+
+/* ----------------------------------------------------------------------
+ We are ready to make a guess. If there are multiple possible choices
+ for this guess, keep track of these.
+------------------------------------------------------------------------- */
+
+void FixBondReact::inner_crosscheck_loop()
+{
+ int **nspecial = atom->nspecial;
+ tagint **special = atom->special;
+ int *type = atom->type;
+ // arbitrarily limited to 5 identical first neighbors
+ tagint tag_choices[5];
+ int nfirst_neighs = onemol->nspecial[pion][0];
+
+ int num_choices = 0;
+ for (int i = 0; i < nfirst_neighs; i++) {
+
+ int already_assigned = 0;
+ for (int j = 0; j < onemol->natoms; j++) {
+ if (glove[j][1] == special[atom->map(glove[pion][1])][i]) {
+ already_assigned = 1;
+ break;
+ }
+ }
+
+ if (already_assigned == 0 &&
+ type[(int)atom->map(special[atom->map(glove[pion][1])][i])] == onemol->type[(int)onemol->special[pion][neigh]-1]) {
+ if (num_choices > 5) { // here failed because too many identical first neighbors. but really no limit if situation arises
+ status = GUESSFAIL;
+ return;
+ }
+ tag_choices[num_choices++] = special[atom->map(glove[pion][1])][i];
+ }
+ }
+
+ // guess branch is for when multiple identical neighbors. then we guess each one in turn
+ // guess_branch must work even when avail_guesses = 0. so index accordingly!
+ // ...actually, avail_guesses should never be zero here anyway
+ if (guess_branch[avail_guesses-1] == 0) guess_branch[avail_guesses-1] = num_choices;
+
+ //std::size_t size = sizeof(tag_choices) / sizeof(tag_choices[0]);
+ std::sort(tag_choices, tag_choices + num_choices); //, std::greater<int>());
+ glove[onemol->special[pion][neigh]-1][0] = onemol->special[pion][neigh];
+ glove[onemol->special[pion][neigh]-1][1] = tag_choices[guess_branch[avail_guesses-1]-1];
+ guess_branch[avail_guesses-1]--;
+
+ //another check for ghost atoms. perhaps remove the one in make_a_guess
+ if (atom->map(glove[(int)onemol->special[pion][neigh]-1][1]) < 0) {
+ error->all(FLERR,"Fix bond/react needs ghost atoms from further away4");
+ }
+
+ if (guess_branch[avail_guesses-1] == 0) avail_guesses--;
+
+ for (int i = 0; i < onemol->nspecial[onemol->special[pion][neigh]-1][0]; i++) {
+ pioneer_count[onemol->special[onemol->special[pion][neigh]-1][i]-1]++;
+ }
+ glove_counter++;
+ if (glove_counter == onemol->natoms) {
+ status = ACCEPT;
+ ring_check();
+ return;
+ }
+ status = CONTINUE;
+
+}
+
+/* ----------------------------------------------------------------------
+ Check that newly assigned atoms have correct bonds
+ Necessary for certain ringed structures
+------------------------------------------------------------------------- */
+
+void FixBondReact::ring_check()
+{
+ // ring_check can be made more efficient by re-introducing 'frozen' atoms
+ // 'frozen' atoms have been assigned and also are no longer pioneers
+ int **nspecial = atom->nspecial;
+ tagint **special = atom->special;
+
+ for (int i = 0; i < onemol->natoms; i++) {
+ for (int j = 0; j < onemol->nspecial[i][0]; j++) {
+ int ring_fail = 1;
+ int ispecial = onemol->special[i][j];
+ for (int k = 0; k < nspecial[atom->map(glove[i][1])][0]; k++) {
+ if (special[atom->map(glove[i][1])][k] == glove[ispecial-1][1]) {
+ ring_fail = 0;
+ break;
+ }
+ }
+ if (ring_fail == 1) {
+ status = GUESSFAIL;
+ return;
+ }
+ }
+ }
+
+}
+
+/* ----------------------------------------------------------------------
+ Determine which pre-reacted template atoms are at least three bonds
+ away from edge atoms.
+------------------------------------------------------------------------- */
+
+void FixBondReact::find_landlocked_atoms(int myrxn)
+{
+ // landlocked_atoms are atoms for which all topology is contained in reacted template
+ // if dihedrals/impropers exist: this means that edge atoms are not in their 1-3 neighbor list
+ // note: due to various usage/defintions of impropers, treated same as dihedrals
+ // if angles exist: this means edge atoms not in their 1-2 neighbors list
+ // if just bonds: this just means that edge atoms are not landlocked
+ // Note: landlocked defined in terms of reacted template
+ // 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
+ 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;
+ }
+ int nspecial_limit = -1;
+ if (force->angle && twomol->angleflag) nspecial_limit = 0;
+
+ 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++) {
+ for (int k = 0; k < onemol->natoms; k++) {
+ if (equivalences[twomol->special[i][j]-1][1][myrxn] == k+1 && edge[k][myrxn] == 1) {
+ landlocked_atoms[i][myrxn] = 0;
+ }
+ }
+ }
+ }
+ }
+
+ // 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");
+ }
+
+}
+
+/* ----------------------------------------------------------------------
+let's dedup global_mega_glove
+allows for same site undergoing different pathways, in parallel
+------------------------------------------------------------------------- */
+
+void FixBondReact::dedup_mega_gloves(int dedup_mode)
+{
+ // dedup_mode == 0 for local_dedup
+ // dedup_mode == 1 for global_mega_glove
+ for (int i = 0; i < nreacts; i++) {
+ if (dedup_mode == 0) local_rxn_count[i] = 0;
+ if (dedup_mode == 1) ghostly_rxn_count[i] = 0;
+ }
+
+ int dedup_size;
+ if (dedup_mode == 0) {
+ dedup_size = local_num_mega;
+ } else if (dedup_mode == 1) {
+ dedup_size = global_megasize;
+ }
+
+ tagint dedup_glove[max_natoms+1][dedup_size];
+
+ if (dedup_mode == 0) {
+ for (int i = 0; i < dedup_size; i++) {
+ for (int j = 0; j < max_natoms+1; j++) {
+ dedup_glove[j][i] = local_mega_glove[j][i];
+ }
+ }
+ } else if (dedup_mode == 1) {
+ for (int i = 0; i < dedup_size; i++) {
+ for (int j = 0; j < max_natoms+1; j++) {
+ dedup_glove[j][i] = global_mega_glove[j][i];
+ }
+ }
+ }
+
+ // 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];
+
+ for (int i = 0; i < dedup_size; i++) {
+ dedup_mask[i] = 0;
+ dup_list[i] = 0;
+ }
+
+ // 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
+ for (int i = dedup_size-1; i > 0; --i) { //dedup_size
+ //choose random entry to swap current one with
+ int k = rand() % (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];
+
+ for (int j = 0; j < max_natoms+1; j++) {
+ dedup_glove[j][i] = dedup_glove[j][k];
+ dedup_glove[j][k] = temp_rxn[j];
+ }
+ }
+
+ for (int i = 0; i < dedup_size; i++) {
+ if (dedup_mask[i] == 0) {
+ int num_dups = 0;
+ int myrxnid1 = dedup_glove[0][i];
+ onemol = atom->molecules[unreacted_mol[myrxnid1]];
+ for (int j = 0; j < onemol->natoms; j++) {
+ int check1 = dedup_glove[j+1][i];
+ for (int ii = i + 1; ii < dedup_size; ii++) {
+ int already_listed = 0;
+ for (int jj = 0; jj < num_dups; jj++) {
+ if (dup_list[jj] == ii) {
+ already_listed = 1;
+ break;
+ }
+ }
+ if (dedup_mask[ii] == 0 && already_listed == 0) {
+ int myrxnid2 = dedup_glove[0][ii];
+ twomol = atom->molecules[unreacted_mol[myrxnid2]];
+ for (int jj = 0; jj < twomol->natoms; jj++) {
+ int check2 = dedup_glove[jj+1][ii];
+ if (check2 == check1) {
+ // add this rxn instance as well
+ if (num_dups == 0) dup_list[num_dups++] = i;
+ dup_list[num_dups++] = ii;
+ break;
+ }
+ }
+ }
+ }
+ }
+ // here we choose random number and therefore reaction instance
+ int myrand = 1;
+ if (num_dups > 0) {
+ myrand = floor(random[0]->uniform()*num_dups);
+ for (int iii = 0; iii < num_dups; iii++) {
+ if (iii != myrand) dedup_mask[dup_list[iii]] = 1;
+ }
+ }
+ }
+ }
+
+ // we must update local_mega_glove and local_megasize
+ // we can simply overwrite local_mega_glove column by column
+ if (dedup_mode == 0) {
+ int new_local_megasize = 0;
+ for (int i = 0; i < local_num_mega; i++) {
+ if (dedup_mask[i] == 0) {
+ local_rxn_count[(int) dedup_glove[0][i]]++;
+ for (int j = 0; j < max_natoms+1; j++) {
+ local_mega_glove[j][new_local_megasize] = dedup_glove[j][i];
+ }
+ new_local_megasize++;
+ }
+ }
+
+ local_num_mega = new_local_megasize;
+ }
+
+ // we must update global_mega_glove and global_megasize
+ // we can simply overwrite global_mega_glove column by column
+ if (dedup_mode == 1) {
+ int new_global_megasize = 0;
+ for (int i = 0; i < global_megasize; i++) {
+ if (dedup_mask[i] == 0) {
+ ghostly_rxn_count[dedup_glove[0][i]]++;
+ for (int j = 0; j < max_natoms + 1; j++) {
+ global_mega_glove[j][new_global_megasize] = dedup_glove[j][i];
+ }
+ new_global_megasize++;
+ }
+ }
+ global_megasize = new_global_megasize;
+ }
+
+}
+
+/* ----------------------------------------------------------------------
+let's limit movement of newly bonded atoms
+and exclude them from other thermostats via exclude_group
+------------------------------------------------------------------------- */
+
+void FixBondReact::limit_bond(int limit_bond_mode)
+{
+ //two types of passes: 1) while superimpose algorithm is iterating (only local atoms)
+ // 2) once more for global_mega_glove [after de-duplicating rxn instances]
+ //in second case, only add local atoms to group
+ //as with update_everything, we can pre-prepare these arrays, then run generic limit_bond code
+
+ //create local, generic variables for onemol->natoms and glove
+ //to be filled differently on respective passes
+
+ int nlocal = atom->nlocal;
+ int temp_limit_num = 0;
+ tagint *temp_limit_glove;
+ if (limit_bond_mode == 0) {
+ int max_temp = local_num_mega * (max_natoms + 1);
+ temp_limit_glove = new tagint[max_temp];
+ for (int j = 0; j < local_num_mega; j++) {
+ rxnID = local_mega_glove[0][j];
+ onemol = atom->molecules[unreacted_mol[rxnID]];
+ for (int i = 0; i < onemol->natoms; i++) {
+ temp_limit_glove[temp_limit_num++] = local_mega_glove[i+1][j];
+ }
+ }
+
+ } else if (limit_bond_mode == 1) {
+ int max_temp = global_megasize * (max_natoms + 1);
+ temp_limit_glove = new tagint[max_temp];
+ for (int j = 0; j < global_megasize; j++) {
+ rxnID = global_mega_glove[0][j];
+ onemol = atom->molecules[unreacted_mol[rxnID]];
+ for (int i = 0; i < onemol->natoms; i++) {
+ if (atom->map(global_mega_glove[i+1][j]) >= 0 &&
+ atom->map(global_mega_glove[i+1][j]) < nlocal)
+ temp_limit_glove[temp_limit_num++] = global_mega_glove[i+1][j];
+ }
+ }
+ }
+
+ if (temp_limit_num == 0) {
+ delete [] temp_limit_glove;
+ return;
+ }
+
+ // we must keep our own list of limited atoms
+ // this will be a new per-atom property!
+
+ int flag;
+ int index1 = atom->find_custom(id_fix2,flag);
+ int *i_limit_tags = atom->ivector[index1];
+
+ int index2 = atom->find_custom(id_fix3,flag);
+ int *i_statted_tags = atom->ivector[index2];
+
+ int index3 = atom->find_custom(id_fix4,flag);
+ int *i_react_tags = atom->ivector[index3];
+
+ for (int i = 0; i < temp_limit_num; i++) {
+ // update->ntimestep could be 0. so add 1 throughout
+ i_limit_tags[atom->map(temp_limit_glove[i])] = update->ntimestep + 1;
+ i_statted_tags[atom->map(temp_limit_glove[i])] = 0;
+ i_react_tags[atom->map(temp_limit_glove[i])] = rxnID;
+ }
+
+ delete [] temp_limit_glove;
+
+}
+
+/* ----------------------------------------------------------------------
+let's unlimit movement of newly bonded atoms after n timesteps.
+we give them back to the system thermostat
+------------------------------------------------------------------------- */
+
+void FixBondReact::unlimit_bond()
+{
+ //let's now unlimit in terms of i_limit_tags
+ //we just run through all nlocal, looking for > limit_duration
+ //then we return i_limit_tag to -1 (which removes from dynamic group)
+ int flag;
+ int index1 = atom->find_custom(id_fix2,flag);
+ int *i_limit_tags = atom->ivector[index1];
+
+ int index2 = atom->find_custom(id_fix3,flag); // statted_tags
+ int *i_statted_tags = atom->ivector[index2];
+
+ int index3 = atom->find_custom(id_fix4,flag); // statted_tags
+ int *i_react_tags = atom->ivector[index3];
+
+ int num2unlimit = 0;
+ for (int i = 0; i < atom->nlocal; i++) {
+ if (i_limit_tags[i] != 0 && (update->ntimestep + 1 - i_limit_tags[i]) > limit_duration[i_react_tags[i]]) {
+ i_limit_tags[i] = 0;
+ i_statted_tags[i] = 1;
+ i_react_tags[i] = 0;
+ }
+ }
+
+ //really should only communicate this per-atom property, not entire reneighboring
+ next_reneighbor = update->ntimestep;
+
+}
+
+/* ----------------------------------------------------------------------
+check mega_glove for ghosts
+if so, flag for broadcasting for perusal by all processors
+------------------------------------------------------------------------- */
+
+void FixBondReact::glove_ghostcheck()
+{
+ // it appears this little loop was deemed important enough for its own function!
+ // noteworthy: it's only relevant for parallel
+
+ // here we add glove to either local_mega_glove or ghostly_mega_glove
+int ghostly = 0;
+ for (int i = 0; i < onemol->natoms; i++) {
+ if (atom->map(glove[i][1]) >= atom->nlocal) {
+ ghostly = 1;
+ break;
+ }
+ }
+
+ if (ghostly == 1) {
+ ghostly_mega_glove[0][ghostly_num_mega] = rxnID;
+ ghostly_rxn_count[rxnID]++; //for debuginng
+ for (int i = 0; i < onemol->natoms; i++) {
+ ghostly_mega_glove[i+1][ghostly_num_mega] = glove[i][1];
+ }
+ ghostly_num_mega++;
+ } else {
+ local_mega_glove[0][local_num_mega] = rxnID;
+ local_rxn_count[rxnID]++; //for debuginng
+ for (int i = 0; i < onemol->natoms; i++) {
+ local_mega_glove[i+1][local_num_mega] = glove[i][1];
+ }
+ local_num_mega++;
+ }
+
+}
+
+/* ----------------------------------------------------------------------
+broadcast entries of mega_glove which contain nonlocal atoms for perusal by all processors
+------------------------------------------------------------------------- */
+
+void FixBondReact::ghost_glovecast()
+{
+#if !defined(MPI_STUBS)
+
+ global_megasize = 0;
+
+ int *allncols = new int[nprocs]();
+ for (int i = 0; i < nprocs; i++)
+ allncols[i] = 0;
+ MPI_Allgather(&ghostly_num_mega, 1, MPI_INT, allncols, 1, MPI_INT, world);
+ for (int i = 0; i < nprocs; i++)
+ global_megasize = global_megasize + allncols[i];
+
+ if (global_megasize == 0) return;
+
+ int *allstarts = new int[nprocs]();
+
+ int start = 0;
+ for (int i = 0; i < me; i++) {
+ start += allncols[i];
+ }
+ MPI_Allgather(&start, 1, MPI_INT, allstarts, 1, MPI_INT, world);
+ MPI_Datatype columnunsized, column;
+ int sizes[2] = {max_natoms+1, global_megasize};
+ int subsizes[2] = {max_natoms+1, 1};
+ int starts[2] = {0,0};
+ MPI_Type_create_subarray (2, sizes, subsizes, starts, MPI_ORDER_C,
+ MPI_LMP_TAGINT, &columnunsized);
+ MPI_Type_create_resized (columnunsized, 0, sizeof(MPI_LMP_TAGINT), &column);
+ MPI_Type_commit(&column);
+
+ if (ghostcheck_flag == 1) memory->destroy(global_mega_glove);
+
+ memory->create(global_mega_glove,max_natoms+1,global_megasize,"bond/react:global_mega_glove");
+
+ for (int i = 0; i < max_natoms+1; i++)
+ 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++) {
+ global_mega_glove[i][j+start] = ghostly_mega_glove[i][j];
+ }
+ }
+ }
+ // let's send to root, dedup, then broadcast
+ MPI_Gatherv(&(global_mega_glove[0][start]), ghostly_num_mega, column,
+ &(global_mega_glove[0][0]), allncols, allstarts,
+ column, 0, world);
+
+ if (me == 0) dedup_mega_gloves(1); // global_mega_glove mode
+ MPI_Bcast(&global_megasize,1,MPI_INT,0,world);
+ MPI_Bcast(&(global_mega_glove[0][0]), global_megasize, column, 0, world);
+
+ delete [] allstarts;
+ delete [] allncols;
+
+#endif
+
+}
+
+/* ----------------------------------------------------------------------
+update charges, types, special lists and all topology
+------------------------------------------------------------------------- */
+
+void FixBondReact::update_everything()
+{
+ tagint *tag = atom->tag;
+ int *type = atom->type;
+
+ int nlocal = atom->nlocal;
+ int nall = nlocal + atom->nghost;
+ int **nspecial = atom->nspecial;
+ tagint **special = atom->special;
+
+ int **bond_type = atom->bond_type;
+ tagint **bond_atom = atom->bond_atom;
+ int *num_bond = atom->num_bond;
+
+ // update atom->nbonds, etc.
+ // TODO: correctly tally with 'newton off'
+ int delta_bonds = 0;
+ int delta_angle = 0;
+ int delta_dihed = 0;
+ int delta_imprp = 0;
+
+ // pass through twice
+ // redefining 'update_num_mega' and 'update_mega_glove' each time
+ // first pass: when glove is all local atoms
+ // second pass: search for local atoms in global_mega_glove
+ // add check for local atoms as well
+
+ int update_num_mega;
+ tagint update_mega_glove[max_natoms+1][MAX(local_num_mega,global_megasize)];
+
+ for (int pass = 0; pass < 2; pass++) {
+
+ if (pass == 0) {
+ update_num_mega = local_num_mega;
+ for (int i = 0; i < update_num_mega; i++) {
+ for (int j = 0; j < max_natoms+1; j++)
+ update_mega_glove[j][i] = local_mega_glove[j][i];
+ }
+ } else if (pass == 1) {
+ update_num_mega = global_megasize;
+ for (int i = 0; i < global_megasize; i++) {
+ for (int j = 0; j < max_natoms+1; j++)
+ update_mega_glove[j][i] = global_mega_glove[j][i];
+ }
+ }
+
+ // update charges and types of landlocked atoms
+ // here, add check for charge instead of requiring it
+ for (int i = 0; i < update_num_mega; i++) {
+ rxnID = update_mega_glove[0][i];
+ twomol = atom->molecules[reacted_mol[rxnID]];
+ for (int j = 0; j < twomol->natoms; j++) {
+ int jj = equivalences[j][1][rxnID]-1;
+ if (landlocked_atoms[j][rxnID] == 1 && atom->map(update_mega_glove[jj+1][i]) >= 0 &&
+ atom->map(update_mega_glove[jj+1][i]) < nlocal) {
+ type[atom->map(update_mega_glove[jj+1][i])] = twomol->type[j];
+ if (twomol->qflag && atom->q_flag) {
+ double *q = atom->q;
+ q[atom->map(update_mega_glove[jj+1][i])] = twomol->q[j];
+ }
+ }
+ }
+ }
+
+ //maybe add check that all 1-3 neighbors of a local atoms are at least ghosts -> unneeded --jg
+ //okay, here goes:
+ for (int i = 0; i < update_num_mega; i++) {
+ rxnID = update_mega_glove[0][i];
+ twomol = atom->molecules[reacted_mol[rxnID]];
+ for (int j = 0; j < twomol->natoms; j++) {
+ int jj = equivalences[j][1][rxnID]-1;
+ if (atom->map(update_mega_glove[jj+1][i]) < nlocal && atom->map(update_mega_glove[jj+1][i]) >= 0) {
+ if (landlocked_atoms[j][rxnID] == 1) {
+ for (int k = 0; k < nspecial[atom->map(update_mega_glove[jj+1][i])][2]; k++) {
+ if (atom->map(special[atom->map(update_mega_glove[jj+1][i])][k]) < 0) {
+ error->all(FLERR,"Fix bond/react needs ghost atoms from further away - most likely too many processors");
+ }
+ }
+ }
+ }
+ }
+ }
+
+ int insert_num;
+ // very nice and easy to completely overwrite special bond info for landlocked atoms
+ for (int i = 0; i < update_num_mega; i++) {
+ rxnID = update_mega_glove[0][i];
+ twomol = atom->molecules[reacted_mol[rxnID]];
+ for (int j = 0; j < twomol->natoms; j++) {
+ int jj = equivalences[j][1][rxnID]-1;
+ if (atom->map(update_mega_glove[jj+1][i]) < nlocal && atom->map(update_mega_glove[jj+1][i]) >= 0) {
+ if (landlocked_atoms[j][rxnID] == 1) {
+ for (int k = 0; k < 3; k++) {
+ nspecial[atom->map(update_mega_glove[jj+1][i])][k] = twomol->nspecial[j][k];
+ }
+ for (int p = 0; p < twomol->nspecial[j][2]; p++) {
+ special[atom->map(update_mega_glove[jj+1][i])][p] = update_mega_glove[equivalences[twomol->special[j][p]-1][1][rxnID]][i];
+ }
+ }
+ // now delete and replace landlocked atoms from non-landlocked atoms' special info
+ if (landlocked_atoms[j][rxnID] == 0) {
+ for (int k = nspecial[atom->map(update_mega_glove[jj+1][i])][2]-1; k > -1; k--) {
+ for (int p = 0; p < twomol->natoms; p++) {
+ int pp = equivalences[p][1][rxnID]-1;
+ if (p!=j && special[atom->map(update_mega_glove[jj+1][i])][k] == update_mega_glove[pp+1][i]
+ && landlocked_atoms[p][rxnID] == 1 ) {
+ for (int n = k; n < nspecial[atom->map(update_mega_glove[jj+1][i])][2]-1; n++) {
+ special[atom->map(update_mega_glove[jj+1][i])][n] = special[atom->map(update_mega_glove[jj+1][i])][n+1];
+ }
+ if (k+1 > nspecial[atom->map(update_mega_glove[jj+1][i])][1]) {
+ nspecial[atom->map(update_mega_glove[jj+1][i])][2]--;
+ } else if (k+1 > nspecial[atom->map(update_mega_glove[jj+1][i])][0]) {
+ nspecial[atom->map(update_mega_glove[jj+1][i])][1]--;
+ nspecial[atom->map(update_mega_glove[jj+1][i])][2]--;
+ } else {
+ nspecial[atom->map(update_mega_glove[jj+1][i])][0]--;
+ nspecial[atom->map(update_mega_glove[jj+1][i])][1]--;
+ nspecial[atom->map(update_mega_glove[jj+1][i])][2]--;
+ }
+ break;
+ }
+ }
+ }
+ // now reassign from reacted template
+ for (int k = 0; k < twomol->nspecial[j][2]; k++) {
+ if (landlocked_atoms[twomol->special[j][k]-1][rxnID] == 1) {
+ if (k > twomol->nspecial[j][1] - 1) {
+ insert_num = nspecial[atom->map(update_mega_glove[jj+1][i])][2]++;
+ } else if (k > twomol->nspecial[j][0] - 1) {
+ insert_num = nspecial[atom->map(update_mega_glove[jj+1][i])][1]++;
+ nspecial[atom->map(update_mega_glove[jj+1][i])][2]++;
+ } else {
+ insert_num = nspecial[atom->map(update_mega_glove[jj+1][i])][0]++;
+ nspecial[atom->map(update_mega_glove[jj+1][i])][1]++;
+ nspecial[atom->map(update_mega_glove[jj+1][i])][2]++;
+ }
+ for (int n = nspecial[atom->map(update_mega_glove[jj+1][i])][2]; n > insert_num; n--) {
+ special[atom->map(update_mega_glove[jj+1][i])][n] = special[atom->map(update_mega_glove[jj+1][i])][n-1];
+ }
+ special[atom->map(update_mega_glove[jj+1][i])][insert_num] = update_mega_glove[equivalences[twomol->special[j][k]-1][1][rxnID]][i];
+ }
+ }
+ }
+ }
+ }
+ }
+
+ // next let's update bond info
+ // cool thing is, newton_bond issues are already taken care of in templates
+ // same with class2 improper issues, which is why this fix started in the first place
+ for (int i = 0; i < update_num_mega; i++) {
+ rxnID = update_mega_glove[0][i];
+ twomol = atom->molecules[reacted_mol[rxnID]];
+ // let's first delete all bond info about landlocked atoms
+ for (int j = 0; j < twomol->natoms; j++) {
+ int jj = equivalences[j][1][rxnID]-1;
+ if (atom->map(update_mega_glove[jj+1][i]) < nlocal && atom->map(update_mega_glove[jj+1][i]) >= 0) {
+ if (landlocked_atoms[j][rxnID] == 1) {
+ delta_bonds -= num_bond[atom->map(update_mega_glove[jj+1][i])];
+ num_bond[atom->map(update_mega_glove[jj+1][i])] = 0;
+ }
+ if (landlocked_atoms[j][rxnID] == 0) {
+ for (int p = num_bond[atom->map(update_mega_glove[jj+1][i])]-1; p > -1 ; p--) {
+ for (int n = 0; n < twomol->natoms; n++) {
+ int nn = equivalences[n][1][rxnID]-1;
+ if (n!=j && bond_atom[atom->map(update_mega_glove[jj+1][i])][p] == update_mega_glove[nn+1][i] && landlocked_atoms[n][rxnID] == 1) {
+ for (int m = p; m < num_bond[atom->map(update_mega_glove[jj+1][i])]-1; m++) {
+ bond_atom[atom->map(update_mega_glove[jj+1][i])][m] = bond_atom[atom->map(update_mega_glove[jj+1][i])][m+1];
+ }
+ num_bond[atom->map(update_mega_glove[jj+1][i])]--;
+ delta_bonds--;
+ }
+ }
+ }
+ }
+ }
+ }
+ // now let's add the new bond info.
+ for (int j = 0; j < twomol->natoms; j++) {
+ int jj = equivalences[j][1][rxnID]-1;
+ if (atom->map(update_mega_glove[jj+1][i]) < nlocal && atom->map(update_mega_glove[jj+1][i]) >= 0) {
+ if (landlocked_atoms[j][rxnID] == 1) {
+ num_bond[atom->map(update_mega_glove[jj+1][i])] = twomol->num_bond[j];
+ delta_bonds += twomol->num_bond[j];
+ for (int p = 0; p < twomol->num_bond[j]; p++) {
+ bond_type[atom->map(update_mega_glove[jj+1][i])][p] = twomol->bond_type[j][p];
+ bond_atom[atom->map(update_mega_glove[jj+1][i])][p] = update_mega_glove[equivalences[twomol->bond_atom[j][p]-1][1][rxnID]][i];
+ }
+ }
+ if (landlocked_atoms[j][rxnID] == 0) {
+ for (int p = 0; p < twomol->num_bond[j]; p++) {
+ if (landlocked_atoms[twomol->bond_atom[j][p]-1][rxnID] == 1) {
+ insert_num = num_bond[atom->map(update_mega_glove[jj+1][i])];
+ bond_type[atom->map(update_mega_glove[jj+1][i])][insert_num] = twomol->bond_type[j][p];
+ bond_atom[atom->map(update_mega_glove[jj+1][i])][insert_num] = update_mega_glove[equivalences[twomol->bond_atom[j][p]-1][1][rxnID]][i];
+ num_bond[atom->map(update_mega_glove[jj+1][i])]++;
+ delta_bonds++;
+ }
+ }
+ }
+ }
+ }
+ }
+
+ if (force->angle && twomol->angleflag) {
+ int *num_angle = atom->num_angle;
+ int **angle_type = atom->angle_type;
+ tagint **angle_atom1 = atom->angle_atom1;
+ tagint **angle_atom2 = atom->angle_atom2;
+ tagint **angle_atom3 = atom->angle_atom3;
+
+ for (int i = 0; i < update_num_mega; i++) {
+ rxnID = update_mega_glove[0][i];
+ twomol = atom->molecules[reacted_mol[rxnID]];
+ // Angles! First let's delete all angle info:
+ for (int j = 0; j < twomol->natoms; j++) {
+ int jj = equivalences[j][1][rxnID]-1;
+ if (atom->map(update_mega_glove[jj+1][i]) < nlocal && atom->map(update_mega_glove[jj+1][i]) >= 0) {
+ if (landlocked_atoms[j][rxnID] == 1) {
+ delta_angle -= num_angle[atom->map(update_mega_glove[jj+1][i])];
+ num_angle[atom->map(update_mega_glove[jj+1][i])] = 0;
+ }
+ if (landlocked_atoms[j][rxnID] == 0) {
+ for (int p = num_angle[atom->map(update_mega_glove[jj+1][i])]-1; p > -1; p--) {
+ for (int n = 0; n < twomol->natoms; n++) {
+ int nn = equivalences[n][1][rxnID]-1;
+ if (n!=j && landlocked_atoms[n][rxnID] == 1 &&
+ (angle_atom1[atom->map(update_mega_glove[jj+1][i])][p] == update_mega_glove[nn+1][i] ||
+ angle_atom2[atom->map(update_mega_glove[jj+1][i])][p] == update_mega_glove[nn+1][i] ||
+ angle_atom3[atom->map(update_mega_glove[jj+1][i])][p] == update_mega_glove[nn+1][i])) {
+ for (int m = p; m < num_angle[atom->map(update_mega_glove[jj+1][i])]-1; m++) {
+ angle_atom1[atom->map(update_mega_glove[jj+1][i])][m] = angle_atom1[atom->map(update_mega_glove[jj+1][i])][m+1];
+ angle_atom2[atom->map(update_mega_glove[jj+1][i])][m] = angle_atom2[atom->map(update_mega_glove[jj+1][i])][m+1];
+ angle_atom3[atom->map(update_mega_glove[jj+1][i])][m] = angle_atom3[atom->map(update_mega_glove[jj+1][i])][m+1];
+ }
+ num_angle[atom->map(update_mega_glove[jj+1][i])]--;
+ delta_angle--;
+ break;
+ }
+ }
+ }
+ }
+ }
+ }
+ // now let's add the new angle info.
+ for (int j = 0; j < twomol->natoms; j++) {
+ int jj = equivalences[j][1][rxnID]-1;
+ if (atom->map(update_mega_glove[jj+1][i]) < nlocal && atom->map(update_mega_glove[jj+1][i]) >= 0) {
+ if (landlocked_atoms[j][rxnID] == 1) {
+ num_angle[atom->map(update_mega_glove[jj+1][i])] = twomol->num_angle[j];
+ delta_angle += twomol->num_angle[j];
+ for (int p = 0; p < twomol->num_angle[j]; p++) {
+ angle_type[atom->map(update_mega_glove[jj+1][i])][p] = twomol->angle_type[j][p];
+ angle_atom1[atom->map(update_mega_glove[jj+1][i])][p] = update_mega_glove[equivalences[twomol->angle_atom1[j][p]-1][1][rxnID]][i];
+ angle_atom2[atom->map(update_mega_glove[jj+1][i])][p] = update_mega_glove[equivalences[twomol->angle_atom2[j][p]-1][1][rxnID]][i];
+ angle_atom3[atom->map(update_mega_glove[jj+1][i])][p] = update_mega_glove[equivalences[twomol->angle_atom3[j][p]-1][1][rxnID]][i];
+ }
+ }
+ if (landlocked_atoms[j][rxnID] == 0) {
+ for (int p = 0; p < twomol->num_angle[j]; p++) {
+ if (landlocked_atoms[twomol->angle_atom1[j][p]-1][rxnID] == 1 ||
+ landlocked_atoms[twomol->angle_atom2[j][p]-1][rxnID] == 1 ||
+ landlocked_atoms[twomol->angle_atom3[j][p]-1][rxnID] == 1) {
+ insert_num = num_angle[atom->map(update_mega_glove[jj+1][i])];
+ angle_type[atom->map(update_mega_glove[jj+1][i])][insert_num] = twomol->angle_type[j][p];
+ angle_atom1[atom->map(update_mega_glove[jj+1][i])][insert_num] = update_mega_glove[equivalences[twomol->angle_atom1[j][p]-1][1][rxnID]][i];
+ angle_atom2[atom->map(update_mega_glove[jj+1][i])][insert_num] = update_mega_glove[equivalences[twomol->angle_atom2[j][p]-1][1][rxnID]][i];
+ angle_atom3[atom->map(update_mega_glove[jj+1][i])][insert_num] = update_mega_glove[equivalences[twomol->angle_atom3[j][p]-1][1][rxnID]][i];
+ num_angle[atom->map(update_mega_glove[jj+1][i])]++;
+ delta_angle++;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ // Dihedrals! first let's delete all dihedral info for landlocked atoms
+ if (force->dihedral && twomol->dihedralflag) {
+ int *num_dihedral = atom->num_dihedral;
+ int **dihedral_type = atom->dihedral_type;
+ tagint **dihedral_atom1 = atom->dihedral_atom1;
+ tagint **dihedral_atom2 = atom->dihedral_atom2;
+ tagint **dihedral_atom3 = atom->dihedral_atom3;
+ tagint **dihedral_atom4 = atom->dihedral_atom4;
+
+ for (int i = 0; i < update_num_mega; i++) {
+ rxnID = update_mega_glove[0][i];
+ twomol = atom->molecules[reacted_mol[rxnID]];
+ for (int j = 0; j < twomol->natoms; j++) {
+ int jj = equivalences[j][1][rxnID]-1;
+ if (atom->map(update_mega_glove[jj+1][i]) < nlocal && atom->map(update_mega_glove[jj+1][i]) >= 0) {
+ if (landlocked_atoms[j][rxnID] == 1) {
+ delta_dihed -= num_dihedral[atom->map(update_mega_glove[jj+1][i])];
+ num_dihedral[atom->map(update_mega_glove[jj+1][i])] = 0;
+ }
+ if (landlocked_atoms[j][rxnID] == 0) {
+ for (int p = num_dihedral[atom->map(update_mega_glove[jj+1][i])]-1; p > -1; p--) {
+ for (int n = 0; n < twomol->natoms; n++) {
+ int nn = equivalences[n][1][rxnID]-1;
+ if (n!=j && landlocked_atoms[n][rxnID] == 1 &&
+ (dihedral_atom1[atom->map(update_mega_glove[jj+1][i])][p] == update_mega_glove[nn+1][i] ||
+ dihedral_atom2[atom->map(update_mega_glove[jj+1][i])][p] == update_mega_glove[nn+1][i] ||
+ dihedral_atom3[atom->map(update_mega_glove[jj+1][i])][p] == update_mega_glove[nn+1][i] ||
+ dihedral_atom4[atom->map(update_mega_glove[jj+1][i])][p] == update_mega_glove[nn+1][i])) {
+ for (int m = p; m < num_dihedral[atom->map(update_mega_glove[jj+1][i])]-1; m++) {
+ dihedral_atom1[atom->map(update_mega_glove[jj+1][i])][m] = dihedral_atom1[atom->map(update_mega_glove[jj+1][i])][m+1];
+ dihedral_atom2[atom->map(update_mega_glove[jj+1][i])][m] = dihedral_atom2[atom->map(update_mega_glove[jj+1][i])][m+1];
+ dihedral_atom3[atom->map(update_mega_glove[jj+1][i])][m] = dihedral_atom3[atom->map(update_mega_glove[jj+1][i])][m+1];
+ dihedral_atom4[atom->map(update_mega_glove[jj+1][i])][m] = dihedral_atom4[atom->map(update_mega_glove[jj+1][i])][m+1];
+ }
+ num_dihedral[atom->map(update_mega_glove[jj+1][i])]--;
+ delta_dihed--;
+ break;
+ }
+ }
+ }
+ }
+ }
+ }
+ // now let's add new dihedral info
+ for (int j = 0; j < twomol->natoms; j++) {
+ int jj = equivalences[j][1][rxnID]-1;
+ if (atom->map(update_mega_glove[jj+1][i]) < nlocal && atom->map(update_mega_glove[jj+1][i]) >= 0) {
+ if (landlocked_atoms[j][rxnID] == 1) {
+ num_dihedral[atom->map(update_mega_glove[jj+1][i])] = twomol->num_dihedral[j];
+ delta_dihed += twomol->num_dihedral[j];
+ for (int p = 0; p < twomol->num_dihedral[j]; p++) {
+ dihedral_type[atom->map(update_mega_glove[jj+1][i])][p] = twomol->dihedral_type[j][p];
+ dihedral_atom1[atom->map(update_mega_glove[jj+1][i])][p] = update_mega_glove[equivalences[twomol->dihedral_atom1[j][p]-1][1][rxnID]][i];
+ dihedral_atom2[atom->map(update_mega_glove[jj+1][i])][p] = update_mega_glove[equivalences[twomol->dihedral_atom2[j][p]-1][1][rxnID]][i];
+ dihedral_atom3[atom->map(update_mega_glove[jj+1][i])][p] = update_mega_glove[equivalences[twomol->dihedral_atom3[j][p]-1][1][rxnID]][i];
+ dihedral_atom4[atom->map(update_mega_glove[jj+1][i])][p] = update_mega_glove[equivalences[twomol->dihedral_atom4[j][p]-1][1][rxnID]][i];
+ }
+ }
+ if (landlocked_atoms[j][rxnID] == 0) {
+ for (int p = 0; p < twomol->num_dihedral[j]; p++) {
+ if (landlocked_atoms[twomol->dihedral_atom1[j][p]-1][rxnID] == 1 ||
+ landlocked_atoms[twomol->dihedral_atom2[j][p]-1][rxnID] == 1 ||
+ landlocked_atoms[twomol->dihedral_atom3[j][p]-1][rxnID] == 1 ||
+ landlocked_atoms[twomol->dihedral_atom4[j][p]-1][rxnID] == 1) {
+ insert_num = num_dihedral[atom->map(update_mega_glove[jj+1][i])];
+ dihedral_type[atom->map(update_mega_glove[jj+1][i])][insert_num] = twomol->dihedral_type[j][p];
+ dihedral_atom1[atom->map(update_mega_glove[jj+1][i])][insert_num] = update_mega_glove[equivalences[twomol->dihedral_atom1[j][p]-1][1][rxnID]][i];
+ dihedral_atom2[atom->map(update_mega_glove[jj+1][i])][insert_num] = update_mega_glove[equivalences[twomol->dihedral_atom2[j][p]-1][1][rxnID]][i];
+ dihedral_atom3[atom->map(update_mega_glove[jj+1][i])][insert_num] = update_mega_glove[equivalences[twomol->dihedral_atom3[j][p]-1][1][rxnID]][i];
+ dihedral_atom4[atom->map(update_mega_glove[jj+1][i])][insert_num] = update_mega_glove[equivalences[twomol->dihedral_atom4[j][p]-1][1][rxnID]][i];
+ num_dihedral[atom->map(update_mega_glove[jj+1][i])]++;
+ delta_dihed++;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ // finally IMPROPERS!!!! first let's delete all improper info for landlocked atoms
+ if (force->improper && twomol->improperflag) {
+ int *num_improper = atom->num_improper;
+ int **improper_type = atom->improper_type;
+ tagint **improper_atom1 = atom->improper_atom1;
+ tagint **improper_atom2 = atom->improper_atom2;
+ tagint **improper_atom3 = atom->improper_atom3;
+ tagint **improper_atom4 = atom->improper_atom4;
+
+ for (int i = 0; i < update_num_mega; i++) {
+ rxnID = update_mega_glove[0][i];
+ twomol = atom->molecules[reacted_mol[rxnID]];
+ for (int j = 0; j < twomol->natoms; j++) {
+ int jj = equivalences[j][1][rxnID]-1;
+ if (atom->map(update_mega_glove[jj+1][i]) < nlocal && atom->map(update_mega_glove[jj+1][i]) >= 0) {
+ if (landlocked_atoms[j][rxnID] == 1) {
+ delta_imprp -= num_improper[atom->map(update_mega_glove[jj+1][i])];
+ num_improper[atom->map(update_mega_glove[jj+1][i])] = 0;
+ }
+ if (landlocked_atoms[j][rxnID] == 0) {
+ for (int p = num_improper[atom->map(update_mega_glove[jj+1][i])]-1; p > -1; p--) {
+ for (int n = 0; n < twomol->natoms; n++) {
+ int nn = equivalences[n][1][rxnID]-1;
+ if (n!=j && landlocked_atoms[n][rxnID] == 1 &&
+ (improper_atom1[atom->map(update_mega_glove[jj+1][i])][p] == update_mega_glove[nn+1][i] ||
+ improper_atom2[atom->map(update_mega_glove[jj+1][i])][p] == update_mega_glove[nn+1][i] ||
+ improper_atom3[atom->map(update_mega_glove[jj+1][i])][p] == update_mega_glove[nn+1][i] ||
+ improper_atom4[atom->map(update_mega_glove[jj+1][i])][p] == update_mega_glove[nn+1][i])) {
+ for (int m = p; m < num_improper[atom->map(update_mega_glove[jj+1][i])]-1; m++) {
+ improper_atom1[atom->map(update_mega_glove[jj+1][i])][m] = improper_atom1[atom->map(update_mega_glove[jj+1][i])][m+1];
+ improper_atom2[atom->map(update_mega_glove[jj+1][i])][m] = improper_atom2[atom->map(update_mega_glove[jj+1][i])][m+1];
+ improper_atom3[atom->map(update_mega_glove[jj+1][i])][m] = improper_atom3[atom->map(update_mega_glove[jj+1][i])][m+1];
+ improper_atom4[atom->map(update_mega_glove[jj+1][i])][m] = improper_atom4[atom->map(update_mega_glove[jj+1][i])][m+1];
+ }
+ num_improper[atom->map(update_mega_glove[jj+1][i])]--;
+ delta_imprp--;
+ break;
+ }
+ }
+ }
+ }
+ }
+ }
+ // now let's add new improper info
+ for (int j = 0; j < twomol->natoms; j++) {
+ int jj = equivalences[j][1][rxnID]-1;
+ if (atom->map(update_mega_glove[jj+1][i]) < nlocal && atom->map(update_mega_glove[jj+1][i]) >= 0) {
+ if (landlocked_atoms[j][rxnID] == 1) {
+ num_improper[atom->map(update_mega_glove[jj+1][i])] = twomol->num_improper[j];
+ delta_imprp += twomol->num_improper[j];
+ for (int p = 0; p < twomol->num_improper[j]; p++) {
+ improper_type[atom->map(update_mega_glove[jj+1][i])][p] = twomol->improper_type[j][p];
+ improper_atom1[atom->map(update_mega_glove[jj+1][i])][p] = update_mega_glove[equivalences[twomol->improper_atom1[j][p]-1][1][rxnID]][i];
+ improper_atom2[atom->map(update_mega_glove[jj+1][i])][p] = update_mega_glove[equivalences[twomol->improper_atom2[j][p]-1][1][rxnID]][i];
+ improper_atom3[atom->map(update_mega_glove[jj+1][i])][p] = update_mega_glove[equivalences[twomol->improper_atom3[j][p]-1][1][rxnID]][i];
+ improper_atom4[atom->map(update_mega_glove[jj+1][i])][p] = update_mega_glove[equivalences[twomol->improper_atom4[j][p]-1][1][rxnID]][i];
+ }
+ }
+ if (landlocked_atoms[j][rxnID] == 0) {
+ for (int p = 0; p < twomol->num_improper[j]; p++) {
+ if (landlocked_atoms[twomol->improper_atom1[j][p]-1][rxnID] == 1 ||
+ landlocked_atoms[twomol->improper_atom2[j][p]-1][rxnID] == 1 ||
+ landlocked_atoms[twomol->improper_atom3[j][p]-1][rxnID] == 1 ||
+ landlocked_atoms[twomol->improper_atom4[j][p]-1][rxnID] == 1) {
+ insert_num = num_improper[atom->map(update_mega_glove[jj+1][i])];
+ improper_type[atom->map(update_mega_glove[jj+1][i])][insert_num] = twomol->improper_type[j][p];
+ improper_atom1[atom->map(update_mega_glove[jj+1][i])][insert_num] = update_mega_glove[equivalences[twomol->improper_atom1[j][p]-1][1][rxnID]][i];
+ improper_atom2[atom->map(update_mega_glove[jj+1][i])][insert_num] = update_mega_glove[equivalences[twomol->improper_atom2[j][p]-1][1][rxnID]][i];
+ improper_atom3[atom->map(update_mega_glove[jj+1][i])][insert_num] = update_mega_glove[equivalences[twomol->improper_atom3[j][p]-1][1][rxnID]][i];
+ improper_atom4[atom->map(update_mega_glove[jj+1][i])][insert_num] = update_mega_glove[equivalences[twomol->improper_atom4[j][p]-1][1][rxnID]][i];
+ num_improper[atom->map(update_mega_glove[jj+1][i])]++;
+ delta_imprp++;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ }
+ // 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
+
+ int Tdelta_bonds;
+ MPI_Allreduce(&delta_bonds,&Tdelta_bonds,1,MPI_LMP_BIGINT,MPI_SUM,world);
+ atom->nbonds += Tdelta_bonds;
+
+ int Tdelta_angle;
+ MPI_Allreduce(&delta_angle,&Tdelta_angle,1,MPI_LMP_BIGINT,MPI_SUM,world);
+ atom->nangles += Tdelta_angle;
+
+ int Tdelta_dihed;
+ MPI_Allreduce(&delta_dihed,&Tdelta_dihed,1,MPI_LMP_BIGINT,MPI_SUM,world);
+ atom->ndihedrals += Tdelta_dihed;
+
+ int Tdelta_imprp;
+ MPI_Allreduce(&delta_imprp,&Tdelta_imprp,1,MPI_LMP_BIGINT,MPI_SUM,world);
+ atom->nimpropers += Tdelta_imprp;
+
+
+}
+
+/* ----------------------------------------------------------------------
+read superimpose file
+------------------------------------------------------------------------- */
+
+void FixBondReact::read(int myrxn)
+{
+ int i;
+ char line[MAXLINE],keyword[MAXLINE];
+ char *eof,*ptr;
+
+ // skip 1st line of file
+ eof = fgets(line,MAXLINE,fp);
+ if (eof == NULL) error->one(FLERR,"Unexpected end of superimpose file");
+
+ // read header lines
+ // skip blank lines or lines that start with "#"
+ // stop when read an unrecognized line
+
+ while (1) {
+
+ readline(line);
+
+ // trim anything from '#' onward
+ // if line is blank, continue
+
+ if ((ptr = strchr(line,'#'))) *ptr = '\0';
+ 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 break;
+ }
+
+ //count = NULL;
+
+ // grab keyword and skip next line
+
+ parse_keyword(0,line,keyword);
+ readline(line);
+
+ // loop over sections of superimpose file
+
+ int equivflag = 0, edgeflag = 0, bondflag = 0;
+ while (strlen(keyword)) {
+ if (strcmp(keyword,"BondingIDs") == 0) {
+ bondflag = 1;
+ readline(line);
+ sscanf(line,"%d",&ibonding[myrxn]);
+ readline(line);
+ sscanf(line,"%d",&jbonding[myrxn]);
+ } else if (strcmp(keyword,"EdgeIDs") == 0) {
+ edgeflag = 1;
+ for (int i = 0; i < onemol->natoms; i++) edge[i][myrxn] = 0;
+ EdgeIDs(line, myrxn);
+ } else if (strcmp(keyword,"Equivalences") == 0) {
+ equivflag = 1;
+ Equivalences(line, myrxn);
+ } else error->one(FLERR,"Unknown section in superimpose file");
+
+ parse_keyword(1,line,keyword);
+
+ }
+
+ // error check
+ if (bondflag == 0 || equivflag == 0 || edgeflag == 0)
+ error->all(FLERR,"Superimpose file missing BondingIDs, EdgeIDs, or Equivalences section\n");
+}
+
+void FixBondReact::EdgeIDs(char *line, int myrxn)
+{
+ // puts a 1 at edge(edgeID)
+
+ int tmp;
+ for (int i = 0; i < nedge; i++) {
+ readline(line);
+ sscanf(line,"%d",&tmp);
+ edge[tmp-1][myrxn] = 1;
+ }
+}
+
+void FixBondReact::Equivalences(char *line, int myrxn)
+{
+ int tmp1;
+ int tmp2;
+ for (int i = 0; i < nequivalent; i++) {
+ readline(line);
+ sscanf(line,"%d %d",&tmp1,&tmp2);
+ //equivalences is-> clmn 1: post-reacted, clmn 2: pre-reacted
+ equivalences[tmp2-1][0][myrxn] = tmp2;
+ equivalences[tmp2-1][1][myrxn] = tmp1;
+ //reverse_equiv is-> clmn 1: pre-reacted, clmn 2: post-reacted
+ reverse_equiv[tmp1-1][0][myrxn] = tmp1;
+ reverse_equiv[tmp1-1][1][myrxn] = tmp2;
+ }
+}
+
+void FixBondReact::open(char *file)
+{
+ fp = fopen(file,"r");
+ if (fp == NULL) {
+ char str[128];
+ sprintf(str,"Cannot open superimpose file %s",file);
+ error->one(FLERR,str);
+ }
+}
+
+void FixBondReact::readline(char *line)
+{
+ int n;
+ if (me == 0) {
+ if (fgets(line,MAXLINE,fp) == NULL) n = 0;
+ else n = strlen(line) + 1;
+ }
+ MPI_Bcast(&n,1,MPI_INT,0,world);
+ if (n == 0) error->all(FLERR,"Unexpected end of superimpose file");
+ MPI_Bcast(line,n,MPI_CHAR,0,world);
+}
+
+void FixBondReact::parse_keyword(int flag, char *line, char *keyword)
+{
+ if (flag) {
+
+ // read upto non-blank line plus 1 following line
+ // eof is set to 1 if any read hits end-of-file
+
+ int eof = 0;
+ if (me == 0) {
+ if (fgets(line,MAXLINE,fp) == NULL) eof = 1;
+ while (eof == 0 && strspn(line," \t\n\r") == strlen(line)) {
+ if (fgets(line,MAXLINE,fp) == NULL) eof = 1;
+ }
+ if (fgets(keyword,MAXLINE,fp) == NULL) eof = 1;
+ }
+
+ // if eof, set keyword empty and return
+
+ MPI_Bcast(&eof,1,MPI_INT,0,world);
+ if (eof) {
+ keyword[0] = '\0';
+ return;
+ }
+
+ // bcast keyword line to all procs
+
+ int n;
+ if (me == 0) n = strlen(line) + 1;
+ MPI_Bcast(&n,1,MPI_INT,0,world);
+ MPI_Bcast(line,n,MPI_CHAR,0,world);
+ }
+
+ // copy non-whitespace portion of line into keyword
+
+ int start = strspn(line," \t\n\r");
+ int stop = strlen(line) - 1;
+ while (line[stop] == ' ' || line[stop] == '\t'
+ || line[stop] == '\n' || line[stop] == '\r') stop--;
+ line[stop+1] = '\0';
+ strcpy(keyword,&line[start]);
+}
+
+
+void FixBondReact::skip_lines(int n, char *line)
+{
+ for (int i = 0; i < n; i++) readline(line);
+}
+
+int FixBondReact::parse(char *line, char **words, int max)
+{
+ char *ptr;
+
+ int nwords = 0;
+ words[nwords++] = strtok(line," \t\n\r\f");
+
+ while ((ptr = strtok(NULL," \t\n\r\f"))) {
+ if (nwords < max) words[nwords] = ptr;
+ nwords++;
+ }
+
+ return nwords;
+}
+
+/* ---------------------------------------------------------------------- */
+
+double FixBondReact::compute_vector(int n)
+{
+ // now we print just the totals for each reaction instance
+ return (double) reaction_count_total[n];
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixBondReact::post_integrate_respa(int ilevel, int iloop)
+{
+ if (ilevel == nlevels_respa-1) post_integrate();
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixBondReact::pack_forward_comm(int n, int *list, double *buf,
+int pbc_flag, int *pbc)
+{
+
+ int i,j,k,m,ns;
+
+ m = 0;
+
+ if (commflag == 1) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ printf("hello you shouldn't be here\n");
+ //buf[m++] = ubuf(bondcount[j]).d;
+ }
+ return m;
+ }
+
+ if (commflag == 2) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = ubuf(partner[j]).d;
+ buf[m++] = probability[j];
+ }
+ return m;
+ }
+
+ int **nspecial = atom->nspecial;
+ tagint **special = atom->special;
+
+ m = 0;
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = ubuf(finalpartner[j]).d;
+ ns = nspecial[j][0];
+ buf[m++] = ubuf(ns).d;
+ for (k = 0; k < ns; k++)
+ buf[m++] = ubuf(special[j][k]).d;
+ }
+ return m;
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixBondReact::unpack_forward_comm(int n, int first, double *buf)
+{
+
+ int i,j,m,ns,last;
+
+ m = 0;
+ last = first + n;
+
+ if (commflag == 1) {
+ for (i = first; i < last; i++)
+ printf("hello you shouldn't be here\n");
+ // bondcount[i] = (int) ubuf(buf[m++]).i;
+
+ } else if (commflag == 2) {
+ for (i = first; i < last; i++) {
+ partner[i] = (tagint) ubuf(buf[m++]).i;
+ probability[i] = buf[m++];
+ }
+
+ } else {
+ int **nspecial = atom->nspecial;
+ tagint **special = atom->special;
+
+ m = 0;
+ last = first + n;
+ for (i = first; i < last; i++) {
+ finalpartner[i] = (tagint) ubuf(buf[m++]).i;
+ ns = (int) ubuf(buf[m++]).i;
+ nspecial[i][0] = ns;
+ for (j = 0; j < ns; j++)
+ special[i][j] = (tagint) ubuf(buf[m++]).i;
+ }
+ }
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixBondReact::pack_reverse_comm(int n, int first, double *buf)
+{
+ int i,m,last;
+
+ m = 0;
+ last = first + n;
+
+ for (i = first; i < last; i++) {
+ buf[m++] = ubuf(partner[i]).d;
+ buf[m++] = distsq[i];
+ }
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixBondReact::unpack_reverse_comm(int n, int *list, double *buf)
+{
+ int i,j,m;
+
+ m = 0;
+
+ if (commflag != 1) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ if (buf[m+1] < distsq[j]) {
+ partner[j] = (tagint) ubuf(buf[m++]).i;
+ distsq[j] = buf[m++];
+ } else m += 2;
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+memory usage of local atom-based arrays
+------------------------------------------------------------------------- */
+
+double FixBondReact::memory_usage()
+{
+ int nmax = atom->nmax;
+ double bytes = nmax * sizeof(int);
+ bytes = 2*nmax * sizeof(tagint);
+ bytes += nmax * sizeof(double);
+ return bytes;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixBondReact::print_bb()
+{
+
+ //fix bond/create cargo code. eg nbonds needs to be added
+ /*
+for (int i = 0; i < atom->nlocal; i++) {
+ // printf("TAG " TAGINT_FORMAT ": %d nbonds: ",atom->tag[i],atom->num_bond[i]);
+ for (int j = 0; j < atom->num_bond[i]; j++) {
+ // printf(" " TAGINT_FORMAT,atom->bond_atom[i][j]);
+ }
+ // printf("\n");
+ // printf("TAG " TAGINT_FORMAT ": %d nangles: ",atom->tag[i],atom->num_angle[i]);
+ for (int j = 0; j < atom->num_angle[i]; j++) {
+ // printf(" " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT ",",
+ atom->angle_atom1[i][j], atom->angle_atom2[i][j],
+ atom->angle_atom3[i][j]);
+ }
+ // printf("\n");
+ // printf("TAG " TAGINT_FORMAT ": %d ndihedrals: ",atom->tag[i],atom->num_dihedral[i]);
+ for (int j = 0; j < atom->num_dihedral[i]; j++) {
+ // printf(" " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " "
+ TAGINT_FORMAT ",", atom->dihedral_atom1[i][j],
+ atom->dihedral_atom2[i][j],atom->dihedral_atom3[i][j],
+ atom->dihedral_atom4[i][j]);
+ }
+ // printf("\n");
+ // printf("TAG " TAGINT_FORMAT ": %d nimpropers: ",atom->tag[i],atom->num_improper[i]);
+ for (int j = 0; j < atom->num_improper[i]; j++) {
+ // printf(" " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " "
+ TAGINT_FORMAT ",",atom->improper_atom1[i][j],
+ atom->improper_atom2[i][j],atom->improper_atom3[i][j],
+ atom->improper_atom4[i][j]);
+ }
+ // printf("\n");
+ // printf("TAG " TAGINT_FORMAT ": %d %d %d nspecial: ",atom->tag[i],
+ atom->nspecial[i][0],atom->nspecial[i][1],atom->nspecial[i][2]);
+ for (int j = 0; j < atom->nspecial[i][2]; j++) {
+ // printf(" " TAGINT_FORMAT,atom->special[i][j]);
+ }
+ // printf("\n");
+}
+*/
+}
diff --git a/src/fix_bond_react.h b/src/fix_bond_react.h
new file mode 100644
index 0000000000000000000000000000000000000000..b235f152dfb3c00547bf6934c81ca8eff6f7d9fa
--- /dev/null
+++ b/src/fix_bond_react.h
@@ -0,0 +1,218 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing Author: Jacob Gissinger (jacob.gissinger@colorado.edu)
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(bond/react,FixBondReact)
+
+#else
+
+#ifndef LMP_FIX_BOND_REACT_H
+#define LMP_FIX_BOND_REACT_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixBondReact : public Fix {
+ public:
+ FixBondReact(class LAMMPS *, int, char **);
+ ~FixBondReact();
+ int setmask();
+ void post_constructor();
+ void init();
+ void init_list(int, class NeighList *);
+ void post_integrate();
+ void post_integrate_respa(int, int);
+
+ int pack_forward_comm(int, int *, double *, int, int *);
+ void unpack_forward_comm(int, int, double *);
+ int pack_reverse_comm(int, int, double *);
+ void unpack_reverse_comm(int, int *, double *);
+ double compute_vector(int);
+ double memory_usage();
+
+ private:
+ int me,nprocs;
+ int nreacts;
+ int *nevery;
+ FILE *fp;
+ int *iatomtype,*jatomtype;
+ int *seed;
+ double *cutsq,*fraction;
+ tagint lastcheck;
+ int stabilization_flag;
+ int *stabilize_steps_flag;
+ int status;
+ int *groupbits;
+
+ int rxnID; // integer ID for identifying current bond/react
+ int *reaction_count;
+ int *reaction_count_total;
+ int nmax; // max num local atoms
+ int max_natoms; // max natoms in a molecule template
+ tagint *partner,*finalpartner;
+ double *distsq,*probability;
+ int *ncreate;
+ int maxcreate;
+ int allncreate;
+ tagint ***created;
+ int *local_ncreate;
+
+ class Molecule *onemol; // pre-reacted molecule template
+ class Molecule *twomol; // post-reacted molecule template
+ class FixNVELimit *fix1; // used to relax reaction sites
+ class FixPropertyAtom *fix2; // used to indicate relaxing atoms
+ class FixPropertyAtom *fix3; // used to indicate system-wide thermostat
+ class FixPropertyAtom *fix4; // indicates to which 'react' atom belongs
+ class RanMars **random;
+ class NeighList *list;
+
+ int *reacted_mol,*unreacted_mol;
+ int *limit_duration; // indicates how long to relax
+ 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 property (recently reacted)
+ char *id_fix3; // id of internally created fix per-atom property (system-wide thermostat)
+ char *id_fix4; // id of internally created fix per-atom property (ID of 'react' argument)
+ char *master_group; // group containing relaxing atoms from all fix rxns
+ char *exclude_group; // group for system-wide thermostat
+
+ int countflag,commflag;
+ int nlevels_respa;
+
+ void superimpose_algorithm(); // main function of the superimpose algorithm
+
+ int *ibonding,*jbonding;
+ 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
+ int *guess_branch; // used when there is more than two choices when guessing
+ int **restore_pt; // contains info about restore points
+ tagint **restore; // contaings info about restore points
+ int *pioneer_count; // counts pioneers
+
+ int **edge; // atoms in molecule templates with incorrect valences
+ int ***equivalences; // relation between pre- and post-reacted templates
+ int ***reverse_equiv; // re-ordered equivalences
+ int **landlocked_atoms; // all atoms at least three bonds away from edge atoms
+
+ int pion,neigh,trace; // important indices for various loops. required for restore points
+ int lcl_inst; // reaction instance
+ tagint **glove; // 1st colmn: pre-reacted template, 2nd colmn: global IDs
+ // for all mega_gloves and global_mega_glove: first row is the ID of bond/react
+ tagint **local_mega_glove; // consolidation local of reaction instances
+ tagint **ghostly_mega_glove; // consolidation nonlocal of reaction instances
+ tagint **global_mega_glove; // consolidation (inter-processor) of gloves containing nonlocal atoms
+ int local_num_mega; // num of local reaction instances
+ int ghostly_num_mega; // num of ghostly reaction instances
+ int global_megasize; // num of reaction instances in global_mega_glove
+ int *pioneers; // during Superimpose Algorithm, atoms which have been assigned, but whose first neighbors haven't
+ int glove_counter; // used to determine when to terminate Superimpose Algorithm
+
+ void read(int);
+ void EdgeIDs(char *,int);
+ void Equivalences(char *,int);
+
+ void make_a_guess ();
+ void neighbor_loop();
+ void check_a_neighbor();
+ void crosscheck_the_neighbor();
+ void inner_crosscheck_loop();
+ void ring_check();
+
+ void open(char *);
+ void readline(char *);
+ void parse_keyword(int, char *, char *);
+ void skip_lines(int, char *);
+ int parse(char *, char **, int);
+
+ void find_landlocked_atoms(int);
+ void glove_ghostcheck();
+ void ghost_glovecast();
+ void update_everything();
+ void unlimit_bond();
+ void limit_bond(int);
+ void dedup_mega_gloves(int); //dedup global mega_glove
+
+ // DEBUG (currently obsolete)
+
+ void print_bb();
+
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory. Check the input script syntax and compare to the
+documentation for the command. You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Invalid exclude group name
+
+Exclude group name should not previously be defined.
+
+E: Cannot use fix bond/react with non-molecular systems
+
+Only systems with bonds that can be changed can be used. Atom_style
+template does not qualify.
+
+E: Fix bond/react cutoff is longer than pairwise cutoff
+
+This is not allowed because bond creation is done using the
+pairwise neighbor list.
+
+E: Molecule template ID for fix bond/react does not exist
+
+A valid molecule template must have been created with the molecule command.
+
+E: Superimpose file errors:
+
+Please ensure superimpose file is properly formatted.
+
+E: Atom affected by reaction too close to template edge
+
+This means an atom which changes type during the reaction is too close
+to an 'edge' atom defined in the superimpose file. This could cause incorrect
+assignment of bonds, angle, etc. Generally, this means you must include
+more atoms in your templates, such that there are at least two atoms
+between each atom involved in the reaction and an edge atom.
+
+E: Fix bond/react needs ghost atoms from farther away
+
+This is because a processor needs to superimpose the entire unreacted
+molecule template onto simulation atoms it can 'see.' The comm_modify cutoff
+command can be used to extend the communication range.
+
+E: Excessive iteration of superimpose algorithm
+
+You may have discovered a bug! But first, please double check that your
+molecule template atom types, bond types, etc. are consistent with your simulation,
+and that all atoms affected by a reaction are sufficently separated from edge atoms.
+If this issue persists, please contact the developer.
+
+*/
diff --git a/src/fix_group.cpp b/src/fix_group.cpp
index fe2495d6c0061a71c4a8a576c0a06acc0b5fa169..51abebdf292adc8c89e00489f197a03441c8b6ae 100644
--- a/src/fix_group.cpp
+++ b/src/fix_group.cpp
@@ -33,7 +33,7 @@ using namespace FixConst;
/* ---------------------------------------------------------------------- */
FixGroup::FixGroup(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg),
-idregion(NULL), idvar(NULL)
+idregion(NULL), idvar(NULL), idprop(NULL)
{
// dgroupbit = bitmask of dynamic group
// group ID is last part of fix ID
@@ -49,6 +49,7 @@ idregion(NULL), idvar(NULL)
regionflag = 0;
varflag = 0;
+ propflag = 0;
nevery = 1;
int iarg = 3;
@@ -73,7 +74,17 @@ idregion(NULL), idvar(NULL)
idvar = new char[n];
strcpy(idvar,arg[iarg+1]);
iarg += 2;
- } else if (strcmp(arg[iarg],"every") == 0) {
+ } else if (strcmp(arg[iarg],"property") == 0) {
+ if (iarg+2 > narg) error->all(FLERR,"Illegal group command");
+ if (atom->find_custom(arg[iarg+1],typeflag) < 0)
+ error->all(FLERR,"Per atom property for group dynamic does not exist");
+ propflag = 1;
+ delete [] idprop;
+ int n = strlen(arg[iarg+1]) + 1;
+ idprop = new char[n];
+ strcpy(idprop,arg[iarg+1]);
+ iarg += 2;
+ } else if (strcmp(arg[iarg],"every") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal group command");
nevery = force->inumeric(FLERR,arg[iarg+1]);
if (nevery <= 0) error->all(FLERR,"Illegal group command");
@@ -88,6 +99,7 @@ FixGroup::~FixGroup()
{
delete [] idregion;
delete [] idvar;
+ delete [] idprop;
}
/* ---------------------------------------------------------------------- */
@@ -130,6 +142,12 @@ void FixGroup::init()
error->all(FLERR,"Variable for group dynamic is invalid style");
}
+ if (propflag) {
+ iprop = atom->find_custom(idprop,typeflag);
+ if (iprop < 0)
+ error->all(FLERR,"Per-atom property for group dynamic does not exist");
+ }
+
// warn if any FixGroup is not at tail end of all post_integrate fixes
Fix **fix = modify->fix;
@@ -188,6 +206,9 @@ void FixGroup::set_group()
// invoke atom-style variable if defined
double *var = NULL;
+ int *ivector = NULL;
+ double *dvector = NULL;
+
if (varflag) {
modify->clearstep_compute();
@@ -196,6 +217,12 @@ void FixGroup::set_group()
modify->addstep_compute(update->ntimestep + nevery);
}
+ // invoke per-atom property if defined
+
+ if (propflag && !typeflag) ivector = atom->ivector[iprop]; //check nlocal > 0?
+
+ if (propflag && typeflag) dvector = atom->dvector[iprop]; //check nlocal > 0?
+
// update region in case it has a variable dependence or is dynamic
if (regionflag) region->prematch();
@@ -214,6 +241,8 @@ void FixGroup::set_group()
inflag = 1;
if (regionflag && !region->match(x[i][0],x[i][1],x[i][2])) inflag = 0;
if (varflag && var[i] == 0.0) inflag = 0;
+ if (propflag && !typeflag && ivector[i] == 0) inflag = 0;
+ if (propflag && typeflag && dvector[i] == 0) inflag = 0;
} else inflag = 0;
if (inflag) mask[i] |= gbit;
diff --git a/src/fix_group.h b/src/fix_group.h
index 6ed842578d5ea1834c84c74e3804c2b43f0416b6..662325492bf068aacfb5ff36d725d339b92272f4 100644
--- a/src/fix_group.h
+++ b/src/fix_group.h
@@ -36,9 +36,9 @@ class FixGroup : public Fix {
private:
int gbit,gbitinverse;
- int regionflag,varflag;
- int iregion,ivar;
- char *idregion,*idvar;
+ int regionflag,varflag,propflag,typeflag;
+ int iregion,ivar,iprop;
+ char *idregion,*idvar,*idprop;
class Region *region;
int nlevels_respa;
diff --git a/src/fix_nve_limit.cpp b/src/fix_nve_limit.cpp
index 966fcfbb09ca0530b5893d7e160a519fc86d9384..6f85d107f113c178332c03450d55a77d0ee44747 100644
--- a/src/fix_nve_limit.cpp
+++ b/src/fix_nve_limit.cpp
@@ -38,6 +38,7 @@ FixNVELimit::FixNVELimit(LAMMPS *lmp, int narg, char **arg) :
scalar_flag = 1;
global_freq = 1;
extscalar = 1;
+ dynamic_group_allow = 1;
xlimit = force->numeric(FLERR,arg[3]);
diff --git a/src/fix_temp_rescale.cpp b/src/fix_temp_rescale.cpp
index fe52a317880f371d647d89ee90f8f5b7ea9273eb..cc16a9a6e616cb80e9998db0098826cecd35b30d 100644
--- a/src/fix_temp_rescale.cpp
+++ b/src/fix_temp_rescale.cpp
@@ -37,11 +37,12 @@ enum{CONSTANT,EQUAL};
/* ---------------------------------------------------------------------- */
FixTempRescale::FixTempRescale(LAMMPS *lmp, int narg, char **arg) :
- Fix(lmp, narg, arg),
- tstr(NULL), id_temp(NULL), tflag(0)
+ Fix(lmp, narg, arg)
{
if (narg < 8) error->all(FLERR,"Illegal fix temp/rescale command");
+ dynamic_group_allow = 1;
+
nevery = force->inumeric(FLERR,arg[3]);
if (nevery <= 0) error->all(FLERR,"Illegal fix temp/rescale command");
@@ -139,7 +140,7 @@ void FixTempRescale::end_of_step()
if (temperature->dof < 1) return;
- // protect against division by zero
+ // protect against division by zero.
if (t_current == 0.0)
error->all(FLERR,"Computed temperature for fix temp/rescale cannot be 0.0");
diff --git a/src/write_data.cpp b/src/write_data.cpp
index bf00249e51f2ba13c91b7f9616991b0c293d3f3a..7c8838628eb482f298c20046a40f8f0bab514a30 100644
--- a/src/write_data.cpp
+++ b/src/write_data.cpp
@@ -243,11 +243,11 @@ void WriteData::header()
fprintf(fp,"%d angle types\n",atom->nangletypes);
}
if (atom->ndihedrals || atom->ndihedraltypes) {
- fprintf(fp,BIGINT_FORMAT " dihedrals\n",atom->ndihedrals);
+ fprintf(fp,BIGINT_FORMAT " dihedrals\n",ndihedrals);
fprintf(fp,"%d dihedral types\n",atom->ndihedraltypes);
}
if (atom->nimpropers || atom->nimpropertypes) {
- fprintf(fp,BIGINT_FORMAT " impropers\n",atom->nimpropers);
+ fprintf(fp,BIGINT_FORMAT " impropers\n",nimpropers);
fprintf(fp,"%d improper types\n",atom->nimpropertypes);
}
}