Merge pull request #890 from jrgissing/bond_react_final_touches

[Syntax:]

fix ID group-ID bond/react common_keyword values ...

  react react-ID react-group-ID Nevery Rmin template-ID(pre-reacted) template-ID(post-reacted) map_file individual_keyword values ...

  react react-ID react-group-ID Nevery Rmin template-ID(pre-reacted) template-ID(post-reacted) map_file individual_keyword values ...

  react react-ID react-group-ID Nevery Rmin template-ID(pre-reacted) template-ID(post-reacted) map_file individual_keyword values ...

  react react-ID react-group-ID Nevery Rmin Rmax template-ID(pre-reacted) template-ID(post-reacted) map_file individual_keyword values ...

  react react-ID react-group-ID Nevery Rmin Rmax template-ID(pre-reacted) template-ID(post-reacted) map_file individual_keyword values ...

  react react-ID react-group-ID Nevery Rmin Rmax template-ID(pre-reacted) template-ID(post-reacted) map_file individual_keyword values ...

  ... :pre

ID, group-ID are documented in "fix"_fix.html command. Group-ID is ignored. :ulb,l

  {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

    distance = xmax value that is used by an internally created "nve/limit"_fix_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

  Rmin = bonding pair atoms must be separated by more than Rmin to initiate reaction (distance units) :l

  Rmax = bonding pair atoms must be separated by less than Rmax to initiate reaction (distance units) :l

  template-ID(pre-reacted) = ID of a molecule template containing pre-reaction topology :l

  template-ID(post-reacted) = ID of a molecule template containing post-reaction topology :l

  map_file = name of file specifying corresponding atomIDs in the pre- and post-reacted templates :l

molecule mol1 pre_reacted_topology.txt

molecule mol2 post_reacted_topology.txt

fix 5 all bond/react stabilization no react myrxn1 all 1 3.25 mol1 mol2 map_file.txt

fix 5 all bond/react stabilization no react myrxn1 all 1 0 3.25 mol1 mol2 map_file.txt

molecule mol1 pre_reacted_rxn1.txt

molecule mol2 post_reacted_rxn1.txt

molecule mol3 pre_reacted_rxn2.txt

molecule mol4 post_reacted_rxn2.txt

fix 5 all bond/react stabilization yes nvt_grp .03 &

  react myrxn1 all 1 3.25 mol1 mol2 map_file_rxn1.txt prob 0.50 12345 &

  react myrxn2 all 1 2.75 mol3 mol4 map_file_rxn2.txt prob 0.25 12345

  react myrxn1 all 1 0 3.25 mol1 mol2 map_file_rxn1.txt prob 0.50 12345 &

  react myrxn2 all 1 0 2.75 mol3 mol4 map_file_rxn2.txt prob 0.25 12345

fix 6 nvt_grp nvt temp 300 300 100 # system-wide thermostat must be defined after bond/react :pre

[Description:]

previously unused group-ID. The fix bond/react command creates a

"dynamic group"_group.html of this name that excludes reacting atoms.

This dynamic group-ID should then be used by a subsequent system-wide

time integrator, as shown in the second example above. It is necessary

to place the time integration command after the fix bond/react command

due to the internal dynamic grouping performed by fix bond/react.

time integrator, as shown in the second example above. It is currently

necessary to place the time integration command after the fix

bond/react command due to the internal dynamic grouping performed by

fix bond/react.

NOTE: The internally created group currently applies to all atoms in

the system, i.e. you should generally not have a separate thermostat

which acts on the 'all' group.

The following comments pertain to each 'react' argument:

A bonding atom pair will be identified if several conditions are met.

First, a pair of atoms within the specified react-group-ID of type

typei and typej must be within a distance Rmin of each other. The atom

types typei and typej are specified in the pre- and post-reaction

templates. The distance calculation uses the pair neighbor list,

therefore bonded neighbor exclusions may prevent a reaction between

1st, 2nd or 3rd bonded neighbor atoms. If multiple bonding atom pairs

are identified for an atom, the closest bonding atom partner is set as

its "nearest" bonding partner. Then, if both an atomi and atomj have

each other as their nearest bonding partners, these two atoms are

identified as the bonding atom pair of the reaction site. Once this

unique bonding atom pair is identified for each reaction, there could

two or more reactions that involve a given atom on the same timestep.

If this is the case, only one such reaction is permitted to occur.

This reaction is chosen randomly from all potential reactions. This

capability allows e.g. for different reaction pathways to proceed from

identical reaction sites with user-specified probabilities.

typei and typej must separated by a distance between Rmin and Rmax. It

is possible that multiple bonding atom pairs are identified: if the

bonding atoms in the pre-reacted template are not 1-2, 1-3, or 1-4

neighbors, the closest bonding atom partner is set as its bonding

partner; otherwise, the farthest potential partner is chosen. Then, if

both an atomi and atomj have each other as their nearest bonding

partners, these two atoms are identified as the bonding atom pair of

the reaction site. Once this unique bonding atom pair is identified

for each reaction, there could two or more reactions that involve a

given atom on the same timestep. If this is the case, only one such

reaction is permitted to occur. This reaction is chosen randomly from

all potential reactions. This capability allows e.g. for different

reaction pathways to proceed from identical reaction sites with

user-specified probabilities.

The pre-reacted molecule template is specified by a molecule command.

This molecule template file contains a sample reaction site and its

A few capabilities to note: 1) You may specify as many 'react'

arguments as desired. For example, you could break down a complicated

reaction mechanism into several reaction steps, each defined by its

own 'react' argument. 2) While typically a bond is formed between the

bonding atom pairs specified in the pre-reacted molecule template,

this is not required.

own 'react' argument. 2) While typically a bond is formed or removed

between the bonding atom pairs specified in the pre-reacted molecule

template, this is not required. 3) By reversing the order of the pre-

and post- reacted molecule templates in another 'react' argument, you

can allow for the possibility of one or more reverse reactions.

The optional keywords deal with the probability of a given reaction

occurring as well as the stable equilibration of each reaction site as

fix 1 bond_react_MASTER_group temp/rescale 1 300 300 10 1

NOTE: This command must be added after the fix bond/react command, and

will apply to all reaction steps.

will apply to all reactions.

Computationally, each timestep this fix operates, it loops over

neighbor lists and computes distances between pairs of atoms in the

list. It also communicates between neighboring processors to

coordinate which bonds are created. All of these operations increase

the cost of a timestep. Thus you should be cautious about invoking

this fix too frequently.

neighbor lists (for bond-forming reactions) and computes distances

between pairs of atoms in the list. It also communicates between

neighboring processors to coordinate which bonds  are created and/or

removed. All of these operations increase the cost of a timestep. Thus

you should be cautious about invoking this fix too frequently.

You can dump out snapshots of the current bond topology via the dump

local command.

# dump 1 all xyz 100 test_vis.xyz

fix myrxns all bond/react stabilization yes statted_grp .03 &

  react rxn1 all 1 2.9 mol1 mol2 rxn1_stp1_map &

  react rxn2 all 1 5 mol3 mol4 rxn1_stp2_map

  react rxn1 all 1 0.0 2.9 mol1 mol2 rxn1_stp1_map &

  react rxn2 all 1 0.0 5.0 mol3 mol4 rxn1_stp2_map

# stable at 800K

fix 1 statted_grp nvt temp 800 800 100

# 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

fix myrxns all bond/react stabilization yes statted_grp .03   react rxn1 all 1 0.0 2.9 mol1 mol2 rxn1_stp1_map   react rxn2 all 1 0.0 5.0 mol3 mol4 rxn1_stp2_map

WARNING: An atom in 'react #1' changes bond connectivity but not atom type (../fix_bond_react.cpp:1489)

WARNING: An atom in 'react #2' changes bond connectivity but not atom type (../fix_bond_react.cpp:1489)

dynamic group bond_react_MASTER_group defined

dynamic group statted_grp defined

dynamic group bond_react_MASTER_group defined

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

     100    697.22819   -4624.0512   0.80366765           57           22 

     150    723.60507   -17175.571   0.80366765           76           48 

     200    736.71277   -12961.963   0.80366765           84           64 

Loop time of 102.825 on 1 procs for 200 steps with 35200 atoms

Performance: 0.161 ns/day, 149.151 hours/ns, 1.862 timesteps/s

Performance: 0.168 ns/day, 142.812 hours/ns, 1.945 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

Pair    | 27.193     | 27.193     | 27.193     |   0.0 | 26.45

Bond    | 11.324     | 11.324     | 11.324     |   0.0 | 11.01

Kspace  | 4.1878     | 4.1878     | 4.1878     |   0.0 |  4.07

Neigh   | 54.724     | 54.724     | 54.724     |   0.0 | 53.22

Comm    | 0.40662    | 0.40662    | 0.40662    |   0.0 |  0.40

Output  | 0.0011101  | 0.0011101  | 0.0011101  |   0.0 |  0.00

Modify  | 4.9422     | 4.9422     | 4.9422     |   0.0 |  4.81

Other   |            | 0.04545    |            |       |  0.04

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

Nghost:    38403 ave 38403 max 38403 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

Neighs:    6.9281e+06 ave 6.9281e+06 max 6.9281e+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

Total # of neighbors = 6928101

Ave neighs/atom = 196.821

Ave special neighs/atom = 9.83727

Neighbor list builds = 200

Dangerous builds = 0

# write_restart restart_longrun

# write_data restart_longrun.data

Total wall time: 0:01:48

Please see the log.cite file for references relevant to this simulation

Total wall time: 0:01:43

# 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

fix myrxns all bond/react stabilization yes statted_grp .03   react rxn1 all 1 0.0 2.9 mol1 mol2 rxn1_stp1_map   react rxn2 all 1 0.0 5.0 mol3 mol4 rxn1_stp2_map

WARNING: An atom in 'react #1' changes bond connectivity but not atom type (../fix_bond_react.cpp:1489)

WARNING: An atom in 'react #2' changes bond connectivity but not atom type (../fix_bond_react.cpp:1489)

dynamic group bond_react_MASTER_group defined

dynamic group statted_grp defined

dynamic group bond_react_MASTER_group defined

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

     100    697.22819   -4624.0512   0.80366765           57           22 

     150    724.40407   -17166.729   0.80366765           76           49 

     200    737.28582   -12968.224   0.80366765           84           65 

Loop time of 51.171 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

Performance: 0.338 ns/day, 71.071 hours/ns, 3.908 timesteps/s

98.4% 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

Pair    | 12.926     | 13.247     | 13.493     |   6.7 | 25.89

Bond    | 5.2132     | 5.2733     | 5.3367     |   1.9 | 10.31

Kspace  | 2.3601     | 2.6534     | 3.0067     |  16.0 |  5.19

Neigh   | 25.93      | 25.934     | 25.937     |   0.1 | 50.68

Comm    | 0.73273    | 0.75464    | 0.78505    |   2.3 |  1.47

Output  | 0.00045228 | 0.00067407 | 0.0013323  |   0.0 |  0.00

Modify  | 3.2682     | 3.2686     | 3.2692     |   0.0 |  6.39

Other   |            | 0.03995    |            |       |  0.08

Nlocal:    8800 ave 8913 max 8652 min

Histogram: 1 0 0 0 1 0 0 0 1 1

Nghost:    18366 ave 18461 max 18190 min

Histogram: 1 0 0 0 0 0 0 1 1 1

Neighs:    1.73203e+06 ave 1.77261e+06 max 1.68165e+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

Total # of neighbors = 6928132

Ave neighs/atom = 196.822

Ave special neighs/atom = 9.83608

Neighbor list builds = 200

Dangerous builds = 0

# write_restart restart_longrun

# write_data restart_longrun.data

Total wall time: 0:00:57

Please see the log.cite file for references relevant to this simulation

Total wall time: 0:00:52

# dump 1 all xyz 1 test_vis.xyz

fix myrxns all bond/react stabilization yes statted_grp .03 &

  react rxn1 all 1 2.9 mol1 mol2 rxn1_stp1_map &

  react rxn2 all 1 5 mol3 mol4 rxn1_stp2_map

  react rxn1 all 1 0.0 2.9 mol1 mol2 rxn1_stp1_map &

  react rxn2 all 1 0.0 5.0 mol3 mol4 rxn1_stp2_map

fix 1 statted_grp nvt temp 300 300 100