diff --git a/doc/src/Section_commands.txt b/doc/src/Section_commands.txt
index 76051490e401b9ff5b9ac3c4aeb292a1b86b66b3..cdbdec9ee3ef6e7a8d13cbb28b53d9d41034c487 100644
--- a/doc/src/Section_commands.txt
+++ b/doc/src/Section_commands.txt
@@ -684,6 +684,7 @@ package"_Section_start.html#start_3.
"addtorque"_fix_addtorque.html,
"atc"_fix_atc.html,
"ave/correlate/long"_fix_ave_correlate_long.html,
+"bond/react"_fix_bond_react.html,
"colvars"_fix_colvars.html,
"dpd/energy (k)"_fix_dpd_energy.html,
"drude"_fix_drude.html,
diff --git a/doc/src/fix_bond_react.txt b/doc/src/fix_bond_react.txt
new file mode 100644
index 0000000000000000000000000000000000000000..b5f9a730842b73805d146d4dc07b2ee0b2c00062
--- /dev/null
+++ b/doc/src/fix_bond_react.txt
@@ -0,0 +1,356 @@
+"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.html :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 :pre
+
+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:
+
+:line
+
+# This is a map file :pre
+
+2 edgeIDs
+7 equivalences :pre
+
+BondingIDs :pre
+
+3 5 :pre
+
+EdgeIDs :pre
+
+1 7 :pre
+
+Equivalences :pre
+
+1 1
+2 2
+3 3
+4 4
+5 5
+6 6
+7 7 :pre
+
+:line
+
+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 its
+own 'react' argument. 2) While typically a bond is formed between the
+bonding atom pairs specified in the pre-reacted molecule template,
+this is not required.
+
+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
+reaction step. 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 steps.
+
+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 USER-MISC 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/doc/src/fixes.txt b/doc/src/fixes.txt
index c3fcc06bf1c53e780b9b7bc3ce9885671fc9f690..fca442cc13836eb33ad422ec2e0a7cfe3b88a948 100644
--- a/doc/src/fixes.txt
+++ b/doc/src/fixes.txt
@@ -23,6 +23,7 @@ Fixes :h1
fix_bond_break
fix_bond_create
fix_bond_swap
+ fix_bond_react
fix_box_relax
fix_cmap
fix_colvars
diff --git a/doc/src/lammps.book b/doc/src/lammps.book
index ec34f41872db2aa143628ce3f99c249a8641cda9..73f8d005be530d4039912f8d6c4d171471bb5b92 100644
--- a/doc/src/lammps.book
+++ b/doc/src/lammps.book
@@ -137,6 +137,7 @@ fix_aveforce.html
fix_balance.html
fix_bond_break.html
fix_bond_create.html
+fix_bond_react.html
fix_bond_swap.html
fix_box_relax.html
fix_cmap.html
diff --git a/examples/USER/misc/bond_react/nylon,6-6_melt/in.large_nylon_melt b/examples/USER/misc/bond_react/nylon,6-6_melt/in.large_nylon_melt
new file mode 100644
index 0000000000000000000000000000000000000000..072a8e3c456af5547403b35d859b08a5c512f941
--- /dev/null
+++ b/examples/USER/misc/bond_react/nylon,6-6_melt/in.large_nylon_melt
@@ -0,0 +1,52 @@
+# 35,000 atom nylon melt example
+
+units real
+
+boundary p p p
+
+atom_style full
+
+kspace_style pppm 1.0e-4
+
+pair_style lj/class2/coul/long 8.5
+
+angle_style class2
+
+bond_style class2
+
+dihedral_style class2
+
+improper_style class2
+
+read_data large_nylon_melt.data.gz
+
+velocity all create 800.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 100 test_vis.xyz
+
+fix myrxns all bond/react stabilization yes statted_grp .03 &
+ react rxn1 all 1 2.9 mol1 mol2 rxn1_stp1_map &
+ react rxn2 all 1 5 mol3 mol4 rxn1_stp2_map
+
+# stable at 800K
+fix 1 statted_grp nvt temp 800 800 100
+
+# in order to customize behavior of reacting atoms,
+# you can use the internally created 'bond_react_MASTER_group', like so:
+# fix 2 bond_react_MASTER_group temp/rescale 1 800 800 10 1
+
+thermo_style custom step temp press density f_myrxns[1] f_myrxns[2] # cumulative reaction counts
+
+# restart 100 restart1 restart2
+
+run 200
+
+# write_restart restart_longrun
+# write_data restart_longrun.data
diff --git a/examples/USER/misc/bond_react/nylon,6-6_melt/large_nylon_melt.data.gz b/examples/USER/misc/bond_react/nylon,6-6_melt/large_nylon_melt.data.gz
new file mode 100644
index 0000000000000000000000000000000000000000..c620b879a883817b8eb7abe4495e585ecacfaf02
Binary files /dev/null and b/examples/USER/misc/bond_react/nylon,6-6_melt/large_nylon_melt.data.gz differ
diff --git a/examples/USER/misc/bond_react/nylon,6-6_melt/log.20Apr18.large_nylon_melt.g++.1 b/examples/USER/misc/bond_react/nylon,6-6_melt/log.20Apr18.large_nylon_melt.g++.1
new file mode 100644
index 0000000000000000000000000000000000000000..f3b3840a925573ab1b83db1402cd68c5579dab35
--- /dev/null
+++ b/examples/USER/misc/bond_react/nylon,6-6_melt/log.20Apr18.large_nylon_melt.g++.1
@@ -0,0 +1,170 @@
+LAMMPS (20 Apr 2018)
+OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:90)
+ using 1 OpenMP thread(s) per MPI task
+# 35,000 atom nylon melt example
+
+units real
+
+boundary p p p
+
+atom_style full
+
+kspace_style pppm 1.0e-4
+
+pair_style lj/class2/coul/long 8.5
+
+angle_style class2
+
+bond_style class2
+
+dihedral_style class2
+
+improper_style class2
+
+read_data large_nylon_melt.data.gz
+ orthogonal box = (-2.68344 -2.06791 -2.21988) to (73.4552 73.2448 73.4065)
+ 1 by 1 by 1 MPI processor grid
+ reading atoms ...
+ 35200 atoms
+ reading velocities ...
+ 35200 velocities
+ scanning bonds ...
+ 9 = max bonds/atom
+ scanning angles ...
+ 21 = max angles/atom
+ scanning dihedrals ...
+ 31 = max dihedrals/atom
+ scanning impropers ...
+ 29 = max impropers/atom
+ reading bonds ...
+ 33600 bonds
+ reading angles ...
+ 59200 angles
+ reading dihedrals ...
+ 80000 dihedrals
+ reading impropers ...
+ 35200 impropers
+ 4 = max # of 1-2 neighbors
+ 6 = max # of 1-3 neighbors
+ 12 = max # of 1-4 neighbors
+ 41 = max # of special neighbors
+
+velocity all create 800.0 4928459 dist gaussian
+
+molecule mol1 rxn1_stp1_unreacted.data_template
+Read molecule mol1:
+ 18 atoms with max type 8
+ 16 bonds with max type 12
+ 25 angles with max type 24
+ 23 dihedrals with max type 33
+ 14 impropers with max type 9
+molecule mol2 rxn1_stp1_reacted.data_template
+Read molecule mol2:
+ 18 atoms with max type 9
+ 17 bonds with max type 11
+ 31 angles with max type 23
+ 39 dihedrals with max type 30
+ 20 impropers with max type 1
+molecule mol3 rxn1_stp2_unreacted.data_template
+Read molecule mol3:
+ 15 atoms with max type 9
+ 14 bonds with max type 11
+ 25 angles with max type 23
+ 30 dihedrals with max type 30
+ 16 impropers with max type 1
+molecule mol4 rxn1_stp2_reacted.data_template
+Read molecule mol4:
+ 15 atoms with max type 11
+ 13 bonds with max type 13
+ 19 angles with max type 25
+ 16 dihedrals with max type 29
+ 10 impropers with max type 11
+
+thermo 50
+
+# dump 1 all xyz 100 test_vis.xyz
+
+fix myrxns all bond/react stabilization yes statted_grp .03 react rxn1 all 1 2.9 mol1 mol2 rxn1_stp1_map react rxn2 all 1 5 mol3 mol4 rxn1_stp2_map
+dynamic group bond_react_MASTER_group defined
+dynamic group statted_grp defined
+dynamic group bond_react_MASTER_group defined
+dynamic group statted_grp defined
+
+# stable at 800K
+fix 1 statted_grp nvt temp 800 800 100
+
+# in order to customize behavior of reacting atoms,
+# you can use the internally created 'bond_react_MASTER_group', like so:
+# fix 2 bond_react_MASTER_group temp/rescale 1 800 800 10 1
+
+thermo_style custom step temp press density f_myrxns[1] f_myrxns[2] # cumulative reaction counts
+
+# restart 100 restart1 restart2
+
+run 200
+PPPM initialization ...
+ using 12-bit tables for long-range coulomb (../kspace.cpp:321)
+ G vector (1/distance) = 0.20765
+ grid = 18 18 18
+ stencil order = 5
+ estimated absolute RMS force accuracy = 0.0333156
+ estimated relative force accuracy = 0.000100329
+ using double precision FFTs
+ 3d grid and FFT values/proc = 12167 5832
+Neighbor list info ...
+ update every 1 steps, delay 10 steps, check yes
+ max neighbors/atom: 2000, page size: 100000
+ master list distance cutoff = 10.5
+ ghost atom cutoff = 10.5
+ binsize = 5.25, bins = 15 15 15
+ 2 neighbor lists, perpetual/occasional/extra = 1 1 0
+ (1) pair lj/class2/coul/long, perpetual
+ attributes: half, newton on
+ pair build: half/bin/newton
+ stencil: half/bin/3d/newton
+ bin: standard
+ (2) fix bond/react, occasional, copy from (1)
+ attributes: half, newton on
+ pair build: copy
+ stencil: none
+ bin: none
+Per MPI rank memory allocation (min/avg/max) = 209.1 | 209.1 | 209.1 Mbytes
+Step Temp Press Density f_myrxns[1] f_myrxns[2]
+ 0 800 3666.3948 0.80366765 0 0
+ 50 673.95238 -9670.9169 0.80366765 31 0
+ 100 693.69241 -4696.4359 0.80366765 57 22
+ 150 715.44689 -14740.892 0.80366765 77 50
+ 200 721.16898 -1411.95 0.80366765 84 66
+Loop time of 107.389 on 1 procs for 200 steps with 35200 atoms
+
+Performance: 0.161 ns/day, 149.151 hours/ns, 1.862 timesteps/s
+99.7% CPU use with 1 MPI tasks x 1 OpenMP threads
+
+MPI task timing breakdown:
+Section | min time | avg time | max time |%varavg| %total
+---------------------------------------------------------------
+Pair | 27.191 | 27.191 | 27.191 | 0.0 | 25.32
+Bond | 11.46 | 11.46 | 11.46 | 0.0 | 10.67
+Kspace | 4.2507 | 4.2507 | 4.2507 | 0.0 | 3.96
+Neigh | 55.544 | 55.544 | 55.544 | 0.0 | 51.72
+Comm | 0.41715 | 0.41715 | 0.41715 | 0.0 | 0.39
+Output | 0.0011044 | 0.0011044 | 0.0011044 | 0.0 | 0.00
+Modify | 8.4756 | 8.4756 | 8.4756 | 0.0 | 7.89
+Other | | 0.04897 | | | 0.05
+
+Nlocal: 35200 ave 35200 max 35200 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Nghost: 38406 ave 38406 max 38406 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Neighs: 6.92787e+06 ave 6.92787e+06 max 6.92787e+06 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+
+Total # of neighbors = 6927872
+Ave neighs/atom = 196.815
+Ave special neighs/atom = 9.83489
+Neighbor list builds = 200
+Dangerous builds = 0
+
+# write_restart restart_longrun
+# write_data restart_longrun.data
+Total wall time: 0:01:48
diff --git a/examples/USER/misc/bond_react/nylon,6-6_melt/log.20Apr18.large_nylon_melt.g++.4 b/examples/USER/misc/bond_react/nylon,6-6_melt/log.20Apr18.large_nylon_melt.g++.4
new file mode 100644
index 0000000000000000000000000000000000000000..992a97e77ef920433575efd080df965045460154
--- /dev/null
+++ b/examples/USER/misc/bond_react/nylon,6-6_melt/log.20Apr18.large_nylon_melt.g++.4
@@ -0,0 +1,170 @@
+LAMMPS (20 Apr 2018)
+OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:90)
+ using 1 OpenMP thread(s) per MPI task
+# 35,000 atom nylon melt example
+
+units real
+
+boundary p p p
+
+atom_style full
+
+kspace_style pppm 1.0e-4
+
+pair_style lj/class2/coul/long 8.5
+
+angle_style class2
+
+bond_style class2
+
+dihedral_style class2
+
+improper_style class2
+
+read_data large_nylon_melt.data.gz
+ orthogonal box = (-2.68344 -2.06791 -2.21988) to (73.4552 73.2448 73.4065)
+ 2 by 1 by 2 MPI processor grid
+ reading atoms ...
+ 35200 atoms
+ reading velocities ...
+ 35200 velocities
+ scanning bonds ...
+ 9 = max bonds/atom
+ scanning angles ...
+ 21 = max angles/atom
+ scanning dihedrals ...
+ 31 = max dihedrals/atom
+ scanning impropers ...
+ 29 = max impropers/atom
+ reading bonds ...
+ 33600 bonds
+ reading angles ...
+ 59200 angles
+ reading dihedrals ...
+ 80000 dihedrals
+ reading impropers ...
+ 35200 impropers
+ 4 = max # of 1-2 neighbors
+ 6 = max # of 1-3 neighbors
+ 12 = max # of 1-4 neighbors
+ 41 = max # of special neighbors
+
+velocity all create 800.0 4928459 dist gaussian
+
+molecule mol1 rxn1_stp1_unreacted.data_template
+Read molecule mol1:
+ 18 atoms with max type 8
+ 16 bonds with max type 12
+ 25 angles with max type 24
+ 23 dihedrals with max type 33
+ 14 impropers with max type 9
+molecule mol2 rxn1_stp1_reacted.data_template
+Read molecule mol2:
+ 18 atoms with max type 9
+ 17 bonds with max type 11
+ 31 angles with max type 23
+ 39 dihedrals with max type 30
+ 20 impropers with max type 1
+molecule mol3 rxn1_stp2_unreacted.data_template
+Read molecule mol3:
+ 15 atoms with max type 9
+ 14 bonds with max type 11
+ 25 angles with max type 23
+ 30 dihedrals with max type 30
+ 16 impropers with max type 1
+molecule mol4 rxn1_stp2_reacted.data_template
+Read molecule mol4:
+ 15 atoms with max type 11
+ 13 bonds with max type 13
+ 19 angles with max type 25
+ 16 dihedrals with max type 29
+ 10 impropers with max type 11
+
+thermo 50
+
+# dump 1 all xyz 100 test_vis.xyz
+
+fix myrxns all bond/react stabilization yes statted_grp .03 react rxn1 all 1 2.9 mol1 mol2 rxn1_stp1_map react rxn2 all 1 5 mol3 mol4 rxn1_stp2_map
+dynamic group bond_react_MASTER_group defined
+dynamic group statted_grp defined
+dynamic group bond_react_MASTER_group defined
+dynamic group statted_grp defined
+
+# stable at 800K
+fix 1 statted_grp nvt temp 800 800 100
+
+# in order to customize behavior of reacting atoms,
+# you can use the internally created 'bond_react_MASTER_group', like so:
+# fix 2 bond_react_MASTER_group temp/rescale 1 800 800 10 1
+
+thermo_style custom step temp press density f_myrxns[1] f_myrxns[2] # cumulative reaction counts
+
+# restart 100 restart1 restart2
+
+run 200
+PPPM initialization ...
+ using 12-bit tables for long-range coulomb (../kspace.cpp:321)
+ G vector (1/distance) = 0.20765
+ grid = 18 18 18
+ stencil order = 5
+ estimated absolute RMS force accuracy = 0.0333156
+ estimated relative force accuracy = 0.000100329
+ using double precision FFTs
+ 3d grid and FFT values/proc = 4508 1620
+Neighbor list info ...
+ update every 1 steps, delay 10 steps, check yes
+ max neighbors/atom: 2000, page size: 100000
+ master list distance cutoff = 10.5
+ ghost atom cutoff = 10.5
+ binsize = 5.25, bins = 15 15 15
+ 2 neighbor lists, perpetual/occasional/extra = 1 1 0
+ (1) pair lj/class2/coul/long, perpetual
+ attributes: half, newton on
+ pair build: half/bin/newton
+ stencil: half/bin/3d/newton
+ bin: standard
+ (2) fix bond/react, occasional, copy from (1)
+ attributes: half, newton on
+ pair build: copy
+ stencil: none
+ bin: none
+Per MPI rank memory allocation (min/avg/max) = 81.11 | 81.13 | 81.15 Mbytes
+Step Temp Press Density f_myrxns[1] f_myrxns[2]
+ 0 800 3666.3948 0.80366765 0 0
+ 50 673.95238 -9670.9169 0.80366765 31 0
+ 100 693.69241 -4696.4359 0.80366765 57 22
+ 150 715.43654 -14742.205 0.80366765 77 50
+ 200 721.1906 -1411.4303 0.80366765 84 66
+Loop time of 56.2311 on 4 procs for 200 steps with 35200 atoms
+
+Performance: 0.307 ns/day, 78.099 hours/ns, 3.557 timesteps/s
+99.1% CPU use with 4 MPI tasks x 1 OpenMP threads
+
+MPI task timing breakdown:
+Section | min time | avg time | max time |%varavg| %total
+---------------------------------------------------------------
+Pair | 13.86 | 14.034 | 14.406 | 5.8 | 24.96
+Bond | 5.5592 | 5.5952 | 5.6492 | 1.4 | 9.95
+Kspace | 2.3969 | 2.7523 | 2.9203 | 12.5 | 4.89
+Neigh | 27.265 | 27.268 | 27.271 | 0.0 | 48.49
+Comm | 0.75523 | 0.77355 | 0.79381 | 1.7 | 1.38
+Output | 0.00051904 | 0.0007363 | 0.0013669 | 0.0 | 0.00
+Modify | 5.7629 | 5.7634 | 5.7641 | 0.0 | 10.25
+Other | | 0.04441 | | | 0.08
+
+Nlocal: 8800 ave 8912 max 8666 min
+Histogram: 1 0 0 1 0 0 0 0 1 1
+Nghost: 18358.8 ave 18432 max 18189 min
+Histogram: 1 0 0 0 0 0 0 0 1 2
+Neighs: 1.73197e+06 ave 1.77209e+06 max 1.68475e+06 min
+Histogram: 1 0 1 0 0 0 0 0 0 2
+
+Total # of neighbors = 6927873
+Ave neighs/atom = 196.815
+Ave special neighs/atom = 9.83489
+Neighbor list builds = 200
+Dangerous builds = 0
+
+# write_restart restart_longrun
+# write_data restart_longrun.data
+Total wall time: 0:00:57
diff --git a/examples/USER/misc/bond_react/nylon,6-6_melt/rxn1_stp1_map b/examples/USER/misc/bond_react/nylon,6-6_melt/rxn1_stp1_map
new file mode 100644
index 0000000000000000000000000000000000000000..44f7ad8137549eaa375046e114c59a7aa09c0c65
--- /dev/null
+++ b/examples/USER/misc/bond_react/nylon,6-6_melt/rxn1_stp1_map
@@ -0,0 +1,35 @@
+this is a nominal 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/USER/misc/bond_react/nylon,6-6_melt/rxn1_stp1_reacted.data_template b/examples/USER/misc/bond_react/nylon,6-6_melt/rxn1_stp1_reacted.data_template
new file mode 100644
index 0000000000000000000000000000000000000000..61c0408ce3e57955e736d15e77736b95dccf240b
--- /dev/null
+++ b/examples/USER/misc/bond_react/nylon,6-6_melt/rxn1_stp1_reacted.data_template
@@ -0,0 +1,189 @@
+this is a molecule template for: initial nylon crosslink, post-reacting
+
+18 atoms
+17 bonds
+31 angles
+39 dihedrals
+20 impropers
+
+Types
+
+1 9
+2 1
+3 1
+4 4
+5 4
+6 3
+7 3
+8 1
+9 1
+10 5
+11 8
+12 6
+13 3
+14 3
+15 7
+16 1
+17 3
+18 3
+
+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 9 1 2
+2 10 1 4
+3 10 1 5
+4 11 1 10
+5 1 2 3
+6 2 2 6
+7 2 2 7
+8 1 3 16
+9 2 3 17
+10 2 3 18
+11 1 8 9
+12 6 9 10
+13 2 9 13
+14 2 9 14
+15 7 10 11
+16 5 10 12
+17 8 12 15
+
+Angles
+
+1 14 2 1 4
+2 14 2 1 5
+3 15 2 1 10
+4 16 4 1 5
+5 17 4 1 10
+6 17 5 1 10
+7 18 1 2 3
+8 19 1 2 6
+9 19 1 2 7
+10 1 3 2 6
+11 1 3 2 7
+12 3 6 2 7
+13 2 2 3 16
+14 1 2 3 17
+15 1 2 3 18
+16 1 16 3 17
+17 1 16 3 18
+18 3 17 3 18
+19 12 8 9 10
+20 1 8 9 13
+21 1 8 9 14
+22 13 13 9 10
+23 13 14 9 10
+24 3 13 9 14
+25 10 9 10 11
+26 8 9 10 12
+27 20 1 10 9
+28 21 11 10 12
+29 22 1 10 11
+30 23 1 10 12
+31 11 10 12 15
+
+Dihedrals
+
+1 16 4 1 2 3
+2 17 4 1 2 6
+3 17 4 1 2 7
+4 16 5 1 2 3
+5 17 5 1 2 6
+6 17 5 1 2 7
+7 18 10 1 2 3
+8 19 10 1 2 6
+9 19 10 1 2 7
+10 20 2 1 10 9
+11 21 2 1 10 11
+12 22 2 1 10 12
+13 23 4 1 10 9
+14 24 4 1 10 11
+15 25 4 1 10 12
+16 23 5 1 10 9
+17 24 5 1 10 11
+18 25 5 1 10 12
+19 26 1 2 3 16
+20 27 1 2 3 17
+21 27 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 14 8 9 10 11
+29 12 8 9 10 12
+30 28 8 9 10 1
+31 15 13 9 10 11
+32 13 13 9 10 12
+33 29 13 9 10 1
+34 15 14 9 10 11
+35 13 14 9 10 12
+36 29 14 9 10 1
+37 10 9 10 12 15
+38 11 11 10 12 15
+39 30 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/USER/misc/bond_react/nylon,6-6_melt/rxn1_stp1_unreacted.data_template b/examples/USER/misc/bond_react/nylon,6-6_melt/rxn1_stp1_unreacted.data_template
new file mode 100644
index 0000000000000000000000000000000000000000..944d6918c57e77d2d96b3a949efac38cfce1c5f6
--- /dev/null
+++ b/examples/USER/misc/bond_react/nylon,6-6_melt/rxn1_stp1_unreacted.data_template
@@ -0,0 +1,160 @@
+this is a molecule template for: initial nylon crosslink, pre-reacting
+
+18 atoms
+16 bonds
+25 angles
+23 dihedrals
+14 impropers
+
+Types
+
+1 2
+2 1
+3 1
+4 4
+5 4
+6 3
+7 3
+8 1
+9 1
+10 5
+11 8
+12 6
+13 3
+14 3
+15 7
+16 1
+17 3
+18 3
+
+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 12 1 2
+2 4 1 4
+3 4 1 5
+4 1 2 3
+5 2 2 6
+6 2 2 7
+7 1 3 16
+8 2 3 17
+9 2 3 18
+10 1 8 9
+11 6 9 10
+12 2 9 13
+13 2 9 14
+14 7 10 11
+15 5 10 12
+16 8 12 15
+
+Angles
+
+1 6 2 1 4
+2 6 2 1 5
+3 7 4 1 5
+4 24 1 2 3
+5 5 1 2 6
+6 5 1 2 7
+7 1 3 2 6
+8 1 3 2 7
+9 3 6 2 7
+10 2 2 3 16
+11 1 2 3 17
+12 1 2 3 18
+13 1 16 3 17
+14 1 16 3 18
+15 3 17 3 18
+16 12 8 9 10
+17 1 8 9 13
+18 1 8 9 14
+19 13 13 9 10
+20 13 14 9 10
+21 3 13 9 14
+22 10 9 10 11
+23 8 9 10 12
+24 21 11 10 12
+25 11 10 12 15
+
+Dihedrals
+
+1 31 4 1 2 3
+2 32 4 1 2 6
+3 32 4 1 2 7
+4 31 5 1 2 3
+5 32 5 1 2 6
+6 32 5 1 2 7
+7 33 1 2 3 16
+8 1 1 2 3 17
+9 1 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 14 8 9 10 11
+17 12 8 9 10 12
+18 15 13 9 10 11
+19 13 13 9 10 12
+20 15 14 9 10 11
+21 13 14 9 10 12
+22 10 9 10 12 15
+23 11 11 10 12 15
+
+Impropers
+
+1 1 2 1 4 5
+2 9 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/USER/misc/bond_react/nylon,6-6_melt/rxn1_stp2_map b/examples/USER/misc/bond_react/nylon,6-6_melt/rxn1_stp2_map
new file mode 100644
index 0000000000000000000000000000000000000000..35fe47fdb3bdb0386d65b950d45cbfb39fbe6f7f
--- /dev/null
+++ b/examples/USER/misc/bond_react/nylon,6-6_melt/rxn1_stp2_map
@@ -0,0 +1,32 @@
+this is a nominal 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/USER/misc/bond_react/nylon,6-6_melt/rxn1_stp2_reacted.data_template b/examples/USER/misc/bond_react/nylon,6-6_melt/rxn1_stp2_reacted.data_template
new file mode 100644
index 0000000000000000000000000000000000000000..ffd3ef733c0ff4eb28c59ba4d8c657a2710dd8ad
--- /dev/null
+++ b/examples/USER/misc/bond_react/nylon,6-6_melt/rxn1_stp2_reacted.data_template
@@ -0,0 +1,131 @@
+this is a molecule template for: water condensation, post-reacting
+
+15 atoms
+13 bonds
+19 angles
+16 dihedrals
+10 impropers
+
+Types
+
+1 9
+2 1
+3 1
+4 10
+5 4
+6 3
+7 3
+8 1
+9 1
+10 5
+11 8
+12 11
+13 3
+14 3
+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 9 1 2
+2 10 1 5
+3 11 1 10
+4 1 2 3
+5 2 2 6
+6 2 2 7
+7 13 4 12
+8 1 8 9
+9 6 9 10
+10 2 9 13
+11 2 9 14
+12 7 10 11
+13 13 15 12
+
+Angles
+
+1 14 2 1 5
+2 15 2 1 10
+3 17 5 1 10
+4 18 1 2 3
+5 19 1 2 6
+6 19 1 2 7
+7 1 3 2 6
+8 1 3 2 7
+9 3 6 2 7
+10 12 8 9 10
+11 1 8 9 13
+12 1 8 9 14
+13 13 13 9 10
+14 13 14 9 10
+15 3 13 9 14
+16 10 9 10 11
+17 20 1 10 9
+18 22 1 10 11
+19 25 15 12 4
+
+Dihedrals
+
+1 16 5 1 2 3
+2 17 5 1 2 6
+3 17 5 1 2 7
+4 18 10 1 2 3
+5 19 10 1 2 6
+6 19 10 1 2 7
+7 20 2 1 10 9
+8 21 2 1 10 11
+9 23 5 1 10 9
+10 24 5 1 10 11
+11 14 8 9 10 11
+12 28 8 9 10 1
+13 15 13 9 10 11
+14 29 13 9 10 1
+15 15 14 9 10 11
+16 29 14 9 10 1
+
+Impropers
+
+1 10 2 1 5 10
+2 11 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/USER/misc/bond_react/nylon,6-6_melt/rxn1_stp2_unreacted.data_template b/examples/USER/misc/bond_react/nylon,6-6_melt/rxn1_stp2_unreacted.data_template
new file mode 100644
index 0000000000000000000000000000000000000000..7abe15ada8d2372e4435ef82b5e60b7417865790
--- /dev/null
+++ b/examples/USER/misc/bond_react/nylon,6-6_melt/rxn1_stp2_unreacted.data_template
@@ -0,0 +1,158 @@
+this is a molecule template for: water condensation, pre-reacting
+
+15 atoms
+14 bonds
+25 angles
+30 dihedrals
+16 impropers
+
+Types
+
+1 9
+2 1
+3 1
+4 4
+5 4
+6 3
+7 3
+8 1
+9 1
+10 5
+11 8
+12 6
+13 3
+14 3
+15 7
+
+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 9 1 2
+2 10 1 4
+3 10 1 5
+4 11 1 10
+5 1 2 3
+6 2 2 6
+7 2 2 7
+8 1 8 9
+9 6 9 10
+10 2 9 13
+11 2 9 14
+12 7 10 11
+13 5 10 12
+14 8 12 15
+
+Angles
+
+1 14 2 1 4
+2 14 2 1 5
+3 15 2 1 10
+4 16 4 1 5
+5 17 4 1 10
+6 17 5 1 10
+7 18 1 2 3
+8 19 1 2 6
+9 19 1 2 7
+10 1 3 2 6
+11 1 3 2 7
+12 3 6 2 7
+13 12 8 9 10
+14 1 8 9 13
+15 1 8 9 14
+16 13 13 9 10
+17 13 14 9 10
+18 3 13 9 14
+19 10 9 10 11
+20 8 9 10 12
+21 20 1 10 9
+22 21 11 10 12
+23 22 1 10 11
+24 23 1 10 12
+25 11 10 12 15
+
+Dihedrals
+
+1 16 4 1 2 3
+2 17 4 1 2 6
+3 17 4 1 2 7
+4 16 5 1 2 3
+5 17 5 1 2 6
+6 17 5 1 2 7
+7 18 10 1 2 3
+8 19 10 1 2 6
+9 19 10 1 2 7
+10 20 2 1 10 9
+11 21 2 1 10 11
+12 22 2 1 10 12
+13 23 4 1 10 9
+14 24 4 1 10 11
+15 25 4 1 10 12
+16 23 5 1 10 9
+17 24 5 1 10 11
+18 25 5 1 10 12
+19 14 8 9 10 11
+20 12 8 9 10 12
+21 28 8 9 10 1
+22 15 13 9 10 11
+23 13 13 9 10 12
+24 29 13 9 10 1
+25 15 14 9 10 11
+26 13 14 9 10 12
+27 29 14 9 10 1
+28 10 9 10 12 15
+29 11 11 10 12 15
+30 30 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/USER/misc/bond_react/tiny_nylon/in.tiny_nylon b/examples/USER/misc/bond_react/tiny_nylon/in.tiny_nylon
new file mode 100644
index 0000000000000000000000000000000000000000..3e21f69331c46b9527f224d58b1d03ca941705a3
--- /dev/null
+++ b/examples/USER/misc/bond_react/tiny_nylon/in.tiny_nylon
@@ -0,0 +1,50 @@
+# two monomer nylon example
+# reaction produces a condensed water molecule
+
+units real
+
+boundary p p p
+
+atom_style full
+
+kspace_style pppm 1.0e-4
+
+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/USER/misc/bond_react/tiny_nylon/log.20Apr18.tiny_nylon.g++.1 b/examples/USER/misc/bond_react/tiny_nylon/log.20Apr18.tiny_nylon.g++.1
new file mode 100644
index 0000000000000000000000000000000000000000..a65a5306f3d8fd89ad5d6b2c4fcb757a6a90a395
--- /dev/null
+++ b/examples/USER/misc/bond_react/tiny_nylon/log.20Apr18.tiny_nylon.g++.1
@@ -0,0 +1,365 @@
+LAMMPS (20 Apr 2018)
+OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:90)
+ using 1 OpenMP thread(s) per MPI task
+# two monomer nylon example
+# reaction produces a condensed water molecule
+
+units real
+
+boundary p p p
+
+atom_style full
+
+kspace_style pppm 1.0e-4
+
+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
+ orthogonal box = (-25 -25 -25) to (25 25 25)
+ 1 by 1 by 1 MPI processor grid
+ reading atoms ...
+ 44 atoms
+ reading velocities ...
+ 44 velocities
+ scanning bonds ...
+ 9 = max bonds/atom
+ scanning angles ...
+ 21 = max angles/atom
+ scanning dihedrals ...
+ 29 = max dihedrals/atom
+ scanning impropers ...
+ 29 = max impropers/atom
+ reading bonds ...
+ 42 bonds
+ reading angles ...
+ 74 angles
+ reading dihedrals ...
+ 100 dihedrals
+ reading impropers ...
+ 44 impropers
+ 4 = max # of 1-2 neighbors
+ 6 = max # of 1-3 neighbors
+ 12 = max # of 1-4 neighbors
+ 41 = max # of special neighbors
+
+velocity all create 300.0 4928459 dist gaussian
+
+molecule mol1 rxn1_stp1_unreacted.data_template
+Read molecule mol1:
+ 18 atoms with max type 8
+ 16 bonds with max type 14
+ 25 angles with max type 28
+ 23 dihedrals with max type 36
+ 14 impropers with max type 11
+molecule mol2 rxn1_stp1_reacted.data_template
+Read molecule mol2:
+ 18 atoms with max type 9
+ 17 bonds with max type 13
+ 31 angles with max type 27
+ 39 dihedrals with max type 33
+ 20 impropers with max type 1
+molecule mol3 rxn1_stp2_unreacted.data_template
+Read molecule mol3:
+ 15 atoms with max type 9
+ 14 bonds with max type 13
+ 25 angles with max type 27
+ 30 dihedrals with max type 33
+ 16 impropers with max type 1
+molecule mol4 rxn1_stp2_reacted.data_template
+Read molecule mol4:
+ 15 atoms with max type 11
+ 13 bonds with max type 15
+ 19 angles with max type 29
+ 16 dihedrals with max type 32
+ 10 impropers with max type 13
+
+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
+dynamic group bond_react_MASTER_group defined
+dynamic group statted_grp defined
+dynamic group bond_react_MASTER_group defined
+dynamic group statted_grp defined
+
+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
+PPPM initialization ...
+ using 12-bit tables for long-range coulomb (../kspace.cpp:321)
+ G vector (1/distance) = 0.0534597
+ grid = 2 2 2
+ stencil order = 5
+ estimated absolute RMS force accuracy = 0.0402256
+ estimated relative force accuracy = 0.000121138
+ using double precision FFTs
+ 3d grid and FFT values/proc = 343 8
+Neighbor list info ...
+ update every 1 steps, delay 10 steps, check yes
+ max neighbors/atom: 2000, page size: 100000
+ master list distance cutoff = 10.5
+ ghost atom cutoff = 10.5
+ binsize = 5.25, bins = 10 10 10
+ 2 neighbor lists, perpetual/occasional/extra = 1 1 0
+ (1) pair lj/class2/coul/long, perpetual
+ attributes: half, newton on
+ pair build: half/bin/newton
+ stencil: half/bin/3d/newton
+ bin: standard
+ (2) fix bond/react, occasional, copy from (1)
+ attributes: half, newton on
+ pair build: copy
+ stencil: none
+ bin: none
+WARNING: Inconsistent image flags (../domain.cpp:786)
+Per MPI rank memory allocation (min/avg/max) = 33.34 | 33.34 | 33.34 Mbytes
+Step Temp Press Density f_myrxns[1] f_myrxns[2]
+ 0 300 346.78165 0.0034851739 0 0
+ 50 296.70408 -51.30066 0.0034851739 1 0
+ 100 275.26011 39.120329 0.0034851739 1 1
+ 150 438.68516 35.257539 0.0034851739 1 1
+ 200 394.48971 15.444537 0.0034851739 1 1
+ 250 356.00369 50.185792 0.0034851739 1 1
+ 300 301.25816 -26.891497 0.0034851739 1 1
+ 350 279.17264 12.694513 0.0034851739 1 1
+ 400 248.71641 10.684558 0.0034851739 1 1
+ 450 240.30602 12.963034 0.0034851739 1 1
+ 500 252.71279 0.91620356 0.0034851739 1 1
+ 550 300.56824 18.541436 0.0034851739 1 1
+ 600 306.46441 -1.9736081 0.0034851739 1 1
+ 650 336.4537 21.474831 0.0034851739 1 1
+ 700 323.46217 8.486376 0.0034851739 1 1
+ 750 271.87146 5.9615231 0.0034851739 1 1
+ 800 268.43041 36.676068 0.0034851739 1 1
+ 850 269.02683 7.6295416 0.0034851739 1 1
+ 900 282.03605 -1.4688833 0.0034851739 1 1
+ 950 243.69136 -22.771489 0.0034851739 1 1
+ 1000 285.01348 17.925748 0.0034851739 1 1
+ 1050 383.47985 0.70536985 0.0034851739 1 1
+ 1100 368.97167 -4.3046933 0.0034851739 1 1
+ 1150 373.64459 6.3491837 0.0034851739 1 1
+ 1200 332.90575 -36.501095 0.0034851739 1 1
+ 1250 314.04078 -13.382767 0.0034851739 1 1
+ 1300 305.84166 4.7713641 0.0034851739 1 1
+ 1350 286.22145 37.621803 0.0034851739 1 1
+ 1400 265.52855 23.621002 0.0034851739 1 1
+ 1450 281.807 -31.266828 0.0034851739 1 1
+ 1500 275.33855 33.495565 0.0034851739 1 1
+ 1550 273.04973 -22.913871 0.0034851739 1 1
+ 1600 308.18478 -17.98151 0.0034851739 1 1
+ 1650 333.27664 28.987529 0.0034851739 1 1
+ 1700 296.16091 -1.1440455 0.0034851739 1 1
+ 1750 365.75611 34.574845 0.0034851739 1 1
+ 1800 292.83503 -33.199654 0.0034851739 1 1
+ 1850 261.50282 42.608703 0.0034851739 1 1
+ 1900 315.14188 -31.058803 0.0034851739 1 1
+ 1950 335.12895 12.40597 0.0034851739 1 1
+ 2000 278.08586 -5.3178633 0.0034851739 1 1
+ 2050 283.23847 -3.3974302 0.0034851739 1 1
+ 2100 281.38904 0.70263324 0.0034851739 1 1
+ 2150 302.23197 23.372316 0.0034851739 1 1
+ 2200 337.49259 -4.8716803 0.0034851739 1 1
+ 2250 409.64134 -12.043721 0.0034851739 1 1
+ 2300 309.21764 -21.824645 0.0034851739 1 1
+ 2350 290.97879 18.690281 0.0034851739 1 1
+ 2400 341.816 -16.967741 0.0034851739 1 1
+ 2450 310.28034 27.375518 0.0034851739 1 1
+ 2500 248.89429 17.061192 0.0034851739 1 1
+ 2550 273.10927 0.10481657 0.0034851739 1 1
+ 2600 289.56101 -9.7239939 0.0034851739 1 1
+ 2650 298.99719 -33.140576 0.0034851739 1 1
+ 2700 326.05198 -24.509827 0.0034851739 1 1
+ 2750 319.20612 24.305526 0.0034851739 1 1
+ 2800 304.8715 -15.076941 0.0034851739 1 1
+ 2850 374.38923 2.0874883 0.0034851739 1 1
+ 2900 354.01554 -20.595102 0.0034851739 1 1
+ 2950 289.89296 48.39731 0.0034851739 1 1
+ 3000 312.10013 -8.5105997 0.0034851739 1 1
+ 3050 296.97004 -31.008446 0.0034851739 1 1
+ 3100 251.72228 35.710197 0.0034851739 1 1
+ 3150 315.35895 -43.331536 0.0034851739 1 1
+ 3200 334.67773 13.331058 0.0034851739 1 1
+ 3250 308.1669 37.236121 0.0034851739 1 1
+ 3300 329.47601 30.798598 0.0034851739 1 1
+ 3350 299.40055 2.0785585 0.0034851739 1 1
+ 3400 272.41031 32.744922 0.0034851739 1 1
+ 3450 279.34594 -26.181793 0.0034851739 1 1
+ 3500 288.89969 8.935052 0.0034851739 1 1
+ 3550 253.4967 9.7244709 0.0034851739 1 1
+ 3600 294.83266 19.33305 0.0034851739 1 1
+ 3650 290.23794 -5.4939069 0.0034851739 1 1
+ 3700 332.5222 -29.834229 0.0034851739 1 1
+ 3750 364.63024 20.706191 0.0034851739 1 1
+ 3800 295.3842 -6.9434004 0.0034851739 1 1
+ 3850 346.84424 37.796066 0.0034851739 1 1
+ 3900 265.67286 -0.31628068 0.0034851739 1 1
+ 3950 260.455 -2.2571902 0.0034851739 1 1
+ 4000 259.82636 -2.2286205 0.0034851739 1 1
+ 4050 257.79848 24.520293 0.0034851739 1 1
+ 4100 295.58626 -0.42318936 0.0034851739 1 1
+ 4150 265.81353 -49.092436 0.0034851739 1 1
+ 4200 302.10333 51.715259 0.0034851739 1 1
+ 4250 258.98448 -4.8516655 0.0034851739 1 1
+ 4300 327.83401 33.717282 0.0034851739 1 1
+ 4350 311.59571 23.580382 0.0034851739 1 1
+ 4400 300.64237 -31.866661 0.0034851739 1 1
+ 4450 294.15643 -0.11927421 0.0034851739 1 1
+ 4500 299.83605 -17.560873 0.0034851739 1 1
+ 4550 326.83265 32.818482 0.0034851739 1 1
+ 4600 260.39068 -8.0567902 0.0034851739 1 1
+ 4650 247.93553 19.462991 0.0034851739 1 1
+ 4700 214.22252 -34.118303 0.0034851739 1 1
+ 4750 203.15329 27.356205 0.0034851739 1 1
+ 4800 257.761 -10.407989 0.0034851739 1 1
+ 4850 307.1923 11.71101 0.0034851739 1 1
+ 4900 319.00942 4.7808306 0.0034851739 1 1
+ 4950 282.23989 24.996151 0.0034851739 1 1
+ 5000 311.53284 -3.0012665 0.0034851739 1 1
+ 5050 317.58212 32.567832 0.0034851739 1 1
+ 5100 267.51501 -47.695087 0.0034851739 1 1
+ 5150 260.19048 29.046388 0.0034851739 1 1
+ 5200 239.83552 -5.4890385 0.0034851739 1 1
+ 5250 234.8852 -18.172633 0.0034851739 1 1
+ 5300 236.43277 -39.06212 0.0034851739 1 1
+ 5350 280.90079 -2.6932923 0.0034851739 1 1
+ 5400 316.65266 23.071362 0.0034851739 1 1
+ 5450 345.63226 19.573323 0.0034851739 1 1
+ 5500 384.57334 41.507217 0.0034851739 1 1
+ 5550 317.14278 9.6992897 0.0034851739 1 1
+ 5600 279.93243 -12.076895 0.0034851739 1 1
+ 5650 268.06471 1.6196401 0.0034851739 1 1
+ 5700 271.85714 -40.378455 0.0034851739 1 1
+ 5750 313.88363 10.722639 0.0034851739 1 1
+ 5800 281.54495 31.914889 0.0034851739 1 1
+ 5850 293.34821 -8.3154922 0.0034851739 1 1
+ 5900 249.25216 -17.307353 0.0034851739 1 1
+ 5950 268.18639 -4.7222512 0.0034851739 1 1
+ 6000 302.99398 -52.615528 0.0034851739 1 1
+ 6050 314.57931 34.51318 0.0034851739 1 1
+ 6100 345.70348 30.334721 0.0034851739 1 1
+ 6150 316.59329 31.862519 0.0034851739 1 1
+ 6200 317.85346 -32.235221 0.0034851739 1 1
+ 6250 282.97676 0.2936745 0.0034851739 1 1
+ 6300 267.91814 19.265567 0.0034851739 1 1
+ 6350 226.20967 -13.093547 0.0034851739 1 1
+ 6400 307.73316 17.439598 0.0034851739 1 1
+ 6450 292.16253 -23.275163 0.0034851739 1 1
+ 6500 335.05939 26.936463 0.0034851739 1 1
+ 6550 380.73546 19.532416 0.0034851739 1 1
+ 6600 373.0103 30.879532 0.0034851739 1 1
+ 6650 335.37975 -2.1762828 0.0034851739 1 1
+ 6700 298.94272 -10.578587 0.0034851739 1 1
+ 6750 255.11531 -50.576215 0.0034851739 1 1
+ 6800 222.87459 3.0499548 0.0034851739 1 1
+ 6850 268.57213 -43.675945 0.0034851739 1 1
+ 6900 260.3024 4.7483005 0.0034851739 1 1
+ 6950 289.15855 31.62106 0.0034851739 1 1
+ 7000 289.11874 21.635533 0.0034851739 1 1
+ 7050 361.08776 22.445996 0.0034851739 1 1
+ 7100 368.95003 4.8383881 0.0034851739 1 1
+ 7150 331.47448 -36.200495 0.0034851739 1 1
+ 7200 304.7251 13.982693 0.0034851739 1 1
+ 7250 284.09747 0.53758275 0.0034851739 1 1
+ 7300 269.17023 -41.571482 0.0034851739 1 1
+ 7350 222.07071 25.564662 0.0034851739 1 1
+ 7400 304.09598 15.482955 0.0034851739 1 1
+ 7450 298.78752 -7.4335841 0.0034851739 1 1
+ 7500 328.78697 14.666097 0.0034851739 1 1
+ 7550 347.07038 -37.165295 0.0034851739 1 1
+ 7600 362.85673 20.52268 0.0034851739 1 1
+ 7650 347.15141 2.3383775 0.0034851739 1 1
+ 7700 262.10132 33.898374 0.0034851739 1 1
+ 7750 275.84724 -33.534813 0.0034851739 1 1
+ 7800 281.14075 -18.284372 0.0034851739 1 1
+ 7850 264.83337 -30.580199 0.0034851739 1 1
+ 7900 257.35275 -35.080567 0.0034851739 1 1
+ 7950 286.32446 26.594779 0.0034851739 1 1
+ 8000 248.36889 15.605894 0.0034851739 1 1
+ 8050 292.55015 28.811985 0.0034851739 1 1
+ 8100 312.47888 0.83990451 0.0034851739 1 1
+ 8150 285.58532 -15.258185 0.0034851739 1 1
+ 8200 292.22819 -38.233195 0.0034851739 1 1
+ 8250 321.6208 -19.052143 0.0034851739 1 1
+ 8300 319.41332 54.97437 0.0034851739 1 1
+ 8350 307.95647 32.009591 0.0034851739 1 1
+ 8400 345.58105 8.8535539 0.0034851739 1 1
+ 8450 357.75168 12.416896 0.0034851739 1 1
+ 8500 370.049 4.3288665 0.0034851739 1 1
+ 8550 360.62882 12.618614 0.0034851739 1 1
+ 8600 290.10834 -4.8081765 0.0034851739 1 1
+ 8650 297.7575 -5.1976816 0.0034851739 1 1
+ 8700 286.57505 -31.469549 0.0034851739 1 1
+ 8750 307.77059 19.338001 0.0034851739 1 1
+ 8800 231.68316 12.159459 0.0034851739 1 1
+ 8850 329.13623 -8.7262592 0.0034851739 1 1
+ 8900 286.40715 10.326025 0.0034851739 1 1
+ 8950 339.43101 2.7809618 0.0034851739 1 1
+ 9000 402.53799 19.481869 0.0034851739 1 1
+ 9050 349.56449 -4.8450179 0.0034851739 1 1
+ 9100 307.64739 16.889327 0.0034851739 1 1
+ 9150 276.54451 -34.808372 0.0034851739 1 1
+ 9200 233.18668 4.9409791 0.0034851739 1 1
+ 9250 266.48384 -19.850366 0.0034851739 1 1
+ 9300 289.14808 13.520201 0.0034851739 1 1
+ 9350 295.08335 17.156468 0.0034851739 1 1
+ 9400 338.08757 -31.112278 0.0034851739 1 1
+ 9450 336.64739 -25.697747 0.0034851739 1 1
+ 9500 338.10622 1.9241797 0.0034851739 1 1
+ 9550 294.82158 -12.043972 0.0034851739 1 1
+ 9600 268.9836 12.235553 0.0034851739 1 1
+ 9650 279.6269 28.710734 0.0034851739 1 1
+ 9700 279.88562 -10.865604 0.0034851739 1 1
+ 9750 287.56565 12.975819 0.0034851739 1 1
+ 9800 278.39949 4.2088595 0.0034851739 1 1
+ 9850 307.61259 9.341169 0.0034851739 1 1
+ 9900 317.53581 2.3948493 0.0034851739 1 1
+ 9950 332.52938 -14.809185 0.0034851739 1 1
+ 10000 401.93365 -7.637581 0.0034851739 1 1
+Loop time of 1.94139 on 1 procs for 10000 steps with 44 atoms
+
+Performance: 445.042 ns/day, 0.054 hours/ns, 5150.945 timesteps/s
+99.3% CPU use with 1 MPI tasks x 1 OpenMP threads
+
+MPI task timing breakdown:
+Section | min time | avg time | max time |%varavg| %total
+---------------------------------------------------------------
+Pair | 0.26479 | 0.26479 | 0.26479 | 0.0 | 13.64
+Bond | 0.76875 | 0.76875 | 0.76875 | 0.0 | 39.60
+Kspace | 0.32111 | 0.32111 | 0.32111 | 0.0 | 16.54
+Neigh | 0.41333 | 0.41333 | 0.41333 | 0.0 | 21.29
+Comm | 0.025956 | 0.025956 | 0.025956 | 0.0 | 1.34
+Output | 0.0043445 | 0.0043445 | 0.0043445 | 0.0 | 0.22
+Modify | 0.12526 | 0.12526 | 0.12526 | 0.0 | 6.45
+Other | | 0.01786 | | | 0.92
+
+Nlocal: 44 ave 44 max 44 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Nghost: 62 ave 62 max 62 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Neighs: 812 ave 812 max 812 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+
+Total # of neighbors = 812
+Ave neighs/atom = 18.4545
+Ave special neighs/atom = 9.77273
+Neighbor list builds = 10000
+Dangerous builds = 0
+
+# write_restart restart_longrun
+# write_data restart_longrun.data
+Total wall time: 0:00:01
diff --git a/examples/USER/misc/bond_react/tiny_nylon/log.20Apr18.tiny_nylon.g++.4 b/examples/USER/misc/bond_react/tiny_nylon/log.20Apr18.tiny_nylon.g++.4
new file mode 100644
index 0000000000000000000000000000000000000000..dddc9f2801e6c62482640b4bed730d57e324854b
--- /dev/null
+++ b/examples/USER/misc/bond_react/tiny_nylon/log.20Apr18.tiny_nylon.g++.4
@@ -0,0 +1,365 @@
+LAMMPS (20 Apr 2018)
+OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:90)
+ using 1 OpenMP thread(s) per MPI task
+# two monomer nylon example
+# reaction produces a condensed water molecule
+
+units real
+
+boundary p p p
+
+atom_style full
+
+kspace_style pppm 1.0e-4
+
+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
+ orthogonal box = (-25 -25 -25) to (25 25 25)
+ 1 by 2 by 2 MPI processor grid
+ reading atoms ...
+ 44 atoms
+ reading velocities ...
+ 44 velocities
+ scanning bonds ...
+ 9 = max bonds/atom
+ scanning angles ...
+ 21 = max angles/atom
+ scanning dihedrals ...
+ 29 = max dihedrals/atom
+ scanning impropers ...
+ 29 = max impropers/atom
+ reading bonds ...
+ 42 bonds
+ reading angles ...
+ 74 angles
+ reading dihedrals ...
+ 100 dihedrals
+ reading impropers ...
+ 44 impropers
+ 4 = max # of 1-2 neighbors
+ 6 = max # of 1-3 neighbors
+ 12 = max # of 1-4 neighbors
+ 41 = max # of special neighbors
+
+velocity all create 300.0 4928459 dist gaussian
+
+molecule mol1 rxn1_stp1_unreacted.data_template
+Read molecule mol1:
+ 18 atoms with max type 8
+ 16 bonds with max type 14
+ 25 angles with max type 28
+ 23 dihedrals with max type 36
+ 14 impropers with max type 11
+molecule mol2 rxn1_stp1_reacted.data_template
+Read molecule mol2:
+ 18 atoms with max type 9
+ 17 bonds with max type 13
+ 31 angles with max type 27
+ 39 dihedrals with max type 33
+ 20 impropers with max type 1
+molecule mol3 rxn1_stp2_unreacted.data_template
+Read molecule mol3:
+ 15 atoms with max type 9
+ 14 bonds with max type 13
+ 25 angles with max type 27
+ 30 dihedrals with max type 33
+ 16 impropers with max type 1
+molecule mol4 rxn1_stp2_reacted.data_template
+Read molecule mol4:
+ 15 atoms with max type 11
+ 13 bonds with max type 15
+ 19 angles with max type 29
+ 16 dihedrals with max type 32
+ 10 impropers with max type 13
+
+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
+dynamic group bond_react_MASTER_group defined
+dynamic group statted_grp defined
+dynamic group bond_react_MASTER_group defined
+dynamic group statted_grp defined
+
+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
+PPPM initialization ...
+ using 12-bit tables for long-range coulomb (../kspace.cpp:321)
+ G vector (1/distance) = 0.0534597
+ grid = 2 2 2
+ stencil order = 5
+ estimated absolute RMS force accuracy = 0.0402256
+ estimated relative force accuracy = 0.000121138
+ using double precision FFTs
+ 3d grid and FFT values/proc = 252 2
+Neighbor list info ...
+ update every 1 steps, delay 10 steps, check yes
+ max neighbors/atom: 2000, page size: 100000
+ master list distance cutoff = 10.5
+ ghost atom cutoff = 10.5
+ binsize = 5.25, bins = 10 10 10
+ 2 neighbor lists, perpetual/occasional/extra = 1 1 0
+ (1) pair lj/class2/coul/long, perpetual
+ attributes: half, newton on
+ pair build: half/bin/newton
+ stencil: half/bin/3d/newton
+ bin: standard
+ (2) fix bond/react, occasional, copy from (1)
+ attributes: half, newton on
+ pair build: copy
+ stencil: none
+ bin: none
+WARNING: Inconsistent image flags (../domain.cpp:786)
+Per MPI rank memory allocation (min/avg/max) = 33.34 | 33.69 | 34.37 Mbytes
+Step Temp Press Density f_myrxns[1] f_myrxns[2]
+ 0 300 346.78165 0.0034851739 0 0
+ 50 296.70408 -51.30066 0.0034851739 1 0
+ 100 275.26011 39.120329 0.0034851739 1 1
+ 150 438.68516 35.257539 0.0034851739 1 1
+ 200 394.48971 15.444537 0.0034851739 1 1
+ 250 356.00369 50.185792 0.0034851739 1 1
+ 300 301.25816 -26.891497 0.0034851739 1 1
+ 350 279.17264 12.694513 0.0034851739 1 1
+ 400 248.71641 10.684558 0.0034851739 1 1
+ 450 240.30602 12.963034 0.0034851739 1 1
+ 500 252.71279 0.91620356 0.0034851739 1 1
+ 550 300.56824 18.541436 0.0034851739 1 1
+ 600 306.46441 -1.9736081 0.0034851739 1 1
+ 650 336.4537 21.474831 0.0034851739 1 1
+ 700 323.46217 8.486376 0.0034851739 1 1
+ 750 271.87146 5.9615231 0.0034851739 1 1
+ 800 268.43041 36.676068 0.0034851739 1 1
+ 850 269.02683 7.6295416 0.0034851739 1 1
+ 900 282.03605 -1.4688833 0.0034851739 1 1
+ 950 243.69136 -22.771489 0.0034851739 1 1
+ 1000 285.01348 17.925748 0.0034851739 1 1
+ 1050 383.47985 0.70536985 0.0034851739 1 1
+ 1100 368.97167 -4.3046933 0.0034851739 1 1
+ 1150 373.64459 6.3491837 0.0034851739 1 1
+ 1200 332.90575 -36.501095 0.0034851739 1 1
+ 1250 314.04078 -13.382767 0.0034851739 1 1
+ 1300 305.84166 4.7713641 0.0034851739 1 1
+ 1350 286.22145 37.621803 0.0034851739 1 1
+ 1400 265.52855 23.621002 0.0034851739 1 1
+ 1450 281.807 -31.266828 0.0034851739 1 1
+ 1500 275.33855 33.495565 0.0034851739 1 1
+ 1550 273.04973 -22.913871 0.0034851739 1 1
+ 1600 308.18478 -17.98151 0.0034851739 1 1
+ 1650 333.27664 28.987529 0.0034851739 1 1
+ 1700 296.16091 -1.1440455 0.0034851739 1 1
+ 1750 365.75611 34.574845 0.0034851739 1 1
+ 1800 292.83503 -33.199654 0.0034851739 1 1
+ 1850 261.50282 42.608703 0.0034851739 1 1
+ 1900 315.14188 -31.058803 0.0034851739 1 1
+ 1950 335.12895 12.40597 0.0034851739 1 1
+ 2000 278.08586 -5.3178633 0.0034851739 1 1
+ 2050 283.23847 -3.3974302 0.0034851739 1 1
+ 2100 281.38904 0.70263324 0.0034851739 1 1
+ 2150 302.23197 23.372316 0.0034851739 1 1
+ 2200 337.49259 -4.8716803 0.0034851739 1 1
+ 2250 409.64134 -12.043721 0.0034851739 1 1
+ 2300 309.21764 -21.824645 0.0034851739 1 1
+ 2350 290.97879 18.690281 0.0034851739 1 1
+ 2400 341.816 -16.967741 0.0034851739 1 1
+ 2450 310.28034 27.375518 0.0034851739 1 1
+ 2500 248.89429 17.061192 0.0034851739 1 1
+ 2550 273.10927 0.1048166 0.0034851739 1 1
+ 2600 289.56101 -9.7239939 0.0034851739 1 1
+ 2650 298.99719 -33.140576 0.0034851739 1 1
+ 2700 326.05198 -24.509827 0.0034851739 1 1
+ 2750 319.20612 24.305526 0.0034851739 1 1
+ 2800 304.8715 -15.076941 0.0034851739 1 1
+ 2850 374.38923 2.0874883 0.0034851739 1 1
+ 2900 354.01554 -20.595102 0.0034851739 1 1
+ 2950 289.89296 48.39731 0.0034851739 1 1
+ 3000 312.10013 -8.5105996 0.0034851739 1 1
+ 3050 296.97004 -31.008446 0.0034851739 1 1
+ 3100 251.72228 35.710197 0.0034851739 1 1
+ 3150 315.35895 -43.331536 0.0034851739 1 1
+ 3200 334.67773 13.331058 0.0034851739 1 1
+ 3250 308.1669 37.236121 0.0034851739 1 1
+ 3300 329.47601 30.798598 0.0034851739 1 1
+ 3350 299.40055 2.0785585 0.0034851739 1 1
+ 3400 272.41031 32.744921 0.0034851739 1 1
+ 3450 279.34594 -26.181793 0.0034851739 1 1
+ 3500 288.89969 8.935052 0.0034851739 1 1
+ 3550 253.4967 9.724471 0.0034851739 1 1
+ 3600 294.83266 19.33305 0.0034851739 1 1
+ 3650 290.23794 -5.493907 0.0034851739 1 1
+ 3700 332.5222 -29.83423 0.0034851739 1 1
+ 3750 364.63024 20.706191 0.0034851739 1 1
+ 3800 295.3842 -6.9434003 0.0034851739 1 1
+ 3850 346.84424 37.796066 0.0034851739 1 1
+ 3900 265.67286 -0.3162804 0.0034851739 1 1
+ 3950 260.455 -2.2571901 0.0034851739 1 1
+ 4000 259.82636 -2.2286207 0.0034851739 1 1
+ 4050 257.79848 24.520293 0.0034851739 1 1
+ 4100 295.58626 -0.42318895 0.0034851739 1 1
+ 4150 265.81352 -49.092436 0.0034851739 1 1
+ 4200 302.10333 51.715258 0.0034851739 1 1
+ 4250 258.98448 -4.8516657 0.0034851739 1 1
+ 4300 327.83401 33.717283 0.0034851739 1 1
+ 4350 311.59571 23.580382 0.0034851739 1 1
+ 4400 300.64237 -31.866661 0.0034851739 1 1
+ 4450 294.15642 -0.11927262 0.0034851739 1 1
+ 4500 299.83605 -17.560872 0.0034851739 1 1
+ 4550 326.83265 32.818481 0.0034851739 1 1
+ 4600 260.39068 -8.0567907 0.0034851739 1 1
+ 4650 247.93553 19.462991 0.0034851739 1 1
+ 4700 214.22252 -34.118304 0.0034851739 1 1
+ 4750 203.15329 27.356204 0.0034851739 1 1
+ 4800 257.761 -10.407986 0.0034851739 1 1
+ 4850 307.1923 11.711008 0.0034851739 1 1
+ 4900 319.00942 4.7808342 0.0034851739 1 1
+ 4950 282.2399 24.996151 0.0034851739 1 1
+ 5000 311.53284 -3.0012669 0.0034851739 1 1
+ 5050 317.58213 32.56782 0.0034851739 1 1
+ 5100 267.51502 -47.695103 0.0034851739 1 1
+ 5150 260.19047 29.046394 0.0034851739 1 1
+ 5200 239.8355 -5.4890372 0.0034851739 1 1
+ 5250 234.88522 -18.172649 0.0034851739 1 1
+ 5300 236.43278 -39.062111 0.0034851739 1 1
+ 5350 280.90083 -2.6932604 0.0034851739 1 1
+ 5400 316.65269 23.071363 0.0034851739 1 1
+ 5450 345.6322 19.573305 0.0034851739 1 1
+ 5500 384.57334 41.50729 0.0034851739 1 1
+ 5550 317.14286 9.6992981 0.0034851739 1 1
+ 5600 279.93246 -12.076859 0.0034851739 1 1
+ 5650 268.06471 1.6196502 0.0034851739 1 1
+ 5700 271.85714 -40.378489 0.0034851739 1 1
+ 5750 313.88361 10.722652 0.0034851739 1 1
+ 5800 281.54499 31.914883 0.0034851739 1 1
+ 5850 293.34819 -8.3155887 0.0034851739 1 1
+ 5900 249.25215 -17.307228 0.0034851739 1 1
+ 5950 268.18645 -4.7223601 0.0034851739 1 1
+ 6000 302.99402 -52.615432 0.0034851739 1 1
+ 6050 314.57946 34.513152 0.0034851739 1 1
+ 6100 345.70342 30.33474 0.0034851739 1 1
+ 6150 316.59329 31.862566 0.0034851739 1 1
+ 6200 317.85341 -32.23511 0.0034851739 1 1
+ 6250 282.97674 0.29367434 0.0034851739 1 1
+ 6300 267.91823 19.265617 0.0034851739 1 1
+ 6350 226.2098 -13.093573 0.0034851739 1 1
+ 6400 307.73307 17.439662 0.0034851739 1 1
+ 6450 292.16311 -23.275251 0.0034851739 1 1
+ 6500 335.05972 26.936588 0.0034851739 1 1
+ 6550 380.7351 19.532324 0.0034851739 1 1
+ 6600 373.01041 30.879146 0.0034851739 1 1
+ 6650 335.37897 -2.1766711 0.0034851739 1 1
+ 6700 298.94275 -10.578361 0.0034851739 1 1
+ 6750 255.11449 -50.575851 0.0034851739 1 1
+ 6800 222.87598 3.0488985 0.0034851739 1 1
+ 6850 268.57268 -43.676136 0.0034851739 1 1
+ 6900 260.30442 4.7484508 0.0034851739 1 1
+ 6950 289.15739 31.622589 0.0034851739 1 1
+ 7000 289.11733 21.636361 0.0034851739 1 1
+ 7050 361.08905 22.442487 0.0034851739 1 1
+ 7100 368.95006 4.8393179 0.0034851739 1 1
+ 7150 331.47878 -36.202032 0.0034851739 1 1
+ 7200 304.72518 13.982604 0.0034851739 1 1
+ 7250 284.0996 0.53900966 0.0034851739 1 1
+ 7300 269.17156 -41.572215 0.0034851739 1 1
+ 7350 222.06563 25.56579 0.0034851739 1 1
+ 7400 304.09479 15.48238 0.0034851739 1 1
+ 7450 298.79046 -7.4369454 0.0034851739 1 1
+ 7500 328.78217 14.672853 0.0034851739 1 1
+ 7550 347.06589 -37.168123 0.0034851739 1 1
+ 7600 362.84157 20.514912 0.0034851739 1 1
+ 7650 347.15916 2.3477485 0.0034851739 1 1
+ 7700 262.09822 33.901831 0.0034851739 1 1
+ 7750 275.85921 -33.536059 0.0034851739 1 1
+ 7800 281.16159 -18.288414 0.0034851739 1 1
+ 7850 264.83553 -30.566284 0.0034851739 1 1
+ 7900 257.35224 -35.087067 0.0034851739 1 1
+ 7950 286.30756 26.586163 0.0034851739 1 1
+ 8000 248.38175 15.601961 0.0034851739 1 1
+ 8050 292.59171 28.784541 0.0034851739 1 1
+ 8100 312.52852 0.87995053 0.0034851739 1 1
+ 8150 285.62346 -15.337252 0.0034851739 1 1
+ 8200 292.24175 -38.192576 0.0034851739 1 1
+ 8250 321.61618 -19.030398 0.0034851739 1 1
+ 8300 319.42189 55.078305 0.0034851739 1 1
+ 8350 308.00357 32.050518 0.0034851739 1 1
+ 8400 345.68186 8.7983733 0.0034851739 1 1
+ 8450 358.00849 12.434592 0.0034851739 1 1
+ 8500 370.14359 4.2184721 0.0034851739 1 1
+ 8550 360.6511 12.580836 0.0034851739 1 1
+ 8600 290.04938 -4.8422589 0.0034851739 1 1
+ 8650 297.95745 -5.3185591 0.0034851739 1 1
+ 8700 286.54284 -31.490479 0.0034851739 1 1
+ 8750 308.08791 19.24417 0.0034851739 1 1
+ 8800 231.72534 12.262217 0.0034851739 1 1
+ 8850 329.2349 -8.9133933 0.0034851739 1 1
+ 8900 287.64023 10.525164 0.0034851739 1 1
+ 8950 341.08296 2.80127 0.0034851739 1 1
+ 9000 403.71266 17.88418 0.0034851739 1 1
+ 9050 348.95132 -4.5813611 0.0034851739 1 1
+ 9100 307.98322 16.707575 0.0034851739 1 1
+ 9150 276.75719 -35.563923 0.0034851739 1 1
+ 9200 230.15547 3.8091656 0.0034851739 1 1
+ 9250 264.64479 -20.97438 0.0034851739 1 1
+ 9300 285.70467 13.881735 0.0034851739 1 1
+ 9350 297.515 13.599319 0.0034851739 1 1
+ 9400 341.45203 -28.494544 0.0034851739 1 1
+ 9450 319.29703 -27.907344 0.0034851739 1 1
+ 9500 332.50473 4.9891138 0.0034851739 1 1
+ 9550 293.10499 -8.3244772 0.0034851739 1 1
+ 9600 255.31174 2.2568315 0.0034851739 1 1
+ 9650 250.01932 5.7005159 0.0034851739 1 1
+ 9700 278.37409 -8.5970424 0.0034851739 1 1
+ 9750 294.86737 17.686447 0.0034851739 1 1
+ 9800 277.07345 -2.0856886 0.0034851739 1 1
+ 9850 295.54707 2.8365471 0.0034851739 1 1
+ 9900 311.51074 29.885116 0.0034851739 1 1
+ 9950 296.01363 12.206068 0.0034851739 1 1
+ 10000 341.35187 -3.0045464 0.0034851739 1 1
+Loop time of 3.64332 on 4 procs for 10000 steps with 44 atoms
+
+Performance: 237.146 ns/day, 0.101 hours/ns, 2744.751 timesteps/s
+94.5% CPU use with 4 MPI tasks x 1 OpenMP threads
+
+MPI task timing breakdown:
+Section | min time | avg time | max time |%varavg| %total
+---------------------------------------------------------------
+Pair | 0.0029962 | 0.10426 | 0.34564 | 43.9 | 2.86
+Bond | 0.005929 | 0.29803 | 0.95305 | 71.2 | 8.18
+Kspace | 0.83969 | 1.5992 | 1.9344 | 35.6 | 43.89
+Neigh | 0.65468 | 0.66443 | 0.67431 | 0.9 | 18.24
+Comm | 0.1727 | 0.23754 | 0.2745 | 8.0 | 6.52
+Output | 0.0048738 | 0.010774 | 0.028434 | 9.8 | 0.30
+Modify | 0.62478 | 0.70376 | 0.779 | 6.5 | 19.32
+Other | | 0.02531 | | | 0.69
+
+Nlocal: 11 ave 40 max 0 min
+Histogram: 2 1 0 0 0 0 0 0 0 1
+Nghost: 36 ave 47 max 7 min
+Histogram: 1 0 0 0 0 0 0 0 0 3
+Neighs: 203 ave 809 max 0 min
+Histogram: 3 0 0 0 0 0 0 0 0 1
+
+Total # of neighbors = 812
+Ave neighs/atom = 18.4545
+Ave special neighs/atom = 9.77273
+Neighbor list builds = 10000
+Dangerous builds = 0
+
+# write_restart restart_longrun
+# write_data restart_longrun.data
+Total wall time: 0:00:03
diff --git a/examples/USER/misc/bond_react/tiny_nylon/rxn1_stp1_map b/examples/USER/misc/bond_react/tiny_nylon/rxn1_stp1_map
new file mode 100644
index 0000000000000000000000000000000000000000..44f7ad8137549eaa375046e114c59a7aa09c0c65
--- /dev/null
+++ b/examples/USER/misc/bond_react/tiny_nylon/rxn1_stp1_map
@@ -0,0 +1,35 @@
+this is a nominal 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/USER/misc/bond_react/tiny_nylon/rxn1_stp1_reacted.data_template b/examples/USER/misc/bond_react/tiny_nylon/rxn1_stp1_reacted.data_template
new file mode 100644
index 0000000000000000000000000000000000000000..d7256f43d2b3c659159ef23e3169dc6123d193ec
--- /dev/null
+++ b/examples/USER/misc/bond_react/tiny_nylon/rxn1_stp1_reacted.data_template
@@ -0,0 +1,189 @@
+this is a molecule template for: initial nylon crosslink, post-reacting
+
+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/USER/misc/bond_react/tiny_nylon/rxn1_stp1_unreacted.data_template b/examples/USER/misc/bond_react/tiny_nylon/rxn1_stp1_unreacted.data_template
new file mode 100644
index 0000000000000000000000000000000000000000..ec3f109d7b50c06c642635693f77257c4ce68c7e
--- /dev/null
+++ b/examples/USER/misc/bond_react/tiny_nylon/rxn1_stp1_unreacted.data_template
@@ -0,0 +1,160 @@
+this is a molecule template for: initial nylon crosslink, pre-reacting
+
+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/USER/misc/bond_react/tiny_nylon/rxn1_stp2_map b/examples/USER/misc/bond_react/tiny_nylon/rxn1_stp2_map
new file mode 100644
index 0000000000000000000000000000000000000000..35fe47fdb3bdb0386d65b950d45cbfb39fbe6f7f
--- /dev/null
+++ b/examples/USER/misc/bond_react/tiny_nylon/rxn1_stp2_map
@@ -0,0 +1,32 @@
+this is a nominal 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/USER/misc/bond_react/tiny_nylon/rxn1_stp2_reacted.data_template b/examples/USER/misc/bond_react/tiny_nylon/rxn1_stp2_reacted.data_template
new file mode 100644
index 0000000000000000000000000000000000000000..785363464696967b7938c142043868792546e3cf
--- /dev/null
+++ b/examples/USER/misc/bond_react/tiny_nylon/rxn1_stp2_reacted.data_template
@@ -0,0 +1,131 @@
+this is a molecule template for: water condensation, post-reacting
+
+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/USER/misc/bond_react/tiny_nylon/rxn1_stp2_unreacted.data_template b/examples/USER/misc/bond_react/tiny_nylon/rxn1_stp2_unreacted.data_template
new file mode 100644
index 0000000000000000000000000000000000000000..847f0622e506865117fbe542f89de35fb20ccede
--- /dev/null
+++ b/examples/USER/misc/bond_react/tiny_nylon/rxn1_stp2_unreacted.data_template
@@ -0,0 +1,158 @@
+this is a molecule template for: water condensation, pre-reacting
+
+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/USER/misc/bond_react/tiny_nylon/tiny_nylon.data b/examples/USER/misc/bond_react/tiny_nylon/tiny_nylon.data
new file mode 100644
index 0000000000000000000000000000000000000000..8466e68ea57fed2ddca9851724e4ea1509049435
--- /dev/null
+++ b/examples/USER/misc/bond_react/tiny_nylon/tiny_nylon.data
@@ -0,0 +1,795 @@
+this is LAMMPS data file containing two nylon monomers
+
+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/.gitignore b/src/.gitignore
index 94d55b87052fccae48f4d24b3db41695aa2f77a5..46c8f1ca080b4bb909a3653da1ad827cda6bcf91 100644
--- a/src/.gitignore
+++ b/src/.gitignore
@@ -374,6 +374,8 @@
/fix_bond_break.h
/fix_bond_create.cpp
/fix_bond_create.h
+/fix_bond_react.cpp
+/fix_bond_react.h
/fix_bond_swap.cpp
/fix_bond_swap.h
/fix_cmap.cpp
diff --git a/src/USER-MISC/README b/src/USER-MISC/README
index a8c33fa38056454194ff2346c9bd2d4822019200..7fc931b9625ed9fa6884bb5b1564692cf77564dd 100644
--- a/src/USER-MISC/README
+++ b/src/USER-MISC/README
@@ -39,6 +39,7 @@ dihedral_style spherical, Andrew Jewett, jewett.aij@gmail.com, 15 Jul 16
dihedral_style table, Andrew Jewett, jewett.aij@gmail.com, 10 Jan 12
fix addtorque, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11
fix ave/correlate/long, Jorge Ramirez (UPM Madrid), jorge.ramirez at upm.es, 21 Oct 2015
+fix bond/react, Jacob Gissinger (CU Boulder), info at disarmmd.org, 24 Feb 2018
fix filter/corotate, Lukas Fath (KIT), lukas.fath at kit.edu, 15 Mar 2017
fix flow/gauss, Joel Eaves (CU Boulder), Joel.Eaves@Colorado.edu, 23 Aug 2016
fix gle, Michele Ceriotti (EPFL Lausanne), michele.ceriotti at gmail.com, 24 Nov 2014
diff --git a/src/USER-MISC/fix_bond_react.cpp b/src/USER-MISC/fix_bond_react.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b3c92eed67354745665760ee5b0f8bb29b68b9ef
--- /dev/null
+++ b/src/USER-MISC/fix_bond_react.cpp
@@ -0,0 +1,2606 @@
+/* ----------------------------------------------------------------------
+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 <mpi.h>
+#include <cmath>
+#include <cstring>
+#include <cstdlib>
+#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 "molecule.h"
+#include "group.h"
+#include "citeme.h"
+#include "memory.h"
+#include "error.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)
+{
+ fix1 = NULL;
+ fix2 = 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";
+
+ // by using fixed group names, only one instance of fix bond/react is allowed.
+ if (modify->find_fix_by_style("bond/react") != -1)
+ error->all(FLERR,"Only one instance of fix bond/react allowed at a time");
+
+ // 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;
+}
+
+/* ---------------------------------------------------------------------- */
+
+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);
+ }
+
+ 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;
+
+ 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
+ int len = strlen("per_atom_props") + 1;
+ id_fix2 = new char[len];
+ strcpy(id_fix2,"per_atom_props");
+
+ int ifix = modify->find_fix(id_fix2);
+ if (ifix == -1) {
+ char **newarg = new char*[8];
+ newarg[0] = (char *) "per_atom_props";
+ newarg[1] = (char *) "all"; // group ID is ignored
+ newarg[2] = (char *) "property/atom";
+ newarg[3] = (char *) "i_limit_tags";
+ newarg[4] = (char *) "i_statted_tags";
+ newarg[5] = (char *) "i_react_tags";
+ newarg[6] = (char *) "ghost";
+ newarg[7] = (char *) "yes";
+ modify->add_fix(8,newarg);
+ fix2 = modify->fix[modify->nfix-1];
+ delete [] newarg;
+ }
+
+ // limit_tags: these are recently reacted atoms being relaxed
+ // per-atom properties already initialized to zero (not in group)
+ // let's do it anyway for clarity
+ int flag;
+ int index = atom->find_custom("limit_tags",flag); //here's where error would happen
+ 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;
+ }
+
+ // on to statted_tags (system-wide thermostat)
+ // intialize per-atom statted_flags to 1
+ index = atom->find_custom("statted_tags",flag);
+ int *i_statted_tags = atom->ivector[index];
+
+ for (int i = 0; i < atom->nlocal; i++)
+ i_statted_tags[i] = 1;
+
+ 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 = modify->fix[modify->nfix-1];
+ delete [] newarg;
+ }
+
+ }
+
+ //react_tags: 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
+
+ //per-atom values initalized to 0
+ index = atom->find_custom("react_tags",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,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;
+ int **nspecial = atom->nspecial;
+ tagint **special = atom->special;
+ int *mask = atom->mask;
+ int *type = atom->type;
+
+ neighbor->build_one(list,1);
+ inum = list->inum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ // per-atom property indicating if in bond/react master group
+ int flag;
+ int index1 = atom->find_custom("limit_tags",flag);
+ int *i_limit_tags = atom->ivector[index1];
+
+ int i,j;
+
+ for (int myrxn = 0; myrxn < nreacts; myrxn++) {
+
+ for (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 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
+
+ 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 nfirst_neighs = onemol->nspecial[pion][0];
+
+ // per-atom property indicating if in bond/react master group
+ int flag;
+ int index1 = atom->find_custom("limit_tags",flag);
+ int *i_limit_tags = atom->ivector[index1];
+
+ 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 nfirst_neighs = onemol->nspecial[pion][0];
+
+ if (status == RESTORE) {
+ check_a_neighbor();
+ return;
+ }
+
+ for (neigh = 0; 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 nfirst_neighs = onemol->nspecial[pion][0];
+
+ if (status == RESTORE) {
+ inner_crosscheck_loop();
+ return;
+ }
+
+ for (trace = 0; 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()
+{
+ 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("limit_tags",flag);
+ int *i_limit_tags = atom->ivector[index1];
+
+ int index2 = atom->find_custom("statted_tags",flag);
+ int *i_statted_tags = atom->ivector[index2];
+
+ int index3 = atom->find_custom("react_tags",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 0 (which removes from dynamic group)
+ int flag;
+ int index1 = atom->find_custom("limit_tags",flag);
+ int *i_limit_tags = atom->ivector[index1];
+
+ int index2 = atom->find_custom("statted_tags",flag);
+ int *i_statted_tags = atom->ivector[index2];
+
+ int index3 = atom->find_custom("react_tags",flag);
+ int *i_react_tags = atom->ivector[index3];
+
+ 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) {
+ delete [] allncols;
+ 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(tagint), &column);
+ MPI_Type_commit(&column);
+
+ 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
+ if (me == 0) {
+ MPI_Gatherv(MPI_IN_PLACE, ghostly_num_mega, column, // Note: some values ignored for MPI_IN_PLACE
+ &(global_mega_glove[0][0]), allncols, allstarts,
+ column, 0, world);
+ } else {
+ 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;
+
+ MPI_Type_free(&column);
+ MPI_Type_free(&columnunsized);
+#endif
+}
+
+/* ----------------------------------------------------------------------
+update charges, types, special lists and all topology
+------------------------------------------------------------------------- */
+
+void FixBondReact::update_everything()
+{
+ int *type = atom->type;
+
+ int nlocal = atom->nlocal;
+ 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_INT,MPI_SUM,world);
+ atom->nbonds += Tdelta_bonds;
+
+ int Tdelta_angle;
+ MPI_Allreduce(&delta_angle,&Tdelta_angle,1,MPI_INT,MPI_SUM,world);
+ atom->nangles += Tdelta_angle;
+
+ int Tdelta_dihed;
+ MPI_Allreduce(&delta_dihed,&Tdelta_dihed,1,MPI_INT,MPI_SUM,world);
+ atom->ndihedrals += Tdelta_dihed;
+
+ int Tdelta_imprp;
+ MPI_Allreduce(&delta_imprp,&Tdelta_imprp,1,MPI_INT,MPI_SUM,world);
+ atom->nimpropers += Tdelta_imprp;
+}
+
+/* ----------------------------------------------------------------------
+read superimpose file
+------------------------------------------------------------------------- */
+
+void FixBondReact::read(int myrxn)
+{
+ 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/USER-MISC/fix_bond_react.h b/src/USER-MISC/fix_bond_react.h
new file mode 100644
index 0000000000000000000000000000000000000000..a0f52c02882353de6c9ddf2794b34e6a58f4b91a
--- /dev/null
+++ b/src/USER-MISC/fix_bond_react.h
@@ -0,0 +1,216 @@
+/* -*- 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
+ Fix *fix1; // nve/limit used to relax reaction sites
+ Fix *fix2; // properties/atom used to indicate 1) indicate relaxing atoms
+ // 2) system-wide thermostat
+ // 3) 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 properties
+ 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.
+
+*/