diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst index 97717f59fc..b995239d08 100644 --- a/doc/src/fix_bond_react.rst +++ b/doc/src/fix_bond_react.rst @@ -41,7 +41,7 @@ Syntax * template-ID(post-reacted) = ID of a molecule template containing post-reaction topology * map_file = name of file specifying corresponding atom-IDs in the pre- and post-reacted templates * zero or more individual keyword/value pairs may be appended to each react argument -* individual_keyword = *prob* or *max_rxn* or *stabilize_steps* or *custom_charges* +* individual_keyword = *prob* or *max_rxn* or *stabilize_steps* or *custom_charges* or *molecule* or *modify_create* .. parsed-literal:: @@ -59,6 +59,12 @@ Syntax off = allow both inter- and intramolecular reactions (default) inter = search for reactions between molecules with different IDs intra = search for reactions within the same molecule + *modify_create* keyword values + *fit* value = *all* or *fragmentID* + all = use all eligible atoms for create-atoms fit (default) + fragmentID = ID of molecule fragment used for create-atoms fit + *overlap* value = R + R = only insert atom/molecule if further than R from existing particles (distance units) Examples """""""" @@ -89,7 +95,9 @@ 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. In addition, reaction by-products or -other molecules can be identified and deleted. +other molecules can be identified and deleted. Finally, atoms can be +created and inserted at specific positions relative to the reaction +site. Fix bond/react does not use quantum mechanical (eg. fix qmmm) or pairwise bond-order potential (eg. Tersoff or AIREBO) methods to @@ -262,14 +270,14 @@ 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 -:doc:`molecule ` command page. +must contain the same number of atoms as the pre-reacted template +(unless there are created atoms). 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 :doc:`molecule ` command page. .. note:: @@ -283,7 +291,7 @@ A discussion of correctly handling this is also provided on the The map file is a text document with the following format: A map file has a header and a body. The header of map file the -contains one mandatory keyword and four optional keywords. The +contains one mandatory keyword and five optional keywords. The mandatory keyword is 'equivalences': .. parsed-literal:: @@ -296,11 +304,12 @@ The optional keywords are 'edgeIDs', 'deleteIDs', 'chiralIDs' and .. parsed-literal:: N *edgeIDs* = # of edge atoms N in the pre-reacted molecule template - N *deleteIDs* = # of atoms N that are specified for deletion - N *chiralIDs* = # of specified chiral centers N - N *constraints* = # of specified reaction constraints N + N *deleteIDs* = # of atoms N that are deleted + N *createIDs* = # of atoms N that are created + N *chiralIDs* = # of chiral centers N + N *constraints* = # of reaction constraints N -The body of the map file contains two mandatory sections and four +The body of the map file contains two mandatory sections and five optional sections. The first mandatory section begins with the keyword 'InitiatorIDs' and lists the two atom IDs of the initiator atom pair in the pre-reacted molecule template. The second mandatory section @@ -313,8 +322,10 @@ the keyword 'EdgeIDs' and lists the atom IDs of edge atoms in the pre-reacted molecule template. The second optional section begins with the keyword 'DeleteIDs' and lists the atom IDs of pre-reaction template atoms to delete. The third optional section begins with the +keyword 'CreateIDs' and lists the atom IDs of the post-reaction +template atoms to create. The fourth optional section begins with the keyword 'ChiralIDs' lists the atom IDs of chiral atoms whose -handedness should be enforced. The fourth optional section begins with +handedness should be enforced. The fifth optional section begins with the keyword 'Constraints' and lists additional criteria that must be satisfied in order for the reaction to occur. Currently, there are five types of constraints available, as discussed below: 'distance', @@ -353,6 +364,38 @@ A sample map file is given below: ---------- +A user-specified set of atoms can be deleted by listing their +pre-reaction template IDs in the DeleteIDs section. A deleted atom +must still be included in the post-reaction molecule template, in +which it cannot be bonded to an atom that is not deleted. In addition +to deleting unwanted reaction by-products, this feature can be used to +remove specific topologies, such as small rings, that may be otherwise +indistinguishable. + +Atoms can be created by listing their post-reaction template IDs in +the CreateIDs section. A created atom should not be included in the +pre-reaction template. The inserted positions of created atoms are +determined by the coordinates of the post-reaction template, after +optimal translation and rotation of the post-reaction template to the +reaction site (using a fit with atoms that are neither created nor +deleted). The *modify_create* keyword can be used to modify the +default behavior when creating atoms. The *modify_create* keyword has +two sub-keywords, *fit* and *overlap*. One or more of the sub-keywords +may be used after the *modify_create* keyword. The *fit* sub-keyword +can be used to specify which post-reaction atoms are used for the +optimal translation and rotation of the post-reaction template. The +*fragmentID* value of the *fit* sub-keyword must be the name of a +molecule fragment defined in the post-reaction :doc:`molecule +` template, and only atoms in this fragment are used for the +fit. Atoms are created only if no current atom in the simulation is +within a distance R of any created atom, including the effect of +periodic boundary conditions if applicable. R is defined by the +*overlap* sub-keyword. Note that the default value for R is 0.0, which +will allow atoms to strongly overlap if you are inserting where other +atoms are present. The velocity of each created atom is initialized in +a random direction with a magnitude calculated from the instantaneous +temperature of the reaction site. + The handedness of atoms that are chiral centers can be enforced by listing their IDs in the ChiralIDs section. A chiral atom must be bonded to four atoms with mutually different atom types. This feature @@ -528,15 +571,6 @@ the same molecule ID are considered for the reaction. A few other considerations: -Many reactions result in one or more atoms that are considered -unwanted by-products. Therefore, bond/react provides the option to -delete a user-specified set of atoms. These pre-reaction atoms are -identified in the map file. A deleted atom must still be included in -the post-reaction molecule template, in which it cannot be bonded to -an atom that is not deleted. In addition to deleting unwanted reaction -by-products, this feature can be used to remove specific topologies, -such as small rings, that may be otherwise indistinguishable. - Optionally, you can enforce additional behaviors on reacting atoms. For example, it may be beneficial to force reacting atoms to remain at a certain temperature. For this, you can use the internally-created @@ -610,7 +644,7 @@ Default """"""" The option defaults are stabilization = no, prob = 1.0, stabilize_steps = 60, -reset_mol_ids = yes, custom_charges = no, molecule = off +reset_mol_ids = yes, custom_charges = no, molecule = off, modify_create = no ---------- diff --git a/examples/USER/reaction/create_atoms_polystyrene/grow_styrene.map b/examples/USER/reaction/create_atoms_polystyrene/grow_styrene.map new file mode 100644 index 0000000000..88d282690c --- /dev/null +++ b/examples/USER/reaction/create_atoms_polystyrene/grow_styrene.map @@ -0,0 +1,66 @@ +map file: styrene growth + +1 edgeIDs +30 equivalences +16 createIDs + +InitiatorIDs + +4 +13 + +EdgeIDs + +30 + +CreateIDs + +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 + +Equivalences + +1 45 +2 46 +3 44 +4 43 +5 42 +6 41 +7 40 +8 39 +9 38 +10 37 +11 36 +12 35 +13 34 +14 33 +15 32 +16 31 +17 17 +18 18 +19 19 +20 20 +21 21 +22 22 +23 23 +24 24 +25 25 +26 26 +27 27 +28 28 +29 29 +30 30 diff --git a/examples/USER/reaction/create_atoms_polystyrene/grow_styrene_post.data_template b/examples/USER/reaction/create_atoms_polystyrene/grow_styrene_post.data_template new file mode 100755 index 0000000000..de0c2383bb --- /dev/null +++ b/examples/USER/reaction/create_atoms_polystyrene/grow_styrene_post.data_template @@ -0,0 +1,456 @@ +molecule template: end of chain plus polymerized styrene + +46 atoms +48 bonds +81 angles +121 dihedrals +35 impropers +1 fragments + +Fragments + +create_fit 34 44 + +Types + +1 1 +2 2 +3 1 +4 5 +5 1 +6 2 +7 1 +8 2 +9 1 +10 2 +11 1 +12 2 +13 2 +14 6 +15 2 +16 2 +17 1 +18 2 +19 1 +20 5 +21 1 +22 2 +23 1 +24 2 +25 1 +26 2 +27 1 +28 2 +29 2 +30 6 +31 1 +32 2 +33 1 +34 5 +35 1 +36 2 +37 1 +38 2 +39 1 +40 2 +41 1 +42 2 +43 2 +44 6 +45 2 +46 2 + +Charges + +1 -0.129000 +2 0.123700 +3 0.026600 +4 -0.018200 +5 -0.129000 +6 0.123700 +7 -0.173400 +8 0.140300 +9 -0.113400 +10 0.128800 +11 -0.173400 +12 0.140300 +13 0.051600 +14 -0.069600 +15 0.035400 +16 0.035400 +17 -0.129000 +18 0.123700 +19 0.026600 +20 -0.018200 +21 -0.129000 +22 0.123700 +23 -0.173400 +24 0.140300 +25 -0.113400 +26 0.128800 +27 -0.173400 +28 0.140300 +29 0.051600 +30 -0.069600 +31 -0.129000 +32 0.123700 +33 0.026600 +34 -0.018200 +35 -0.129000 +36 0.123700 +37 -0.173400 +38 0.140300 +39 -0.113400 +40 0.128800 +41 -0.173400 +42 0.140300 +43 0.051600 +44 -0.069600 +45 0.035400 +46 0.035400 + +Coords + +1 24.130699 1.043900 -1.309300 +2 25.062700 1.582900 -1.309300 +3 22.900700 1.753900 -1.309300 +4 22.900700 3.253900 -1.309300 +5 21.670700 1.043900 -1.309300 +6 20.738701 1.582900 -1.309300 +7 21.670700 -0.376100 -1.309300 +8 20.738701 -0.915100 -1.309300 +9 22.900700 -1.086100 -1.309300 +10 22.900700 -2.163100 -1.309300 +11 24.130699 -0.376100 -1.309300 +12 25.062700 -0.915100 -1.309300 +13 23.766701 3.658900 -0.952300 +14 21.622700 3.802900 -1.871300 +15 21.672701 4.544900 -1.970300 +16 20.979700 2.979900 -2.165300 +17 13.465800 0.682500 -1.658900 +18 14.397800 1.221500 -1.658900 +19 12.235800 1.392500 -1.658900 +20 12.235800 2.892500 -1.658900 +21 11.005800 0.682500 -1.658900 +22 10.073800 1.221500 -1.658900 +23 11.005800 -0.737500 -1.658900 +24 10.073800 -1.276500 -1.658900 +25 12.235800 -1.447500 -1.658900 +26 12.235800 -2.524500 -1.658900 +27 13.465800 -0.737500 -1.658900 +28 14.397800 -1.276500 -1.658900 +29 13.101800 3.297500 -1.301900 +30 10.957800 3.441500 -2.220900 +31 18.663500 0.855500 -1.372100 +32 19.595501 1.394500 -1.372100 +33 17.433500 1.565500 -1.372100 +34 17.433500 3.065500 -1.372100 +35 16.203501 0.855500 -1.372100 +36 15.271500 1.394500 -1.372100 +37 16.203501 -0.564500 -1.372100 +38 15.271500 -1.103500 -1.372100 +39 17.433500 -1.274500 -1.372100 +40 17.433500 -2.351500 -1.372100 +41 18.663500 -0.564500 -1.372100 +42 19.595501 -1.103500 -1.372100 +43 18.299500 3.470500 -1.015100 +44 16.155500 3.614500 -1.934100 +45 16.205500 4.356500 -2.033100 +46 15.512500 2.791500 -2.228100 + +Bonds + +1 1 1 2 +2 2 1 3 +3 2 1 11 +4 11 3 4 +5 2 3 5 +6 12 13 4 +7 13 4 14 +8 1 5 6 +9 2 5 7 +10 1 7 8 +11 2 7 9 +12 1 9 10 +13 2 9 11 +14 1 11 12 +15 10 15 14 +16 10 16 14 +17 9 14 34 +18 1 17 18 +19 2 17 19 +20 2 17 27 +21 7 19 20 +22 2 19 21 +23 8 29 20 +24 9 30 20 +25 9 44 20 +26 1 21 22 +27 2 21 23 +28 1 23 24 +29 2 23 25 +30 1 25 26 +31 2 25 27 +32 1 27 28 +33 1 31 32 +34 2 31 33 +35 2 31 41 +36 7 33 34 +37 2 33 35 +38 8 43 34 +39 9 44 34 +40 1 35 36 +41 2 35 37 +42 1 37 38 +43 2 37 39 +44 1 39 40 +45 2 39 41 +46 1 41 42 +47 10 45 44 +48 10 46 44 + +Angles + +1 1 3 1 2 +2 1 11 1 2 +3 2 3 1 11 +4 17 1 3 4 +5 2 1 3 5 +6 17 5 3 4 +7 18 3 4 13 +8 19 3 4 14 +9 20 13 4 14 +10 1 3 5 6 +11 2 3 5 7 +12 1 7 5 6 +13 1 5 7 8 +14 2 5 7 9 +15 1 9 7 8 +16 1 7 9 10 +17 2 7 9 11 +18 1 11 9 10 +19 2 1 11 9 +20 1 1 11 12 +21 1 9 11 12 +22 21 15 14 4 +23 21 16 14 4 +24 22 4 14 34 +25 15 15 14 16 +26 14 15 14 34 +27 14 16 14 34 +28 1 19 17 18 +29 1 27 17 18 +30 2 19 17 27 +31 9 17 19 20 +32 2 17 19 21 +33 9 21 19 20 +34 10 19 20 29 +35 11 19 20 30 +36 11 19 20 44 +37 12 29 20 30 +38 12 29 20 44 +39 13 30 20 44 +40 1 19 21 22 +41 2 19 21 23 +42 1 23 21 22 +43 1 21 23 24 +44 2 21 23 25 +45 1 25 23 24 +46 1 23 25 26 +47 2 23 25 27 +48 1 27 25 26 +49 2 17 27 25 +50 1 17 27 28 +51 1 25 27 28 +52 1 33 31 32 +53 1 41 31 32 +54 2 33 31 41 +55 9 31 33 34 +56 2 31 33 35 +57 9 35 33 34 +58 11 33 34 14 +59 12 43 34 14 +60 13 14 34 44 +61 10 33 34 43 +62 11 33 34 44 +63 12 43 34 44 +64 1 33 35 36 +65 2 33 35 37 +66 1 37 35 36 +67 1 35 37 38 +68 2 35 37 39 +69 1 39 37 38 +70 1 37 39 40 +71 2 37 39 41 +72 1 41 39 40 +73 2 31 41 39 +74 1 31 41 42 +75 1 39 41 42 +76 16 20 44 34 +77 14 45 44 20 +78 14 46 44 20 +79 14 45 44 34 +80 14 46 44 34 +81 15 45 44 46 + +Dihedrals + +1 20 2 1 3 4 +2 2 5 3 1 2 +3 21 11 1 3 4 +4 4 11 1 3 5 +5 2 9 11 1 2 +6 5 2 1 11 12 +7 4 3 1 11 9 +8 2 3 1 11 12 +9 22 1 3 4 13 +10 23 1 3 4 14 +11 22 5 3 4 13 +12 23 5 3 4 14 +13 2 1 3 5 6 +14 4 1 3 5 7 +15 20 6 5 3 4 +16 21 7 5 3 4 +17 24 3 4 14 15 +18 24 3 4 14 16 +19 25 3 4 14 34 +20 26 13 4 14 15 +21 26 13 4 14 16 +22 27 13 4 14 34 +23 2 3 5 7 8 +24 4 3 5 7 9 +25 5 6 5 7 8 +26 2 9 7 5 6 +27 2 5 7 9 10 +28 4 5 7 9 11 +29 5 8 7 9 10 +30 2 11 9 7 8 +31 4 7 9 11 1 +32 2 7 9 11 12 +33 2 1 11 9 10 +34 5 10 9 11 12 +35 28 4 14 34 33 +36 29 4 14 34 43 +37 30 4 14 34 44 +38 31 15 14 34 33 +39 32 15 14 34 43 +40 33 15 14 34 44 +41 31 16 14 34 33 +42 32 16 14 34 43 +43 33 16 14 34 44 +44 10 18 17 19 20 +45 2 21 19 17 18 +46 11 27 17 19 20 +47 4 27 17 19 21 +48 2 25 27 17 18 +49 5 18 17 27 28 +50 4 19 17 27 25 +51 2 19 17 27 28 +52 12 17 19 20 29 +53 13 17 19 20 30 +54 13 17 19 20 44 +55 12 21 19 20 29 +56 13 21 19 20 30 +57 13 21 19 20 44 +58 2 17 19 21 22 +59 4 17 19 21 23 +60 10 22 21 19 20 +61 11 23 21 19 20 +62 34 34 44 20 19 +63 31 45 44 20 19 +64 31 46 44 20 19 +65 35 34 44 20 29 +66 32 45 44 20 29 +67 32 46 44 20 29 +68 36 34 44 20 30 +69 33 45 44 20 30 +70 33 46 44 20 30 +71 2 19 21 23 24 +72 4 19 21 23 25 +73 5 22 21 23 24 +74 2 25 23 21 22 +75 2 21 23 25 26 +76 4 21 23 25 27 +77 5 24 23 25 26 +78 2 27 25 23 24 +79 4 23 25 27 17 +80 2 23 25 27 28 +81 2 17 27 25 26 +82 5 26 25 27 28 +83 10 32 31 33 34 +84 2 35 33 31 32 +85 11 41 31 33 34 +86 4 41 31 33 35 +87 2 39 41 31 32 +88 5 32 31 41 42 +89 4 33 31 41 39 +90 2 33 31 41 42 +91 13 31 33 34 14 +92 12 31 33 34 43 +93 13 31 33 34 44 +94 13 35 33 34 14 +95 12 35 33 34 43 +96 13 35 33 34 44 +97 2 31 33 35 36 +98 4 31 33 35 37 +99 10 36 35 33 34 +100 11 37 35 33 34 +101 36 20 44 34 14 +102 33 45 44 34 14 +103 33 46 44 34 14 +104 34 20 44 34 33 +105 31 45 44 34 33 +106 31 46 44 34 33 +107 35 20 44 34 43 +108 32 45 44 34 43 +109 32 46 44 34 43 +110 2 33 35 37 38 +111 4 33 35 37 39 +112 5 36 35 37 38 +113 2 39 37 35 36 +114 2 35 37 39 40 +115 4 35 37 39 41 +116 5 38 37 39 40 +117 2 41 39 37 38 +118 4 37 39 41 31 +119 2 37 39 41 42 +120 2 31 41 39 40 +121 5 40 39 41 42 + +Impropers + +1 1 3 1 11 2 +2 8 1 3 5 4 +3 9 3 4 13 14 +4 1 3 5 7 6 +5 1 5 7 9 8 +6 1 7 9 11 10 +7 1 1 11 9 12 +8 1 19 17 27 18 +9 5 17 19 21 20 +10 1 19 21 23 22 +11 1 21 23 25 24 +12 1 23 25 27 26 +13 1 17 27 25 28 +14 1 33 31 41 32 +15 5 31 33 35 34 +16 1 33 35 37 36 +17 1 35 37 39 38 +18 1 37 39 41 40 +19 1 31 41 39 42 +20 1 15 14 16 4 +21 1 15 14 4 34 +22 1 16 14 4 34 +23 1 15 14 16 34 +24 1 19 20 29 30 +25 1 19 20 29 44 +26 1 19 20 30 44 +27 1 29 20 30 44 +28 1 33 34 43 14 +29 1 33 34 14 44 +30 1 43 34 14 44 +31 1 33 34 43 44 +32 1 45 44 34 20 +33 1 46 44 34 20 +34 1 45 44 46 20 +35 1 45 44 46 34 diff --git a/examples/USER/reaction/create_atoms_polystyrene/grow_styrene_pre.data_template b/examples/USER/reaction/create_atoms_polystyrene/grow_styrene_pre.data_template new file mode 100644 index 0000000000..d04fefccf5 --- /dev/null +++ b/examples/USER/reaction/create_atoms_polystyrene/grow_styrene_pre.data_template @@ -0,0 +1,294 @@ +molecule template: end of styrene chain + +30 atoms +31 bonds +51 angles +73 dihedrals +21 impropers + +Types + +1 2 +2 2 +3 6 +4 2 +5 2 +6 1 +7 2 +8 1 +9 2 +10 1 +11 2 +12 1 +13 5 +14 1 +15 2 +16 1 +17 1 +18 2 +19 1 +20 5 +21 1 +22 2 +23 1 +24 2 +25 1 +26 2 +27 1 +28 2 +29 2 +30 6 + +Coords + +1 59.89981112372972 62.733697275315585 59.09884284578856 +2 61.41970248324232 63.42116581894993 59.52874545893742 +3 60.864754970096406 62.91724243011892 59.559720865992695 +4 62.139819000186826 61.41011937002877 60.81065044071466 +5 60.036455711425084 57.160029629288026 60.31958663310848 +6 59.734195751174056 58.18706337912225 60.20562410798949 +7 57.64574781117771 57.712432799329 59.860109977091554 +8 58.37408644866664 58.50134169314242 59.94422053768215 +9 56.94300092269842 60.093170109004795 59.5955638127831 +10 57.974275786582744 59.85577775892068 59.793714995577716 +11 58.63231375134033 61.922969938852454 59.79065033121885 +12 58.934573711591355 60.89593618901822 59.904612856337835 +13 61.30908151524225 61.68041745837013 60.28316188676589 +14 60.29468229868386 60.58165855333751 60.16601625920239 +15 61.725768540066994 58.98982945913568 60.51467315154424 +16 60.69449367618267 59.2272218092198 60.31652196874961 +17 56.90935800040509 62.609851248143706 59.150831390216375 +18 57.940632148874506 62.37245957639904 59.3489824055682 +19 56.509546622906285 63.96428799226142 59.00032568066915 +20 57.52394583946467 65.06304689729403 59.11747130823266 +21 55.14943732039887 64.27856630628159 58.738922110361806 +22 54.84717807556275 65.30559937777636 58.62495975268562 +23 54.18913939539026 63.23840787618404 58.62802424960169 +24 53.15786524692084 63.4757995479287 58.42987323424986 +25 54.58895077288906 61.88397113206633 58.77852995914891 +26 53.86061213540014 61.09506223825291 58.69441939855832 +27 55.94906007539648 61.56969281804616 59.039933529456256 +28 56.2513193202326 60.54265974655139 59.15389588713244 +29 58.35468332440925 64.79274880895268 59.64495986218142 +30 57.07961929431883 66.29987186904283 58.394030287459465 + +Charges + +1 0.0354 +2 0.0354 +3 -0.0696 +4 0.0516 +5 0.1403 +6 -0.1734 +7 0.1288 +8 -0.1134 +9 0.1403 +10 -0.1734 +11 0.1237 +12 -0.129 +13 -0.0182 +14 0.0266 +15 0.1237 +16 -0.129 +17 -0.129 +18 0.1237 +19 0.0266 +20 -0.0182 +21 -0.129 +22 0.1237 +23 -0.1734 +24 0.1403 +25 -0.1134 +26 0.1288 +27 -0.1734 +28 0.1403 +29 0.0516 +30 -0.0696 + +Bonds + +1 10 1 3 +2 10 2 3 +3 8 4 13 +4 1 6 5 +5 1 8 7 +6 2 8 6 +7 1 10 9 +8 2 10 8 +9 1 12 11 +10 2 12 10 +11 9 13 3 +12 7 14 13 +13 2 14 12 +14 1 16 15 +15 2 16 14 +16 2 16 6 +17 1 17 18 +18 2 17 19 +19 2 17 27 +20 7 19 20 +21 2 19 21 +22 9 20 30 +23 9 20 3 +24 1 21 22 +25 2 21 23 +26 1 23 24 +27 2 23 25 +28 1 25 26 +29 2 25 27 +30 1 27 28 +31 8 29 20 + +Angles + +1 16 20 3 13 +2 14 2 3 20 +3 14 1 3 20 +4 14 2 3 13 +5 14 1 3 13 +6 15 2 3 1 +7 2 16 6 8 +8 1 16 6 5 +9 1 8 6 5 +10 1 10 8 7 +11 2 10 8 6 +12 1 6 8 7 +13 1 12 10 9 +14 2 12 10 8 +15 1 8 10 9 +16 1 14 12 11 +17 2 14 12 10 +18 1 10 12 11 +19 10 14 13 4 +20 11 14 13 3 +21 12 4 13 3 +22 9 16 14 13 +23 2 16 14 12 +24 9 12 14 13 +25 1 14 16 15 +26 1 6 16 15 +27 2 14 16 6 +28 1 19 17 18 +29 1 27 17 18 +30 2 19 17 27 +31 9 17 19 20 +32 2 17 19 21 +33 9 21 19 20 +34 10 19 20 29 +35 11 19 20 30 +36 11 19 20 3 +37 12 29 20 30 +38 12 29 20 3 +39 13 30 20 3 +40 1 19 21 22 +41 2 19 21 23 +42 1 23 21 22 +43 1 21 23 24 +44 2 21 23 25 +45 1 25 23 24 +46 1 23 25 26 +47 2 23 25 27 +48 1 27 25 26 +49 2 17 27 25 +50 1 17 27 28 +51 1 25 27 28 + +Dihedrals + +1 2 8 6 16 15 +2 2 16 6 8 7 +3 2 6 8 10 9 +4 4 10 8 6 16 +5 2 10 8 6 5 +6 5 7 8 6 5 +7 2 8 10 12 11 +8 2 12 10 8 7 +9 4 12 10 8 6 +10 5 9 10 8 7 +11 10 11 12 14 13 +12 11 10 12 14 13 +13 2 14 12 10 9 +14 4 14 12 10 8 +15 5 11 12 10 9 +16 17 14 13 3 20 +17 14 14 13 3 2 +18 14 14 13 3 1 +19 18 4 13 3 20 +20 15 4 13 3 2 +21 15 4 13 3 1 +22 2 12 14 16 15 +23 12 16 14 13 4 +24 13 16 14 13 3 +25 12 12 14 13 4 +26 13 12 14 13 3 +27 2 16 14 12 11 +28 4 16 14 12 10 +29 10 15 16 14 13 +30 11 6 16 14 13 +31 4 6 16 14 12 +32 5 15 16 6 5 +33 4 14 16 6 8 +34 2 14 16 6 5 +35 10 18 17 19 20 +36 11 27 17 19 20 +37 4 27 17 19 21 +38 5 18 17 27 28 +39 4 19 17 27 25 +40 2 19 17 27 28 +41 2 21 19 17 18 +42 12 17 19 20 29 +43 13 17 19 20 30 +44 13 17 19 20 3 +45 12 21 19 20 29 +46 13 21 19 20 30 +47 13 21 19 20 3 +48 2 17 19 21 22 +49 4 17 19 21 23 +50 17 19 20 3 13 +51 14 19 20 3 2 +52 14 19 20 3 1 +53 18 29 20 3 13 +54 15 29 20 3 2 +55 15 29 20 3 1 +56 19 30 20 3 13 +57 16 30 20 3 2 +58 16 30 20 3 1 +59 10 22 21 19 20 +60 11 23 21 19 20 +61 2 19 21 23 24 +62 4 19 21 23 25 +63 5 22 21 23 24 +64 2 25 23 21 22 +65 2 21 23 25 26 +66 4 21 23 25 27 +67 5 24 23 25 26 +68 2 27 25 23 24 +69 4 23 25 27 17 +70 2 23 25 27 28 +71 5 26 25 27 28 +72 2 25 27 17 18 +73 2 17 27 25 26 + +Impropers + +1 1 2 3 13 20 +2 1 1 3 13 20 +3 1 2 3 1 20 +4 1 2 3 1 13 +5 1 16 6 8 5 +6 1 10 8 6 7 +7 1 12 10 8 9 +8 1 14 12 10 11 +9 7 14 13 4 3 +10 5 16 14 12 13 +11 1 14 16 6 15 +12 1 19 17 27 18 +13 5 17 19 21 20 +14 1 19 20 29 30 +15 1 19 20 29 3 +16 1 19 20 30 3 +17 1 29 20 30 3 +18 1 19 21 23 22 +19 1 21 23 25 24 +20 1 23 25 27 26 +21 1 17 27 25 28 diff --git a/examples/USER/reaction/create_atoms_polystyrene/in.grow_styrene b/examples/USER/reaction/create_atoms_polystyrene/in.grow_styrene new file mode 100755 index 0000000000..8950cec614 --- /dev/null +++ b/examples/USER/reaction/create_atoms_polystyrene/in.grow_styrene @@ -0,0 +1,48 @@ +# use bond/react 'create atoms' feature to add 30 new styrene monomers to chain + +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 + +variable T equal 530 + +read_data trimer.data & + extra/bond/per/atom 5 & + extra/angle/per/atom 15 & + extra/dihedral/per/atom 15 & + extra/improper/per/atom 25 & + extra/special/per/atom 25 + +molecule mol1 grow_styrene_pre.data_template +molecule mol2 grow_styrene_post.data_template + +fix myrxns all bond/react stabilization yes statted_grp .03 & + react rxn1 all 1 0 3.0 mol1 mol2 grow_styrene.map & + modify_create fit create_fit overlap 2.0 & + stabilize_steps 100 max_rxn 30 + +fix 1 statted_grp_REACT nvt temp $T $T 100 + +fix 4 bond_react_MASTER_group temp/rescale 1 $T $T 1 1 + +thermo_style custom step temp press density f_myrxns[1] + +thermo 100 + +run 8000 + +# write_data final.data nofix diff --git a/examples/USER/reaction/create_atoms_polystyrene/log.24Dec20.grow_styrene.g++.1 b/examples/USER/reaction/create_atoms_polystyrene/log.24Dec20.grow_styrene.g++.1 new file mode 100644 index 0000000000..5f1f2c6698 --- /dev/null +++ b/examples/USER/reaction/create_atoms_polystyrene/log.24Dec20.grow_styrene.g++.1 @@ -0,0 +1,196 @@ +LAMMPS (24 Dec 2020) +Reading data file ... + orthogonal box = (50.000000 50.000000 50.000000) to (250.00000 250.00000 250.00000) + 1 by 1 by 1 MPI processor grid + reading atoms ... + 48 atoms + reading velocities ... + 48 velocities + scanning bonds ... + 8 = max bonds/atom + scanning angles ... + 21 = max angles/atom + scanning dihedrals ... + 33 = max dihedrals/atom + scanning impropers ... + 29 = max impropers/atom + reading bonds ... + 50 bonds + reading angles ... + 84 angles + reading dihedrals ... + 127 dihedrals + reading impropers ... + 36 impropers +Finding 1-2 1-3 1-4 neighbors ... + special bond factors lj: 0 0 0 + special bond factors coul: 0 0 0 + 4 = max # of 1-2 neighbors + 8 = max # of 1-3 neighbors + 17 = max # of 1-4 neighbors + 46 = max # of special neighbors + special bonds CPU = 0.000 seconds + read_data CPU = 0.077 seconds +Read molecule template mol1: + 1 molecules + 30 atoms with max type 6 + 31 bonds with max type 10 + 51 angles with max type 16 + 73 dihedrals with max type 19 + 21 impropers with max type 7 +Read molecule template mol2: + 1 molecules + 46 atoms with max type 6 + 48 bonds with max type 13 + 81 angles with max type 22 + 121 dihedrals with max type 36 + 35 impropers with max type 9 +dynamic group bond_react_MASTER_group defined +dynamic group statted_grp_REACT defined +PPPM initialization ... +WARNING: System is not charge neutral, net charge = -0.00060000000 (../kspace.cpp:324) + using 12-bit tables for long-range coulomb (../kspace.cpp:339) + G vector (1/distance) = 0.20144813 + grid = 45 45 45 + stencil order = 5 + estimated absolute RMS force accuracy = 0.00053712952 + estimated relative force accuracy = 1.6175496e-06 + using double precision KISS FFT + 3d grid and FFT values/proc = 125000 91125 +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 = 39 39 39 + 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 +Setting up Verlet run ... + Unit style : real + Current step : 0 + Time step : 1 +Per MPI rank memory allocation (min/avg/max) = 48.02 | 48.02 | 48.02 Mbytes +Step Temp Press Density f_myrxns[1] + 0 496.23742 0.9983211 6.4856516e-05 0 + 100 534.05394 -0.76952227 6.4856516e-05 0 + 200 552.2225 -0.55375493 6.4856516e-05 0 + 300 857.52834 -0.4272061 8.6475354e-05 1 + 400 714.10681 1.5004615 8.6475354e-05 1 + 500 678.19171 0.21965471 8.6475354e-05 1 + 600 572.3234 0.87879933 8.6475354e-05 1 + 700 996.17398 -0.24269717 0.00010809419 2 + 800 904.50395 1.3662054 0.00010809419 2 + 900 1097.1568 -2.2909907 0.00012971303 3 + 1000 954.08892 1.7705672 0.00012971303 3 + 1100 1102.0377 -1.7018446 0.00015133187 4 + 1200 1239.785 -0.30442903 0.00015133187 4 + 1300 1388.4127 1.3301175 0.00017295071 5 + 1400 1559.3853 1.6709729 0.00017295071 5 + 1500 1471.8623 0.8268427 0.00017295071 5 + 1600 1543.6793 2.1987908 0.00019456955 6 + 1700 1694.5595 0.48852817 0.00019456955 6 + 1800 1632.7737 -1.4617692 0.00021618839 7 + 1900 1922.6502 1.1664257 0.00021618839 7 + 2000 2223.503 -0.95799878 0.00023780722 8 + 2100 2142.6035 0.88444463 0.00025942606 9 + 2200 2298.8636 3.4239313 0.00025942606 9 + 2300 2252.4355 0.82167302 0.00025942606 9 + 2400 2321.0788 1.7499714 0.00025942606 9 + 2500 2095.6715 0.55288444 0.00025942606 9 + 2600 2136.0316 -3.833114 0.00025942606 9 + 2700 2466.3134 -2.2519511 0.00025942606 9 + 2800 2294.3454 1.0637304 0.00025942606 9 + 2900 2340.3891 1.3997049 0.0002810449 10 + 3000 2272.0013 -0.27591886 0.0002810449 10 + 3100 2333.9696 -0.11772138 0.0002810449 10 + 3200 2409.0946 -1.025473 0.0002810449 10 + 3300 2148.023 1.6752329 0.0002810449 10 + 3400 2267.636 -0.45297583 0.0002810449 10 + 3500 2457.622 0.35627297 0.0002810449 10 + 3600 2288.008 -15.516626 0.00030266374 11 + 3700 2458.2681 1.4571773 0.00030266374 11 + 3800 2566.7623 -29.140553 0.00032428258 12 + 3900 2839.4062 0.64583638 0.00032428258 12 + 4000 2893.9852 -52.954497 0.00034590142 13 + 4100 3021.3611 -65.03731 0.00036752025 14 + 4200 3002.7136 1.5750081 0.00036752025 14 + 4300 3218.6248 -120.74039 0.00038913909 15 + 4400 3345.1482 -0.96545269 0.00038913909 15 + 4500 3603.2429 1.2438833 0.00038913909 15 + 4600 3129.8814 -249.91806 0.00041075793 16 + 4700 3769.052 -289.24351 0.00043237677 17 + 4800 3560.4714 -3.1655406 0.00043237677 17 + 4900 3452.2717 -2.1270765 0.00043237677 17 + 5000 3594.3247 -523.48506 0.00045399561 18 + 5100 3578.4199 1.0009097 0.00045399561 18 + 5200 3822.1566 1.0526914 0.00047561445 19 + 5300 3901.8883 -0.14607602 0.00047561445 19 + 5400 4059.3644 -1.7789927 0.00049723329 20 + 5500 4163.6847 1.0240127 0.00049723329 20 + 5600 4109.1649 0.80199787 0.00049723329 20 + 5700 4391.2091 2.8730036 0.00049723329 20 + 5800 4279.6579 -0.36499822 0.00051885212 21 + 5900 4296.2695 -1.3064528 0.00051885212 21 + 6000 4065.3758 -2.0483224 0.00051885212 21 + 6100 4772.5362 -2.6814694 0.00054047096 22 + 6200 4627.029 2.999215 0.0005620898 23 + 6300 5120.7881 0.65372968 0.00058370864 24 + 6400 4588.9559 3.7570705 0.00058370864 24 + 6500 5008.7814 2.3595833 0.00060532748 25 + 6600 5195.0053 1.4641612 0.00060532748 25 + 6700 5622.293 -0.33396047 0.00062694632 26 + 6800 5515.1957 -4.234874 0.00062694632 26 + 6900 5156.7455 0.40171954 0.00064856516 27 + 7000 5120.1639 -1.6065245 0.00064856516 27 + 7100 5650.0327 0.94436323 0.00067018399 28 + 7200 5985.1115 -3.8940347 0.00069180283 29 + 7300 5983.197 0.5293568 0.00069180283 29 + 7400 6001.1559 -0.13712834 0.00071342167 30 + 7500 5889.2134 0.17230892 0.00071342167 30 + 7600 5797.31 2.0920058 0.00071342167 30 + 7700 5865.2783 -0.18556395 0.00071342167 30 + 7800 6207.0659 -5.6237083 0.00071342167 30 + 7900 5627.5108 -2.3718942 0.00071342167 30 + 8000 5823.9502 -0.85418578 0.00071342167 30 +Loop time of 184.87 on 1 procs for 8000 steps with 528 atoms + +Performance: 3.739 ns/day, 6.419 hours/ns, 43.274 timesteps/s +99.9% CPU use with 1 MPI tasks x no OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 3.3043 | 3.3043 | 3.3043 | 0.0 | 1.79 +Bond | 8.0003 | 8.0003 | 8.0003 | 0.0 | 4.33 +Kspace | 168.33 | 168.33 | 168.33 | 0.0 | 91.05 +Neigh | 4.6322 | 4.6322 | 4.6322 | 0.0 | 2.51 +Comm | 0.077927 | 0.077927 | 0.077927 | 0.0 | 0.04 +Output | 0.0020548 | 0.0020548 | 0.0020548 | 0.0 | 0.00 +Modify | 0.5005 | 0.5005 | 0.5005 | 0.0 | 0.27 +Other | | 0.02483 | | | 0.01 + +Nlocal: 528.000 ave 528 max 528 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 341.000 ave 341 max 341 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 35111.0 ave 35111 max 35111 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 35111 +Ave neighs/atom = 66.498106 +Ave special neighs/atom = 11.409091 +Neighbor list builds = 8000 +Dangerous builds = 0 + +Please see the log.cite file for references relevant to this simulation + +Total wall time: 0:03:05 diff --git a/examples/USER/reaction/create_atoms_polystyrene/log.24Dec20.grow_styrene.g++.4 b/examples/USER/reaction/create_atoms_polystyrene/log.24Dec20.grow_styrene.g++.4 new file mode 100644 index 0000000000..8daa6d8161 --- /dev/null +++ b/examples/USER/reaction/create_atoms_polystyrene/log.24Dec20.grow_styrene.g++.4 @@ -0,0 +1,196 @@ +LAMMPS (24 Dec 2020) +Reading data file ... + orthogonal box = (50.000000 50.000000 50.000000) to (250.00000 250.00000 250.00000) + 1 by 2 by 2 MPI processor grid + reading atoms ... + 48 atoms + reading velocities ... + 48 velocities + scanning bonds ... + 8 = max bonds/atom + scanning angles ... + 21 = max angles/atom + scanning dihedrals ... + 33 = max dihedrals/atom + scanning impropers ... + 29 = max impropers/atom + reading bonds ... + 50 bonds + reading angles ... + 84 angles + reading dihedrals ... + 127 dihedrals + reading impropers ... + 36 impropers +Finding 1-2 1-3 1-4 neighbors ... + special bond factors lj: 0 0 0 + special bond factors coul: 0 0 0 + 4 = max # of 1-2 neighbors + 8 = max # of 1-3 neighbors + 17 = max # of 1-4 neighbors + 46 = max # of special neighbors + special bonds CPU = 0.000 seconds + read_data CPU = 0.007 seconds +Read molecule template mol1: + 1 molecules + 30 atoms with max type 6 + 31 bonds with max type 10 + 51 angles with max type 16 + 73 dihedrals with max type 19 + 21 impropers with max type 7 +Read molecule template mol2: + 1 molecules + 46 atoms with max type 6 + 48 bonds with max type 13 + 81 angles with max type 22 + 121 dihedrals with max type 36 + 35 impropers with max type 9 +dynamic group bond_react_MASTER_group defined +dynamic group statted_grp_REACT defined +PPPM initialization ... +WARNING: System is not charge neutral, net charge = -0.00060000000 (../kspace.cpp:324) + using 12-bit tables for long-range coulomb (../kspace.cpp:339) + G vector (1/distance) = 0.20144813 + grid = 45 45 45 + stencil order = 5 + estimated absolute RMS force accuracy = 0.00053712952 + estimated relative force accuracy = 1.6175496e-06 + using double precision KISS FFT + 3d grid and FFT values/proc = 39200 24300 +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 = 39 39 39 + 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 +Setting up Verlet run ... + Unit style : real + Current step : 0 + Time step : 1 +Per MPI rank memory allocation (min/avg/max) = 38.70 | 38.92 | 39.43 Mbytes +Step Temp Press Density f_myrxns[1] + 0 496.23742 0.9983211 6.4856516e-05 0 + 100 534.05394 -0.76952227 6.4856516e-05 0 + 200 552.2225 -0.55375493 6.4856516e-05 0 + 300 857.52834 -0.4272061 8.6475354e-05 1 + 400 714.10681 1.5004615 8.6475354e-05 1 + 500 678.19171 0.21965471 8.6475354e-05 1 + 600 572.3234 0.87879933 8.6475354e-05 1 + 700 996.17398 -0.24269717 0.00010809419 2 + 800 904.50395 1.3662054 0.00010809419 2 + 900 1097.1568 -2.2909907 0.00012971303 3 + 1000 954.08892 1.7705672 0.00012971303 3 + 1100 1102.0377 -1.7018446 0.00015133187 4 + 1200 1239.785 -0.30442903 0.00015133187 4 + 1300 1388.4127 1.3301175 0.00017295071 5 + 1400 1559.3853 1.6709729 0.00017295071 5 + 1500 1471.8623 0.8268427 0.00017295071 5 + 1600 1543.6793 2.1987908 0.00019456955 6 + 1700 1694.5595 0.48852817 0.00019456955 6 + 1800 1632.7737 -1.4617692 0.00021618839 7 + 1900 1922.6502 1.1664257 0.00021618839 7 + 2000 2223.503 -0.95799878 0.00023780722 8 + 2100 2142.6035 0.88444463 0.00025942606 9 + 2200 2298.8636 3.4239313 0.00025942606 9 + 2300 2252.4355 0.82167302 0.00025942606 9 + 2400 2321.0788 1.7499714 0.00025942606 9 + 2500 2095.6715 0.55288444 0.00025942606 9 + 2600 2136.0316 -3.833114 0.00025942606 9 + 2700 2466.3134 -2.2519511 0.00025942606 9 + 2800 2294.3454 1.0637304 0.00025942606 9 + 2900 2340.3891 1.3997049 0.0002810449 10 + 3000 2272.0013 -0.27591886 0.0002810449 10 + 3100 2333.9696 -0.11772138 0.0002810449 10 + 3200 2409.0946 -1.025473 0.0002810449 10 + 3300 2148.023 1.6752329 0.0002810449 10 + 3400 2267.636 -0.45297583 0.0002810449 10 + 3500 2457.622 0.35627297 0.0002810449 10 + 3600 2288.008 -15.516626 0.00030266374 11 + 3700 2458.2681 1.4571773 0.00030266374 11 + 3800 2566.7623 -29.140553 0.00032428258 12 + 3900 2839.4062 0.64583638 0.00032428258 12 + 4000 2893.2204 -53.187892 0.00034590142 13 + 4100 3024.6375 -65.068146 0.00036752025 14 + 4200 3004.6784 1.4155214 0.00036752025 14 + 4300 3033.1895 1.8572273 0.00036752025 14 + 4400 3157.2542 -0.92462977 0.00036752025 14 + 4500 3557.7137 -194.46498 0.00038913909 15 + 4600 3096.485 -1.830492 0.00038913909 15 + 4700 3488.088 -286.81055 0.00041075793 16 + 4800 3390.5493 -372.77818 0.00043237677 17 + 4900 3773.7226 -446.58574 0.00045399561 18 + 5000 3703.0159 -0.81188551 0.00045399561 18 + 5100 4051.3067 1.2567439 0.00045399561 18 + 5200 3813.3682 0.92945737 0.00047561445 19 + 5300 4036.0078 -2.5336258 0.00049723329 20 + 5400 4219.803 -0.96928261 0.00049723329 20 + 5500 4433.7447 -0.026762463 0.00051885212 21 + 5600 4477.4505 -1.417316 0.00054047096 22 + 5700 4500.0306 -1.0551443 0.00054047096 22 + 5800 4600.3507 -4.9580056 0.00054047096 22 + 5900 4765.4978 -2.2546941 0.0005620898 23 + 6000 5442.6193 0.91161284 0.00058370864 24 + 6100 5086.8047 -0.9875332 0.00060532748 25 + 6200 5485.3437 -2.8296626 0.00062694632 26 + 6300 4988.0396 -0.15179023 0.00064856516 27 + 6400 5597.3703 4.2941885 0.00067018399 28 + 6500 5677.0263 -2.8611595 0.00069180283 29 + 6600 6058.0009 1.4111778 0.00071342167 30 + 6700 5859.0817 -2.5782466 0.00071342167 30 + 6800 5879.3941 -4.5681807 0.00071342167 30 + 6900 6398.288 2.5259412 0.00071342167 30 + 7000 6250.1096 -2.6049627 0.00071342167 30 + 7100 5849.651 -0.44062578 0.00071342167 30 + 7200 5778.6532 -0.27299118 0.00071342167 30 + 7300 5977.6661 4.2483639 0.00071342167 30 + 7400 5862.4231 1.0289519 0.00071342167 30 + 7500 6482.376 7.5412373 0.00071342167 30 + 7600 5810.4325 1.0343075 0.00071342167 30 + 7700 5916.7304 2.304302 0.00071342167 30 + 7800 5869.9504 -0.5946555 0.00071342167 30 + 7900 5804.0522 -4.1207689 0.00071342167 30 + 8000 6077.1704 0.52211243 0.00071342167 30 +Loop time of 60.5603 on 4 procs for 8000 steps with 528 atoms + +Performance: 11.413 ns/day, 2.103 hours/ns, 132.100 timesteps/s +99.9% CPU use with 4 MPI tasks x no OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 0.0041695 | 0.90113 | 2.3423 | 102.8 | 1.49 +Bond | 0.011606 | 2.1188 | 5.8107 | 163.9 | 3.50 +Kspace | 47.987 | 52.817 | 55.679 | 43.7 | 87.21 +Neigh | 3.5961 | 3.6262 | 3.6496 | 1.2 | 5.99 +Comm | 0.11097 | 0.16569 | 0.26369 | 15.3 | 0.27 +Output | 0.0020366 | 0.0023427 | 0.0032469 | 1.1 | 0.00 +Modify | 0.62302 | 0.91659 | 1.1227 | 21.5 | 1.51 +Other | | 0.0126 | | | 0.02 + +Nlocal: 132.000 ave 295 max 0 min +Histogram: 2 0 0 0 0 0 0 1 0 1 +Nghost: 133.000 ave 349 max 0 min +Histogram: 2 0 0 0 0 1 0 0 0 1 +Neighs: 8383.50 ave 20143 max 0 min +Histogram: 2 0 0 0 0 0 1 0 0 1 + +Total # of neighbors = 33534 +Ave neighs/atom = 63.511364 +Ave special neighs/atom = 11.409091 +Neighbor list builds = 8000 +Dangerous builds = 0 + +Please see the log.cite file for references relevant to this simulation + +Total wall time: 0:01:00 diff --git a/examples/USER/reaction/create_atoms_polystyrene/trimer.data b/examples/USER/reaction/create_atoms_polystyrene/trimer.data new file mode 100644 index 0000000000..b3ad132f03 --- /dev/null +++ b/examples/USER/reaction/create_atoms_polystyrene/trimer.data @@ -0,0 +1,796 @@ +polystyrene trimer + +48 atoms +7 atom types +50 bonds +13 bond types +84 angles +22 angle types +127 dihedrals +36 dihedral types +36 impropers +9 improper types + +50 250 xlo xhi +50 250 ylo yhi +50 250 zlo zhi + +Masses + +1 12.0112 +2 1.00797 +3 12.0112 +4 12.0112 +5 12.0112 +6 12.0112 +7 12.0112 + +Pair Coeffs # lj/class2/coul/long + +1 0.064 4.01 +2 0.02 2.7 +3 0.064 4.01 +4 0.064 3.9 +5 0.054 4.01 +6 0.054 4.01 +7 0.054 4.01 + +Bond Coeffs # class2 + +1 1.0982 372.825 -803.453 894.317 +2 1.417 470.836 -627.618 1327.63 +3 1.501 321.902 -521.821 572.163 +4 1.0883 365.768 -725.54 781.662 +5 1.34 543.99 -1238.2 1644.03 +6 1.0883 365.768 -725.54 781.662 +7 1.501 321.902 -521.821 572.163 +8 1.101 345 -691.89 844.6 +9 1.53 299.67 -501.77 679.81 +10 1.101 345 -691.89 844.6 +11 1.501 321.902 -521.821 572.163 +12 1.101 345 -691.89 844.6 +13 1.53 299.67 -501.77 679.81 + +Angle Coeffs # class2 + +1 117.94 35.1558 -12.4682 0 +2 118.9 61.0226 -34.9931 0 +3 120.05 44.7148 -22.7352 0 +4 111 44.3234 -9.4454 0 +5 108.4 43.9594 -8.3924 -9.3379 +6 124.88 35.2766 -17.774 -1.6215 +7 124.88 35.2766 -17.774 -1.6215 +8 115.49 29.6363 -12.4853 -6.2218 +9 120.05 44.7148 -22.7352 0 +10 111 44.3234 -9.4454 0 +11 108.4 43.9594 -8.3924 -9.3379 +12 110.77 41.453 -10.604 5.129 +13 112.67 39.516 -7.443 -9.5583 +14 110.77 41.453 -10.604 5.129 +15 107.66 39.641 -12.921 -2.4318 +16 112.67 39.516 -7.443 -9.5583 +17 120.05 44.7148 -22.7352 0 +18 111 44.3234 -9.4454 0 +19 108.4 43.9594 -8.3924 -9.3379 +20 110.77 41.453 -10.604 5.129 +21 110.77 41.453 -10.604 5.129 +22 112.67 39.516 -7.443 -9.5583 + +BondBond Coeffs + +1 1.0795 1.417 1.0982 +2 68.2856 1.417 1.417 +3 12.0676 1.417 1.501 +4 2.9168 1.501 1.0883 +5 0 1.501 1.34 +6 10.1047 1.0883 1.34 +7 10.1047 1.0883 1.34 +8 4.8506 1.0883 1.0883 +9 12.0676 1.417 1.501 +10 2.9168 1.501 1.101 +11 0 1.501 1.53 +12 3.3872 1.101 1.53 +13 0 1.53 1.53 +14 3.3872 1.101 1.53 +15 5.3316 1.101 1.101 +16 0 1.53 1.53 +17 12.0676 1.417 1.501 +18 2.9168 1.501 1.101 +19 0 1.501 1.53 +20 3.3872 1.101 1.53 +21 3.3872 1.101 1.53 +22 0 1.53 1.53 + +BondAngle Coeffs + +1 20.0033 24.2183 1.417 1.0982 +2 28.8708 28.8708 1.417 1.417 +3 31.0771 47.0579 1.417 1.501 +4 26.4608 11.7717 1.501 1.0883 +5 0 0 1.501 1.34 +6 19.0592 23.3588 1.0883 1.34 +7 19.0592 23.3588 1.0883 1.34 +8 17.9795 17.9795 1.0883 1.0883 +9 31.0771 47.0579 1.417 1.501 +10 26.4608 11.7717 1.501 1.101 +11 0 0 1.501 1.53 +12 11.421 20.754 1.101 1.53 +13 8.016 8.016 1.53 1.53 +14 11.421 20.754 1.101 1.53 +15 18.103 18.103 1.101 1.101 +16 8.016 8.016 1.53 1.53 +17 31.0771 47.0579 1.417 1.501 +18 26.4608 11.7717 1.501 1.101 +19 0 0 1.501 1.53 +20 11.421 20.754 1.101 1.53 +21 11.421 20.754 1.101 1.53 +22 8.016 8.016 1.53 1.53 + +Dihedral Coeffs # class2 + +1 0 0 1.559 0 0 0 +2 0 0 3.9661 0 0 0 +3 0 0 4.4072 0 0 0 +4 8.3667 0 1.1932 0 0 0 +5 0 0 1.8769 0 0 0 +6 0 0 0 0 0 0 +7 0 0 0 0 0 0 +8 0 0 0 0 0 0 +9 0 0 4.8974 0 0 0 +10 0 0 1.559 0 0 0 +11 0 0 4.4072 0 0 0 +12 -0.2801 0 -0.0678 0 -0.0122 0 +13 -0.2802 0 -0.0678 0 -0.0122 0 +14 -0.0228 0 0.028 0 -0.1863 0 +15 -0.1432 0 0.0617 0 -0.1083 0 +16 0 0 0.0316 0 -0.1681 0 +17 0 0 0 0 0 0 +18 0 0 0.0316 0 -0.1681 0 +19 0 0 0.0514 0 -0.143 0 +20 0 0 1.559 0 0 0 +21 0 0 4.4072 0 0 0 +22 -0.2801 0 -0.0678 0 -0.0122 0 +23 -0.2802 0 -0.0678 0 -0.0122 0 +24 -0.0228 0 0.028 0 -0.1863 0 +25 0 0 0 0 0 0 +26 -0.1432 0 0.0617 0 -0.1083 0 +27 0 0 0.0316 0 -0.1681 0 +28 0 0 0 0 0 0 +29 0 0 0.0316 0 -0.1681 0 +30 0 0 0.0514 0 -0.143 0 +31 -0.0228 0 0.028 0 -0.1863 0 +32 -0.1432 0 0.0617 0 -0.1083 0 +33 0 0 0.0316 0 -0.1681 0 +34 0 0 0 0 0 0 +35 0 0 0.0316 0 -0.1681 0 +36 0 0 0.0514 0 -0.143 0 + +AngleAngleTorsion Coeffs + +1 4.4444 117.94 120.05 +2 -4.8141 118.9 117.94 +3 -14.4097 118.9 120.05 +4 0 118.9 118.9 +5 0.3598 117.94 117.94 +6 0 120.05 111 +7 0 120.05 108.4 +8 0 108.4 124.88 +9 -7.0058 124.88 124.88 +10 4.4444 117.94 120.05 +11 -14.4097 118.9 120.05 +12 -5.8888 120.05 111 +13 0 120.05 108.4 +14 0 108.4 110.77 +15 -12.564 110.77 110.77 +16 -16.164 112.67 110.77 +17 0 108.4 112.67 +18 -16.164 110.77 112.67 +19 -22.045 112.67 112.67 +20 4.4444 117.94 120.05 +21 -14.4097 118.9 120.05 +22 -5.8888 120.05 111 +23 0 120.05 108.4 +24 0 108.4 110.77 +25 0 108.4 112.67 +26 -12.564 110.77 110.77 +27 -16.164 110.77 112.67 +28 0 112.67 108.4 +29 -16.164 112.67 110.77 +30 -22.045 112.67 112.67 +31 0 110.77 108.4 +32 -12.564 110.77 110.77 +33 -16.164 110.77 112.67 +34 0 112.67 108.4 +35 -16.164 112.67 110.77 +36 -22.045 112.67 112.67 + +EndBondTorsion Coeffs + +1 0 -0.4879 0 0 -1.797 0 1.0982 1.501 +2 0 -6.8958 0 0 -0.4669 0 1.417 1.0982 +3 0 -0.6918 0 0 0.2421 0 1.417 1.501 +4 -0.1185 6.3204 0 -0.1185 6.3204 0 1.417 1.417 +5 0 -0.689 0 0 -0.689 0 1.0982 1.0982 +6 0 0 0 0 0 0 1.417 1.0883 +7 0 0 0 0 0 0 1.417 1.34 +8 0 0 0 0 0 0 1.501 1.0883 +9 0.7129 0.5161 0 0.7129 0.5161 0 1.0883 1.0883 +10 0 -0.4879 0 0 -1.797 0 1.0982 1.501 +11 0 -0.6918 0 0 0.2421 0 1.417 1.501 +12 -0.5835 1.122 0.3978 1.3997 0.7756 0 1.417 1.101 +13 0 0 0 0 0 0 1.417 1.53 +14 0 0 0 0 0 0 1.501 1.101 +15 0.213 0.312 0.0777 0.213 0.312 0.0777 1.101 1.101 +16 0.2486 0.2422 -0.0925 0.0814 0.0591 0.2219 1.53 1.101 +17 0 0 0 0 0 0 1.501 1.53 +18 0.0814 0.0591 0.2219 0.2486 0.2422 -0.0925 1.101 1.53 +19 -0.0732 0 0 -0.0732 0 0 1.53 1.53 +20 0 -0.4879 0 0 -1.797 0 1.0982 1.501 +21 0 -0.6918 0 0 0.2421 0 1.417 1.501 +22 -0.5835 1.122 0.3978 1.3997 0.7756 0 1.417 1.101 +23 0 0 0 0 0 0 1.417 1.53 +24 0 0 0 0 0 0 1.501 1.101 +25 0 0 0 0 0 0 1.501 1.53 +26 0.213 0.312 0.0777 0.213 0.312 0.0777 1.101 1.101 +27 0.0814 0.0591 0.2219 0.2486 0.2422 -0.0925 1.101 1.53 +28 0 0 0 0 0 0 1.53 1.501 +29 0.2486 0.2422 -0.0925 0.0814 0.0591 0.2219 1.53 1.101 +30 -0.0732 0 0 -0.0732 0 0 1.53 1.53 +31 0 0 0 0 0 0 1.101 1.501 +32 0.213 0.312 0.0777 0.213 0.312 0.0777 1.101 1.101 +33 0.0814 0.0591 0.2219 0.2486 0.2422 -0.0925 1.101 1.53 +34 0 0 0 0 0 0 1.53 1.501 +35 0.2486 0.2422 -0.0925 0.0814 0.0591 0.2219 1.53 1.101 +36 -0.0732 0 0 -0.0732 0 0 1.53 1.53 + +MiddleBondTorsion Coeffs + +1 0 3.9421 0 1.417 +2 0 -1.1521 0 1.417 +3 0 9.1792 0 1.417 +4 27.5989 -2.312 0 1.417 +5 0 4.8228 0 1.417 +6 0 0 0 1.501 +7 0 0 0 1.501 +8 0 0 0 1.34 +9 0.8558 6.3911 0 1.34 +10 0 3.9421 0 1.417 +11 0 9.1792 0 1.417 +12 -5.5679 1.4083 0.301 1.501 +13 0 0 0 1.501 +14 0 0 0 1.53 +15 -14.261 -0.5322 -0.4864 1.53 +16 -14.879 -3.6581 -0.3138 1.53 +17 0 0 0 1.53 +18 -14.879 -3.6581 -0.3138 1.53 +19 -17.787 -7.1877 0 1.53 +20 0 3.9421 0 1.417 +21 0 9.1792 0 1.417 +22 -5.5679 1.4083 0.301 1.501 +23 0 0 0 1.501 +24 0 0 0 1.53 +25 0 0 0 1.53 +26 -14.261 -0.5322 -0.4864 1.53 +27 -14.879 -3.6581 -0.3138 1.53 +28 0 0 0 1.53 +29 -14.879 -3.6581 -0.3138 1.53 +30 -17.787 -7.1877 0 1.53 +31 0 0 0 1.53 +32 -14.261 -0.5322 -0.4864 1.53 +33 -14.879 -3.6581 -0.3138 1.53 +34 0 0 0 1.53 +35 -14.879 -3.6581 -0.3138 1.53 +36 -17.787 -7.1877 0 1.53 + +BondBond13 Coeffs + +1 0.8743 1.0982 1.501 +2 -6.2741 1.417 1.0982 +3 2.5085 1.417 1.501 +4 53 1.417 1.417 +5 -1.7077 1.0982 1.0982 +6 0 1.417 1.0883 +7 0 1.417 1.34 +8 0 1.501 1.0883 +9 0 1.0883 1.0883 +10 0.8743 1.0982 1.501 +11 2.5085 1.417 1.501 +12 -3.4826 1.417 1.101 +13 0 1.417 1.53 +14 0 1.501 1.101 +15 0 1.101 1.101 +16 0 1.53 1.101 +17 0 1.501 1.53 +18 0 1.101 1.53 +19 0 1.53 1.53 +20 0.8743 1.0982 1.501 +21 2.5085 1.417 1.501 +22 -3.4826 1.417 1.101 +23 0 1.417 1.53 +24 0 1.501 1.101 +25 0 1.501 1.53 +26 0 1.101 1.101 +27 0 1.101 1.53 +28 0 1.53 1.501 +29 0 1.53 1.101 +30 0 1.53 1.53 +31 0 1.101 1.501 +32 0 1.101 1.101 +33 0 1.101 1.53 +34 0 1.53 1.501 +35 0 1.53 1.101 +36 0 1.53 1.53 + +AngleTorsion Coeffs + +1 0 3.4601 0 0 -0.1242 0 117.94 120.05 +2 0 2.5014 0 0 2.7147 0 118.9 117.94 +3 0 3.8987 0 0 -4.4683 0 118.9 120.05 +4 1.9767 1.0239 0 1.9767 1.0239 0 118.9 118.9 +5 0 2.4501 0 0 2.4501 0 117.94 117.94 +6 0 0 0 0 0 0 120.05 111 +7 0 0 0 0 0 0 120.05 108.4 +8 0 0 0 0 0 0 108.4 124.88 +9 -1.8911 3.254 0 -1.8911 3.254 0 124.88 124.88 +10 0 3.4601 0 0 -0.1242 0 117.94 120.05 +11 0 3.8987 0 0 -4.4683 0 118.9 120.05 +12 0.2251 0.6548 0.1237 4.6266 0.1632 0.0461 120.05 111 +13 0 0 0 0 0 0 120.05 108.4 +14 0 0 0 0 0 0 108.4 110.77 +15 -0.8085 0.5569 -0.2466 -0.8085 0.5569 -0.2466 110.77 110.77 +16 -0.2454 0 -0.1136 0.3113 0.4516 -0.1988 112.67 110.77 +17 0 0 0 0 0 0 108.4 112.67 +18 0.3113 0.4516 -0.1988 -0.2454 0 -0.1136 110.77 112.67 +19 0.3886 -0.3139 0.1389 0.3886 -0.3139 0.1389 112.67 112.67 +20 0 3.4601 0 0 -0.1242 0 117.94 120.05 +21 0 3.8987 0 0 -4.4683 0 118.9 120.05 +22 0.2251 0.6548 0.1237 4.6266 0.1632 0.0461 120.05 111 +23 0 0 0 0 0 0 120.05 108.4 +24 0 0 0 0 0 0 108.4 110.77 +25 0 0 0 0 0 0 108.4 112.67 +26 -0.8085 0.5569 -0.2466 -0.8085 0.5569 -0.2466 110.77 110.77 +27 0.3113 0.4516 -0.1988 -0.2454 0 -0.1136 110.77 112.67 +28 0 0 0 0 0 0 112.67 108.4 +29 -0.2454 0 -0.1136 0.3113 0.4516 -0.1988 112.67 110.77 +30 0.3886 -0.3139 0.1389 0.3886 -0.3139 0.1389 112.67 112.67 +31 0 0 0 0 0 0 110.77 108.4 +32 -0.8085 0.5569 -0.2466 -0.8085 0.5569 -0.2466 110.77 110.77 +33 0.3113 0.4516 -0.1988 -0.2454 0 -0.1136 110.77 112.67 +34 0 0 0 0 0 0 112.67 108.4 +35 -0.2454 0 -0.1136 0.3113 0.4516 -0.1988 112.67 110.77 +36 0.3886 -0.3139 0.1389 0.3886 -0.3139 0.1389 112.67 112.67 + +Improper Coeffs # class2 + +1 4.8912 0 +2 7.8153 0 +3 0 0 +4 2.8561 0 +5 7.8153 0 +6 0 0 +7 0 0 +8 7.8153 0 +9 0 0 + +AngleAngle Coeffs + +1 0 0 0 118.9 117.94 117.94 +2 0 0 0 118.9 120.05 120.05 +3 0 0 0 111 124.88 108.4 +4 0 0 0 115.49 124.88 124.88 +5 0 0 0 118.9 120.05 120.05 +6 0 0 0 107.66 110.77 110.77 +7 0 0 0 111 110.77 108.4 +8 0 0 0 118.9 120.05 120.05 +9 0 0 0 111 110.77 108.4 + +Atoms # full + +44 1 2 3.5400000000000001e-02 6.1476397222913839e+01 8.2376490601205234e+01 6.0906939115836181e+01 +45 1276 2 3.5400000000000001e-02 5.8398688202244472e+01 8.0172948526664996e+01 6.2115536813582672e+01 +46 1276 6 -6.9599999999999995e-02 5.9489073989392523e+01 8.0264057167571636e+01 6.1984002598976552e+01 +48 1276 2 3.5400000000000001e-02 5.9675170230342431e+01 8.0048052449390738e+01 6.0920159395372401e+01 +47 1276 2 1.2370000000000000e-01 5.9297455513100488e+01 8.3187777608476154e+01 5.9645157256520122e+01 +18 1 5 -1.8200000000000001e-02 6.2426251430535707e+01 8.2055473568260709e+01 6.2971661388612958e+01 +19 1 6 -6.9599999999999995e-02 6.1399255844467369e+01 8.1794665295860213e+01 6.1821819828185660e+01 +21 1 1 -1.2900000000000000e-01 6.4032918371445831e+01 8.0190179089286701e+01 6.3021564712316334e+01 +22 1 1 2.6599999999999999e-02 6.3672975135915053e+01 8.1418558650051665e+01 6.2448012627881994e+01 +23 1 2 3.5400000000000001e-02 6.1545198223694939e+01 8.0836309422842305e+01 6.1349823957467130e+01 +27 1276 2 5.1600000000000000e-02 5.9809503696580933e+01 8.1831265916389881e+01 6.3253745193271065e+01 +28 1276 5 -1.8200000000000001e-02 5.9900307947967441e+01 8.1677453781363639e+01 6.2190757403657820e+01 +31 1276 2 1.2370000000000000e-01 5.8050043823867973e+01 8.2698312265456622e+01 6.3667111329534436e+01 +38 1 2 1.2370000000000000e-01 6.3754126973935612e+01 7.9931147303963002e+01 6.4022259163067275e+01 +20 1 2 1.2370000000000000e-01 6.4070158368422781e+01 8.2950071388392274e+01 6.1042631212883315e+01 +24 1 1 -1.2900000000000000e-01 6.4337973861569580e+01 8.1916618276489871e+01 6.1387866780102470e+01 +37 1 2 1.4030000000000001e-01 6.5360115866618415e+01 7.8586112104863830e+01 6.3004997314380716e+01 +39 1 1 -1.7340000000000000e-01 6.5018338085325610e+01 7.9478260591306125e+01 6.2440745569712817e+01 +40 1 1 -1.1340000000000000e-01 6.5628759887796605e+01 7.9941156332165264e+01 6.1248476296558067e+01 +41 1 1 -1.7340000000000000e-01 6.5247995680260402e+01 8.1172439250598345e+01 6.0753045571239831e+01 +42 1 2 1.2880000000000000e-01 6.6569600059599281e+01 7.9514748976494360e+01 6.0810611807135601e+01 +43 1 2 1.4030000000000001e-01 6.5780165393063371e+01 8.1570974991007958e+01 5.9850915261812396e+01 +9 1276 2 1.2880000000000000e-01 5.5651795605743445e+01 8.5074472139235127e+01 6.1094480497979262e+01 +30 1276 2 1.4030000000000001e-01 5.6082982679196888e+01 8.3912863624076010e+01 6.3351889697403472e+01 +33 1276 1 -1.7340000000000000e-01 5.6718133911388506e+01 8.3758479063002000e+01 6.2493293749545209e+01 +34 1276 1 -1.1340000000000000e-01 5.6498352105218459e+01 8.4426576393179090e+01 6.1290147608586011e+01 +6 3822 1 -1.7340000000000000e-01 6.3308103537340351e+01 8.7713509787622499e+01 6.4643082313868433e+01 +7 3822 1 -1.2900000000000000e-01 6.3010291684764312e+01 8.6423650045069493e+01 6.4252844241495922e+01 +8 3822 2 1.2370000000000000e-01 6.2089199187020355e+01 8.6309198636296912e+01 6.3711263099850854e+01 +10 1276 2 1.4030000000000001e-01 5.7266131308654970e+01 8.4599328362003035e+01 5.9281511478144402e+01 +11 3822 2 3.5400000000000001e-02 6.1694306618059791e+01 8.3823470438280594e+01 6.3778953909925114e+01 +12 3822 5 -1.8200000000000001e-02 6.3814926998838651e+01 8.3900077798460728e+01 6.4108991789590448e+01 +13 3822 6 -6.9599999999999995e-02 6.2604540882379787e+01 8.3491998603381077e+01 6.3249610918984622e+01 +14 3822 2 1.2370000000000000e-01 6.5739253131027880e+01 8.4813736128157771e+01 6.5351692111169555e+01 +15 3822 1 -1.2900000000000000e-01 6.5071144269009466e+01 8.5646783550482454e+01 6.5086813218945636e+01 +16 3822 1 2.6599999999999999e-02 6.3957099792282079e+01 8.5375816595044753e+01 6.4385073943729708e+01 +17 1 2 5.1600000000000000e-02 6.2256484483973310e+01 8.1576962161157596e+01 6.3963984654065122e+01 +26 3822 2 5.1600000000000000e-02 6.4196825763126355e+01 8.3291442832977836e+01 6.4907094488854057e+01 +29 1276 1 2.6599999999999999e-02 5.8784742332505303e+01 8.2766055380197670e+01 6.1667239692876961e+01 +32 1276 1 -1.2900000000000000e-01 5.7836199787435064e+01 8.3005060229118428e+01 6.2669788306756018e+01 +35 1276 1 -1.2900000000000000e-01 5.8572661840325132e+01 8.3404075689965083e+01 6.0443288532625175e+01 +36 1276 1 -1.7340000000000000e-01 5.7380616699226330e+01 8.4134680429976896e+01 6.0248710539932475e+01 +25 3822 2 3.5400000000000001e-02 6.2750675036816460e+01 8.3891633300878468e+01 6.2249429178485677e+01 +5 3822 2 1.4030000000000001e-01 6.2626160082050376e+01 8.8416565740835182e+01 6.4093918967496805e+01 +1 3822 2 1.2880000000000000e-01 6.4863557606529355e+01 8.9096029197548390e+01 6.5342927535537825e+01 +2 3822 1 -1.1340000000000000e-01 6.4627442641031166e+01 8.8047381925321190e+01 6.5138073202291650e+01 +3 3822 2 1.4030000000000001e-01 6.6470254992065406e+01 8.6991893750821745e+01 6.5857474890608984e+01 +4 3822 1 -1.7340000000000000e-01 6.5416488888088338e+01 8.6963894801200169e+01 6.5357641085394278e+01 + +Velocities + +44 -1.1274099342391698e-02 2.8614364731871914e-02 7.8116535486555949e-03 +45 2.3164382404151666e-03 3.9815732957733160e-03 -2.9971878581527899e-02 +46 -7.1653099619954563e-03 4.5491360587300133e-04 4.9898614093692017e-03 +48 9.8069086061434527e-03 4.0008139512159270e-03 6.2934259772882122e-03 +47 2.2646445306743783e-03 1.3029071608409702e-03 4.2232440120174040e-02 +18 7.0040064100195757e-03 3.2877451206009701e-03 -3.5376010407568422e-04 +19 -1.3998188760009689e-02 7.2238210565990146e-03 7.7956220633332383e-03 +21 3.1954292320462373e-03 -2.9717583309420764e-03 -3.1753395094325522e-03 +22 5.2997643939121201e-03 -2.9646963088534335e-03 -4.1351926198204894e-03 +23 7.6443400078766528e-03 4.0358953976530103e-02 -2.6684706183248367e-02 +27 1.9261652416455359e-03 -1.1632914130150688e-02 1.0061732021630769e-02 +28 -8.2251676802878315e-03 -1.5111873066969876e-04 1.3808893565582731e-02 +31 5.2475840572179860e-03 1.8266996572138715e-02 2.3453280610166885e-03 +38 -2.0343905130199073e-02 3.2815536859276281e-02 3.6511922534330152e-03 +20 2.2914549087537126e-02 1.4424503744223915e-02 2.1708279654645496e-03 +24 -2.4717233344142471e-03 1.2966123098719246e-02 8.1261459853411936e-03 +37 -2.4547379584186218e-02 -3.0213966592845171e-02 -3.1437442951939183e-02 +39 2.5476117829076835e-03 1.2743160680987653e-03 1.8775880208113892e-03 +40 -6.9216508143939990e-03 1.0986173624795060e-02 8.4543093049661480e-03 +41 -6.9641432145561661e-03 3.4497795547843439e-03 -6.5914679936187716e-03 +42 -1.6682931637687005e-02 -7.9952140358728052e-03 -5.4993265930488526e-02 +43 -1.2747392921213267e-03 -8.9033092043203244e-03 -1.4285400545629027e-02 +9 -4.6235166357676289e-03 -1.3071850427027999e-02 -1.4097407987100977e-02 +30 -1.0949617396609294e-02 2.8255703113196974e-03 1.7171748232322353e-02 +33 -6.1375812469323665e-03 -2.4748644899411924e-03 -9.4761978149296138e-03 +34 1.3676079846441525e-03 5.6076140293943458e-03 4.3217204641336267e-03 +6 -1.0264635053701928e-02 6.5278337056107680e-03 7.0056151148588212e-04 +7 -8.7519451205145676e-03 -4.6476440106580945e-03 2.5970484253527112e-03 +8 2.1377395557690311e-02 -3.3261274153819453e-03 -1.0112266596677577e-02 +10 -3.5793767912309253e-02 -4.7139872292323019e-02 -1.6709528481405608e-02 +11 8.5071485795589590e-03 9.9402848610678270e-03 -3.8088596341056854e-03 +12 -7.1678159384257103e-04 -6.9164463557228907e-04 -6.4073519808107186e-03 +13 -4.8443902657902991e-03 -1.1919190682985097e-03 6.3946846087726637e-03 +14 1.4810157483257907e-02 1.9829623839419017e-03 -2.7393844990063056e-02 +15 2.4171850935506777e-03 8.5003135180758520e-03 -1.4373227798951704e-03 +16 2.7567342910947553e-03 4.7168484476890456e-03 -5.5131873288712992e-03 +17 -3.8456662730386774e-02 2.0220106671151108e-02 -1.3822049134399602e-02 +26 2.7415414728694614e-02 1.4392155257037418e-03 -6.7281635499082748e-03 +29 2.8284983560440745e-03 2.8809942505517976e-03 -9.0489583066552114e-04 +32 -3.8543634697614316e-03 4.6751647301899795e-03 4.2171867397204537e-03 +35 -8.6957974827209118e-03 -4.4615282666186267e-04 -2.6571026120482824e-03 +36 9.4881057996863086e-04 -7.5665878069688429e-03 2.0333670960646154e-03 +25 1.8105924111310519e-02 -8.6933495274689535e-03 -1.9695291360338044e-04 +5 -5.0447438383189585e-03 -4.5665146331657552e-02 1.0653751333175230e-02 +1 -1.7372868398038824e-02 -2.3625357536259349e-03 1.2220266128368908e-02 +2 3.7050246021929395e-03 -1.0236943515935205e-03 7.2206774682170580e-03 +3 2.3669435799326944e-02 2.7891427939155597e-02 -6.7091036888174346e-03 +4 3.4910623999263577e-03 2.6370880132825258e-03 -6.4694788112864129e-03 + +Bonds + +1 10 44 19 +2 10 45 46 +3 10 48 46 +4 9 19 18 +5 1 21 38 +6 2 21 22 +7 2 21 39 +8 7 22 18 +9 2 22 24 +10 10 23 19 +11 8 27 28 +12 9 28 46 +13 9 28 19 +14 1 24 20 +15 2 24 41 +16 1 39 37 +17 1 40 42 +18 2 40 39 +19 1 41 43 +20 2 41 40 +21 1 33 30 +22 1 34 9 +23 2 34 33 +24 1 6 5 +25 2 6 2 +26 1 7 8 +27 2 7 6 +28 10 11 13 +29 13 12 13 +30 9 13 18 +31 1 15 14 +32 2 15 16 +33 2 15 4 +34 11 16 12 +35 2 16 7 +36 8 17 18 +37 12 26 12 +38 7 29 28 +39 2 29 35 +40 1 32 31 +41 2 32 29 +42 2 32 33 +43 1 35 47 +44 2 35 36 +45 1 36 10 +46 2 36 34 +47 10 25 13 +48 1 2 1 +49 2 2 4 +50 1 4 3 + +Angles + +1 14 45 46 28 +2 14 48 46 28 +3 15 45 46 48 +4 11 22 18 13 +5 12 17 18 13 +6 13 13 18 19 +7 10 22 18 17 +8 11 22 18 19 +9 12 17 18 19 +10 16 28 19 18 +11 14 44 19 28 +12 14 23 19 28 +13 14 44 19 18 +14 14 23 19 18 +15 15 44 19 23 +16 1 22 21 38 +17 1 39 21 38 +18 2 22 21 39 +19 9 21 22 18 +20 2 21 22 24 +21 9 24 22 18 +22 10 29 28 27 +23 11 29 28 46 +24 11 29 28 19 +25 12 27 28 46 +26 12 27 28 19 +27 13 46 28 19 +28 1 22 24 20 +29 2 22 24 41 +30 1 41 24 20 +31 2 21 39 40 +32 1 21 39 37 +33 1 40 39 37 +34 1 41 40 42 +35 2 41 40 39 +36 1 39 40 42 +37 1 24 41 43 +38 2 24 41 40 +39 1 40 41 43 +40 2 32 33 34 +41 1 32 33 30 +42 1 34 33 30 +43 1 36 34 9 +44 2 36 34 33 +45 1 33 34 9 +46 1 7 6 5 +47 2 7 6 2 +48 1 2 6 5 +49 1 16 7 8 +50 2 16 7 6 +51 1 6 7 8 +52 18 16 12 26 +53 19 16 12 13 +54 20 26 12 13 +55 21 25 13 12 +56 21 11 13 12 +57 22 12 13 18 +58 15 25 13 11 +59 14 25 13 18 +60 14 11 13 18 +61 1 16 15 14 +62 1 4 15 14 +63 2 16 15 4 +64 17 15 16 12 +65 2 15 16 7 +66 17 7 16 12 +67 9 32 29 28 +68 2 32 29 35 +69 9 35 29 28 +70 1 29 32 31 +71 1 33 32 31 +72 2 29 32 33 +73 1 29 35 47 +74 2 29 35 36 +75 1 36 35 47 +76 1 35 36 10 +77 2 35 36 34 +78 1 34 36 10 +79 1 6 2 1 +80 2 6 2 4 +81 1 4 2 1 +82 2 15 4 2 +83 1 15 4 3 +84 1 2 4 3 + +Dihedrals + +1 34 18 19 28 29 +2 31 44 19 28 29 +3 31 23 19 28 29 +4 35 18 19 28 27 +5 32 44 19 28 27 +6 32 23 19 28 27 +7 36 18 19 28 46 +8 33 44 19 28 46 +9 33 23 19 28 46 +10 36 28 19 18 13 +11 33 44 19 18 13 +12 33 23 19 18 13 +13 34 28 19 18 22 +14 31 44 19 18 22 +15 31 23 19 18 22 +16 35 28 19 18 17 +17 32 44 19 18 17 +18 32 23 19 18 17 +19 10 38 21 22 18 +20 11 39 21 22 18 +21 4 39 21 22 24 +22 5 38 21 39 37 +23 4 22 21 39 40 +24 2 22 21 39 37 +25 2 24 22 21 38 +26 13 21 22 18 13 +27 12 21 22 18 17 +28 13 21 22 18 19 +29 13 24 22 18 13 +30 12 24 22 18 17 +31 13 24 22 18 19 +32 2 21 22 24 20 +33 4 21 22 24 41 +34 14 29 28 46 45 +35 14 29 28 46 48 +36 15 27 28 46 45 +37 15 27 28 46 48 +38 16 19 28 46 45 +39 16 19 28 46 48 +40 10 20 24 22 18 +41 11 41 24 22 18 +42 2 22 24 41 43 +43 4 22 24 41 40 +44 5 20 24 41 43 +45 2 40 39 21 38 +46 2 21 39 40 42 +47 2 39 40 41 43 +48 4 41 40 39 21 +49 2 41 40 39 37 +50 5 42 40 39 37 +51 2 40 41 24 20 +52 2 24 41 40 42 +53 4 24 41 40 39 +54 5 43 41 40 42 +55 2 34 33 32 31 +56 2 32 33 34 9 +57 2 33 34 36 10 +58 4 36 34 33 32 +59 2 36 34 33 30 +60 5 9 34 33 30 +61 2 2 6 7 8 +62 2 7 6 2 1 +63 4 7 6 2 4 +64 5 5 6 2 1 +65 20 8 7 16 12 +66 21 6 7 16 12 +67 2 16 7 6 5 +68 4 16 7 6 2 +69 5 8 7 6 5 +70 24 16 12 13 25 +71 24 16 12 13 11 +72 25 16 12 13 18 +73 26 26 12 13 25 +74 26 26 12 13 11 +75 27 26 12 13 18 +76 28 12 13 18 22 +77 29 12 13 18 17 +78 30 12 13 18 19 +79 31 25 13 18 22 +80 32 25 13 18 17 +81 33 25 13 18 19 +82 31 11 13 18 22 +83 32 11 13 18 17 +84 33 11 13 18 19 +85 20 14 15 16 12 +86 21 4 15 16 12 +87 4 4 15 16 7 +88 5 14 15 4 3 +89 4 16 15 4 2 +90 2 16 15 4 3 +91 2 7 16 15 14 +92 22 15 16 12 26 +93 23 15 16 12 13 +94 22 7 16 12 26 +95 23 7 16 12 13 +96 2 15 16 7 8 +97 4 15 16 7 6 +98 2 35 29 32 31 +99 12 32 29 28 27 +100 13 32 29 28 46 +101 13 32 29 28 19 +102 12 35 29 28 27 +103 13 35 29 28 46 +104 13 35 29 28 19 +105 2 32 29 35 47 +106 4 32 29 35 36 +107 10 31 32 29 28 +108 11 33 32 29 28 +109 4 33 32 29 35 +110 5 31 32 33 30 +111 4 29 32 33 34 +112 2 29 32 33 30 +113 10 47 35 29 28 +114 11 36 35 29 28 +115 2 29 35 36 10 +116 4 29 35 36 34 +117 5 47 35 36 10 +118 2 34 36 35 47 +119 2 35 36 34 9 +120 4 35 36 34 33 +121 5 10 36 34 9 +122 2 4 2 6 5 +123 4 6 2 4 15 +124 2 6 2 4 3 +125 5 1 2 4 3 +126 2 2 4 15 14 +127 2 15 4 2 1 + +Impropers + +1 6 45 46 48 28 +2 1 22 18 17 13 +3 1 22 18 13 19 +4 1 17 18 13 19 +5 1 22 18 17 19 +6 1 44 19 18 28 +7 1 23 19 18 28 +8 1 44 19 23 28 +9 1 44 19 23 18 +10 1 22 21 39 38 +11 5 21 22 24 18 +12 1 29 28 27 46 +13 1 29 28 27 19 +14 1 29 28 46 19 +15 1 27 28 46 19 +16 1 22 24 41 20 +17 1 21 39 40 37 +18 1 41 40 39 42 +19 1 24 41 40 43 +20 1 32 33 34 30 +21 1 36 34 33 9 +22 1 7 6 2 5 +23 1 16 7 6 8 +24 9 16 12 26 13 +25 1 25 13 11 12 +26 1 25 13 12 18 +27 1 11 13 12 18 +28 1 25 13 11 18 +29 1 16 15 4 14 +30 8 15 16 7 12 +31 5 32 29 35 28 +32 1 29 32 33 31 +33 1 29 35 36 47 +34 1 35 36 34 10 +35 1 6 2 4 1 +36 1 15 4 2 3 diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp index 86e7977d9b..657caf1e68 100644 --- a/src/USER-REACTION/fix_bond_react.cpp +++ b/src/USER-REACTION/fix_bond_react.cpp @@ -88,6 +88,9 @@ enum{ATOM,FRAG}; // keyword values that accept variables as input enum{NEVERY,RMIN,RMAX,PROB}; +// flag for one-proc vs shared reaction sites +enum{LOCAL,GLOBAL}; + // values for molecule_keyword enum{OFF,INTER,INTRA}; @@ -208,6 +211,9 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : memory->create(limit_duration,nreacts,"bond/react:limit_duration"); memory->create(stabilize_steps_flag,nreacts,"bond/react:stabilize_steps_flag"); memory->create(custom_charges_fragid,nreacts,"bond/react:custom_charges_fragid"); + memory->create(create_atoms_flag,nreacts,"bond/react:create_atoms_flag"); + memory->create(modify_create_fragid,nreacts,"bond/react:modify_create_fragid"); + memory->create(overlapsq,nreacts,"bond/react:overlapsq"); memory->create(molecule_keyword,nreacts,"bond/react:molecule_keyword"); memory->create(nconstraints,nreacts,"bond/react:nconstraints"); memory->create(constraintstr,nreacts,MAXLINE,"bond/react:constraintstr"); @@ -231,6 +237,9 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : max_rxn[i] = INT_MAX; stabilize_steps_flag[i] = 0; custom_charges_fragid[i] = -1; + create_atoms_flag[i] = 0; + modify_create_fragid[i] = -1; + overlapsq[i] = 0; molecule_keyword[i] = OFF; nconstraints[i] = 0; // set default limit duration to 60 timesteps @@ -396,6 +405,28 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : else if (strcmp(arg[iarg+1],"intra") == 0) molecule_keyword[rxn] = INTRA; else error->one(FLERR,"Bond/react: Illegal option for 'molecule' keyword"); iarg += 2; + } else if (strcmp(arg[iarg],"modify_create") == 0) { + if (iarg++ > narg) error->all(FLERR,"Illegal fix bond/react command: " + "'modify_create' has too few arguments"); + while (iarg < narg && strcmp(arg[iarg],"react") != 0 ) { + if (strcmp(arg[iarg],"fit") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/react command: " + "'modify_create' has too few arguments"); + if (strcmp(arg[iarg+1],"all") == 0) modify_create_fragid[rxn] = -1; //default + else { + modify_create_fragid[rxn] = atom->molecules[reacted_mol[rxn]]->findfragment(arg[iarg+1]); + if (modify_create_fragid[rxn] < 0) error->one(FLERR,"Bond/react: Molecule fragment for " + "'modify_create' keyword does not exist"); + } + iarg += 2; + } else if (strcmp(arg[iarg],"overlap") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/react command: " + "'modify_create' has too few arguments"); + overlapsq[rxn] = utils::numeric(FLERR,arg[iarg+1],false,lmp); + overlapsq[rxn] *= overlapsq[rxn]; + iarg += 2; + } else break; + } } else error->all(FLERR,"Illegal fix bond/react command: unknown keyword"); } } @@ -412,6 +443,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : memory->create(landlocked_atoms,max_natoms,nreacts,"bond/react:landlocked_atoms"); memory->create(custom_charges,max_natoms,nreacts,"bond/react:custom_charges"); memory->create(delete_atoms,max_natoms,nreacts,"bond/react:delete_atoms"); + memory->create(create_atoms,max_natoms,nreacts,"bond/react:create_atoms"); memory->create(chiral_atoms,max_natoms,6,nreacts,"bond/react:chiral_atoms"); for (int j = 0; j < nreacts; j++) @@ -419,9 +451,15 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : edge[i][j] = 0; custom_charges[i][j] = 1; // update all partial charges by default delete_atoms[i][j] = 0; + create_atoms[i][j] = 0; for (int k = 0; k < 6; k++) { chiral_atoms[i][k][j] = 0; } + // default equivalences to their own mol index + // all but created atoms will be updated + for (int m = 0; m < 2; m++) { + equivalences[i][m][j] = i+1; + } } // read all map files afterward @@ -431,11 +469,13 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : twomol = atom->molecules[reacted_mol[i]]; onemol->check_attributes(0); twomol->check_attributes(0); - if (onemol->natoms != twomol->natoms) - error->all(FLERR,"Bond/react: Reaction templates must contain the same number of atoms"); get_molxspecials(); read(i); fclose(fp); + if (ncreate == 0 && onemol->natoms != twomol->natoms) + error->all(FLERR,"Bond/react: Reaction templates must contain the same number of atoms"); + else if (ncreate > 0 && onemol->natoms + ncreate != twomol->natoms) + error->all(FLERR,"Bond/react: Incorrect number of created atoms"); iatomtype[i] = onemol->type[ibonding[i]-1]; jatomtype[i] = onemol->type[jbonding[i]-1]; find_landlocked_atoms(i); @@ -498,10 +538,10 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : partner = finalpartner = nullptr; distsq = nullptr; probability = nullptr; - maxcreate = 0; - created = nullptr; - ncreate = nullptr; - allncreate = 0; + maxattempt = 0; + attempt = nullptr; + nattempt = nullptr; + allnattempt = 0; local_num_mega = 0; ghostly_num_mega = 0; restore = nullptr; @@ -543,18 +583,20 @@ FixBondReact::~FixBondReact() memory->destroy(partner); memory->destroy(finalpartner); - memory->destroy(ncreate); + memory->destroy(nattempt); memory->destroy(distsq); memory->destroy(probability); - memory->destroy(created); + memory->destroy(attempt); memory->destroy(edge); memory->destroy(equivalences); memory->destroy(reverse_equiv); memory->destroy(landlocked_atoms); memory->destroy(custom_charges); memory->destroy(delete_atoms); + memory->destroy(create_atoms); memory->destroy(chiral_atoms); + memory->destroy(rxn_name); memory->destroy(nevery); memory->destroy(cutsq); memory->destroy(unreacted_mol); @@ -572,7 +614,10 @@ FixBondReact::~FixBondReact() memory->destroy(molecule_keyword); memory->destroy(constraints); memory->destroy(nconstraints); - // need to delete rxn_name and constraintstr + memory->destroy(constraintstr); + memory->destroy(create_atoms_flag); + memory->destroy(modify_create_fragid); + memory->destroy(overlapsq); memory->destroy(iatomtype); memory->destroy(jatomtype); @@ -777,7 +822,7 @@ void FixBondReact::init() // check cutoff for iatomtype,jatomtype for (int i = 0; i < nreacts; i++) { - if (closeneigh[i] == -1) // indicates will search for non-bonded bonding atoms + if (!utils::strmatch(force->pair_style,"^hybrid")) if (force->pair == nullptr || cutsq[i][1] > force->pair->cutsq[iatomtype[i]][jatomtype[i]]) error->all(FLERR,"Bond/react: Fix bond/react cutoff is longer than pairwise cutoff"); } @@ -844,19 +889,19 @@ void FixBondReact::post_integrate() memory->destroy(partner); memory->destroy(finalpartner); memory->destroy(distsq); - memory->destroy(ncreate); + memory->destroy(nattempt); memory->destroy(probability); nmax = atom->nmax; memory->create(partner,nmax,"bond/react:partner"); memory->create(finalpartner,nmax,"bond/react:finalpartner"); memory->create(distsq,nmax,2,"bond/react:distsq"); - memory->create(ncreate,nreacts,"bond/react:ncreate"); + memory->create(nattempt,nreacts,"bond/react:nattempt"); memory->create(probability,nmax,"bond/react:probability"); } // reset create counts for (int i = 0; i < nreacts; i++) { - ncreate[i] = 0; + nattempt[i] = 0; } int nlocal = atom->nlocal; @@ -938,7 +983,7 @@ void FixBondReact::post_integrate() // and probability constraint is satisfied // if other atom is owned by another proc, it should do same thing - int temp_ncreate = 0; + int temp_nattempt = 0; for (int i = 0; i < nlocal; i++) { if (partner[i] == 0) { continue; @@ -959,17 +1004,17 @@ void FixBondReact::post_integrate() } } - // store final created bond partners and count the rxn possibility once + // store final bond partners and count the rxn possibility once finalpartner[i] = tag[j]; finalpartner[j] = tag[i]; - if (tag[i] < tag[j]) temp_ncreate++; + if (tag[i] < tag[j]) temp_nattempt++; } // cycle 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); + MPI_Allreduce(&temp_nattempt,&some_chance,1,MPI_INT,MPI_SUM,world); if (!some_chance) continue; // communicate final partner @@ -977,7 +1022,7 @@ void FixBondReact::post_integrate() commflag = 3; comm->forward_comm_fix(this); - // add instance to 'created' only if this processor + // add instance to 'attempt' 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) @@ -988,21 +1033,21 @@ void FixBondReact::post_integrate() j = atom->map(finalpartner[i]); // if (j < 0 || tag[i] < tag[j]) { if (tag[i] < tag[j]) { //atom->map(std::min(tag[i],tag[j])) <= nlocal && - if (ncreate[rxnID] == maxcreate) { - maxcreate += DELTA; - // third column of 'created': bond/react integer ID - memory->grow(created,maxcreate,2,nreacts,"bond/react:created"); + if (nattempt[rxnID] == maxattempt) { + maxattempt += DELTA; + // third column of 'attempt': bond/react integer ID + memory->grow(attempt,maxattempt,2,nreacts,"bond/react:attempt"); } // to ensure types remain in same order // unnecessary now taken from reaction map file if (iatomtype[rxnID] == type[i]) { - created[ncreate[rxnID]][0][rxnID] = tag[i]; - created[ncreate[rxnID]][1][rxnID] = finalpartner[i]; + attempt[nattempt[rxnID]][0][rxnID] = tag[i]; + attempt[nattempt[rxnID]][1][rxnID] = finalpartner[i]; } else { - created[ncreate[rxnID]][0][rxnID] = finalpartner[i]; - created[ncreate[rxnID]][1][rxnID] = tag[i]; + attempt[nattempt[rxnID]][0][rxnID] = finalpartner[i]; + attempt[nattempt[rxnID]][1][rxnID] = tag[i]; } - ncreate[rxnID]++; + nattempt[rxnID]++; } } } @@ -1010,11 +1055,11 @@ void FixBondReact::post_integrate() // break loop if no even eligible bonding atoms were found (on any proc) int some_chance; - allncreate = 0; + allnattempt = 0; for (int i = 0; i < nreacts; i++) - allncreate += ncreate[i]; + allnattempt += nattempt[i]; - MPI_Allreduce(&allncreate,&some_chance,1,MPI_INT,MPI_SUM,world); + MPI_Allreduce(&allnattempt,&some_chance,1,MPI_INT,MPI_SUM,world); if (!some_chance) { unlimit_bond(); return; @@ -1242,13 +1287,13 @@ void FixBondReact::superimpose_algorithm() memory->create(restore_pt,MAXGUESS,4,"bond/react:restore_pt"); memory->create(pioneers,max_natoms,"bond/react:pioneers"); memory->create(restore,max_natoms,MAXGUESS*4,"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"); + memory->create(local_mega_glove,max_natoms+1,allnattempt,"bond/react:local_mega_glove"); + memory->create(ghostly_mega_glove,max_natoms+1,allnattempt,"bond/react:ghostly_mega_glove"); attempted_rxn = 1; for (int i = 0; i < max_natoms+1; i++) { - for (int j = 0; j < allncreate; j++) { + for (int j = 0; j < allnattempt; j++) { local_mega_glove[i][j] = 0; ghostly_mega_glove[i][j] = 0; } @@ -1256,7 +1301,7 @@ void FixBondReact::superimpose_algorithm() // let's finally begin the superimpose loop for (rxnID = 0; rxnID < nreacts; rxnID++) { - for (lcl_inst = 0; lcl_inst < ncreate[rxnID]; lcl_inst++) { + for (lcl_inst = 0; lcl_inst < nattempt[rxnID]; lcl_inst++) { onemol = atom->molecules[unreacted_mol[rxnID]]; twomol = atom->molecules[reacted_mol[rxnID]]; @@ -1279,10 +1324,10 @@ void FixBondReact::superimpose_algorithm() int myjbonding = jbonding[rxnID]; glove[myibonding-1][0] = myibonding; - glove[myibonding-1][1] = created[lcl_inst][0][rxnID]; + glove[myibonding-1][1] = attempt[lcl_inst][0][rxnID]; glove_counter++; glove[myjbonding-1][0] = myjbonding; - glove[myjbonding-1][1] = created[lcl_inst][1][rxnID]; + glove[myjbonding-1][1] = attempt[lcl_inst][1][rxnID]; glove_counter++; // special case, only two atoms in reaction templates @@ -1353,7 +1398,7 @@ void FixBondReact::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 + dedup_mega_gloves(LOCAL); // 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); @@ -1405,8 +1450,8 @@ void FixBondReact::superimpose_algorithm() 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); + limit_bond(LOCAL); // add reacting atoms to nve/limit + limit_bond(GLOBAL); update_everything(); // change topology } @@ -1912,7 +1957,7 @@ int FixBondReact::check_constraints() } if (ANDgate != 1) satisfied[i] = 0; } else if (constraints[i][rxnID].type == ARRHENIUS) { - t = get_temperature(); + t = get_temperature(glove,0,1); prrhob = constraints[i][rxnID].par[1]*pow(t,constraints[i][rxnID].par[2])* exp(-constraints[i][rxnID].par[3]/(force->boltz*t)); if (prrhob < rrhandom[(int) constraints[i][rxnID].par[0]]->uniform()) satisfied[i] = 0; @@ -2044,7 +2089,7 @@ void FixBondReact::get_IDcoords(int mode, int myID, double *center) compute local temperature: average over all atoms in reaction template ------------------------------------------------------------------------- */ -double FixBondReact::get_temperature() +double FixBondReact::get_temperature(tagint **myglove, int row_offset, int col) { int i,ilocal; double adof = domain->dimension; @@ -2058,13 +2103,13 @@ double FixBondReact::get_temperature() if (rmass) { for (i = 0; i < onemol->natoms; i++) { - ilocal = atom->map(glove[i][1]); + ilocal = atom->map(myglove[i+row_offset][col]); t += (v[ilocal][0]*v[ilocal][0] + v[ilocal][1]*v[ilocal][1] + v[ilocal][2]*v[ilocal][2]) * rmass[ilocal]; } } else { for (i = 0; i < onemol->natoms; i++) { - ilocal = atom->map(glove[i][1]); + ilocal = atom->map(myglove[i+row_offset][col]); t += (v[ilocal][0]*v[ilocal][0] + v[ilocal][1]*v[ilocal][1] + v[ilocal][2]*v[ilocal][2]) * mass[type[ilocal]]; } @@ -2169,7 +2214,8 @@ void FixBondReact::find_landlocked_atoms(int myrxn) // always remove edge atoms from landlocked list for (int i = 0; i < twomol->natoms; i++) { - if (edge[equivalences[i][1][myrxn]-1][myrxn] == 1) landlocked_atoms[i][myrxn] = 0; + if (create_atoms[i][myrxn] == 0 && edge[equivalences[i][1][myrxn]-1][myrxn] == 1) + landlocked_atoms[i][myrxn] = 0; else landlocked_atoms[i][myrxn] = 1; } int nspecial_limit = -1; @@ -2193,44 +2239,48 @@ void FixBondReact::find_landlocked_atoms(int myrxn) // bad molecule templates check // if atoms change types, but aren't landlocked, that's bad for (int i = 0; i < twomol->natoms; i++) { - if (twomol->type[i] != onemol->type[equivalences[i][1][myrxn]-1] && landlocked_atoms[i][myrxn] == 0) { - char str[128]; - snprintf(str,128,"Bond/react: Atom type affected by reaction %s too close to template edge",rxn_name[myrxn]); - error->all(FLERR,str); + if (create_atoms[i][myrxn] == 0) { + if (twomol->type[i] != onemol->type[equivalences[i][1][myrxn]-1] && landlocked_atoms[i][myrxn] == 0) { + char str[128]; + snprintf(str,128,"Bond/react: Atom type affected by reaction %s too close to template edge",rxn_name[myrxn]); + error->all(FLERR,str); + } } } // additionally, if a bond changes type, but neither involved atom is landlocked, bad // would someone want to change an angle type but not bond or atom types? (etc.) ...hopefully not yet for (int i = 0; i < twomol->natoms; i++) { - if (landlocked_atoms[i][myrxn] == 0) { - for (int j = 0; j < twomol->num_bond[i]; j++) { - int twomol_atomj = twomol->bond_atom[i][j]; - if (landlocked_atoms[twomol_atomj-1][myrxn] == 0) { - int onemol_atomi = equivalences[i][1][myrxn]; - int onemol_batom; - for (int m = 0; m < onemol->num_bond[onemol_atomi-1]; m++) { - onemol_batom = onemol->bond_atom[onemol_atomi-1][m]; - if (onemol_batom == equivalences[twomol_atomj-1][1][myrxn]) { - if (twomol->bond_type[i][j] != onemol->bond_type[onemol_atomi-1][m]) { - char str[128]; - snprintf(str,128,"Bond/react: Bond type affected by reaction %s too close to template edge",rxn_name[myrxn]); - error->all(FLERR,str); - } - } - } - if (newton_bond) { - int onemol_atomj = equivalences[twomol_atomj-1][1][myrxn]; - for (int m = 0; m < onemol->num_bond[onemol_atomj-1]; m++) { - onemol_batom = onemol->bond_atom[onemol_atomj-1][m]; - if (onemol_batom == equivalences[i][1][myrxn]) { - if (twomol->bond_type[i][j] != onemol->bond_type[onemol_atomj-1][m]) { + if (create_atoms[i][myrxn] == 0) { + if (landlocked_atoms[i][myrxn] == 0) { + for (int j = 0; j < twomol->num_bond[i]; j++) { + int twomol_atomj = twomol->bond_atom[i][j]; + if (landlocked_atoms[twomol_atomj-1][myrxn] == 0) { + int onemol_atomi = equivalences[i][1][myrxn]; + int onemol_batom; + for (int m = 0; m < onemol->num_bond[onemol_atomi-1]; m++) { + onemol_batom = onemol->bond_atom[onemol_atomi-1][m]; + if (onemol_batom == equivalences[twomol_atomj-1][1][myrxn]) { + if (twomol->bond_type[i][j] != onemol->bond_type[onemol_atomi-1][m]) { char str[128]; snprintf(str,128,"Bond/react: Bond type affected by reaction %s too close to template edge",rxn_name[myrxn]); error->all(FLERR,str); } } } + if (newton_bond) { + int onemol_atomj = equivalences[twomol_atomj-1][1][myrxn]; + for (int m = 0; m < onemol->num_bond[onemol_atomj-1]; m++) { + onemol_batom = onemol->bond_atom[onemol_atomj-1][m]; + if (onemol_batom == equivalences[i][1][myrxn]) { + if (twomol->bond_type[i][j] != onemol->bond_type[onemol_atomj-1][m]) { + char str[128]; + snprintf(str,128,"Bond/react: Bond type affected by reaction %s too close to template edge",rxn_name[myrxn]); + error->all(FLERR,str); + } + } + } + } } } } @@ -2252,13 +2302,22 @@ void FixBondReact::find_landlocked_atoms(int myrxn) // also, if atoms change number of bonds, but aren't landlocked, that could be bad if (me == 0) for (int i = 0; i < twomol->natoms; i++) { - if (twomol_nxspecial[i][0] != onemol_nxspecial[equivalences[i][1][myrxn]-1][0] && landlocked_atoms[i][myrxn] == 0) { - char str[128]; - snprintf(str,128,"Bond/react: Atom affected by reaction %s too close to template edge",rxn_name[myrxn]); - error->warning(FLERR,str); - break; + if (create_atoms[i][myrxn] == 0) { + if (twomol_nxspecial[i][0] != onemol_nxspecial[equivalences[i][1][myrxn]-1][0] && landlocked_atoms[i][myrxn] == 0) { + char str[128]; + snprintf(str,128,"Bond/react: Atom affected by reaction %s too close to template edge",rxn_name[myrxn]); + error->warning(FLERR,str); + break; + } } } + + // finally, if a created atom is not landlocked, bad! + for (int i = 0; i < twomol->natoms; i++) { + if (create_atoms[i][myrxn] == 1 && landlocked_atoms[i][myrxn] == 0) { + error->one(FLERR,"Bond/react: Created atom too close to template edge"); + } + } } /* ---------------------------------------------------------------------- @@ -2268,30 +2327,30 @@ 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 + // dedup_mode == LOCAL for local_dedup + // dedup_mode == GLOBAL 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; + if (dedup_mode == LOCAL) local_rxn_count[i] = 0; + if (dedup_mode == GLOBAL) ghostly_rxn_count[i] = 0; } int dedup_size = 0; - if (dedup_mode == 0) { + if (dedup_mode == LOCAL) { dedup_size = local_num_mega; - } else if (dedup_mode == 1) { + } else if (dedup_mode == GLOBAL) { dedup_size = global_megasize; } tagint **dedup_glove; memory->create(dedup_glove,max_natoms+1,dedup_size,"bond/react:dedup_glove"); - if (dedup_mode == 0) { + if (dedup_mode == LOCAL) { 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) { + } else if (dedup_mode == GLOBAL) { 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]; @@ -2371,7 +2430,7 @@ void FixBondReact::dedup_mega_gloves(int dedup_mode) // we must update local_mega_glove and local_megasize // we can simply overwrite local_mega_glove column by column - if (dedup_mode == 0) { + if (dedup_mode == LOCAL) { int new_local_megasize = 0; for (int i = 0; i < local_num_mega; i++) { if (dedup_mask[i] == 0) { @@ -2388,7 +2447,7 @@ void FixBondReact::dedup_mega_gloves(int dedup_mode) // we must update global_mega_glove and global_megasize // we can simply overwrite global_mega_glove column by column - if (dedup_mode == 1) { + if (dedup_mode == GLOBAL) { int new_global_megasize = 0; for (int i = 0; i < global_megasize; i++) { if (dedup_mask[i] == 0) { @@ -2425,7 +2484,7 @@ void FixBondReact::limit_bond(int limit_bond_mode) int nlocal = atom->nlocal; int temp_limit_num = 0; tagint *temp_limit_glove; - if (limit_bond_mode == 0) { + if (limit_bond_mode == LOCAL) { 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++) { @@ -2436,7 +2495,7 @@ void FixBondReact::limit_bond(int limit_bond_mode) } } - } else if (limit_bond_mode == 1) { + } else if (limit_bond_mode == GLOBAL) { int max_temp = global_megasize * (max_natoms + 1); temp_limit_glove = new tagint[max_temp]; for (int j = 0; j < global_megasize; j++) { @@ -2532,11 +2591,15 @@ void FixBondReact::glove_ghostcheck() int ghostly = 0; #if !defined(MPI_STUBS) if (comm->style == 0) { - for (int i = 0; i < onemol->natoms; i++) { - int ilocal = atom->map(glove[i][1]); - if (ilocal >= atom->nlocal || localsendlist[ilocal] == 1) { - ghostly = 1; - break; + if (create_atoms_flag[rxnID] == 1) { + ghostly = 1; + } else { + for (int i = 0; i < onemol->natoms; i++) { + int ilocal = atom->map(glove[i][1]); + if (ilocal >= atom->nlocal || localsendlist[ilocal] == 1) { + ghostly = 1; + break; + } } } } else { @@ -2624,7 +2687,7 @@ void FixBondReact::ghost_glovecast() column, 0, world); } - if (me == 0) dedup_mega_gloves(1); // global_mega_glove mode + if (me == 0) dedup_mega_gloves(GLOBAL); // 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); @@ -2642,9 +2705,8 @@ update molecule IDs, charges, types, special lists and all topology void FixBondReact::update_everything() { + int nlocal; // must be defined after create_atoms int *type = atom->type; - - int nlocal = atom->nlocal; int **nspecial = atom->nspecial; tagint **special = atom->special; @@ -2654,8 +2716,7 @@ void FixBondReact::update_everything() // used when deleting atoms int ndel,ndelone; - int *mark = new int[nlocal]; - for (int i = 0; i < nlocal; i++) mark[i] = 0; + int *mark; tagint *tag = atom->tag; AtomVec *avec = atom->avec; @@ -2694,6 +2755,20 @@ void FixBondReact::update_everything() rxnID = global_mega_glove[0][i]; // reactions already shuffled from dedup procedure, so can skip first N if (iskip[rxnID]++ < nghostlyskips[rxnID]) continue; + + // we can insert atoms here, now that reactions are finalized + // can't do it any earlier, due to skipped reactions (max_rxn) + // reactions that create atoms are always treated as 'global' + if (create_atoms_flag[rxnID] == 1) { + onemol = atom->molecules[unreacted_mol[rxnID]]; + twomol = atom->molecules[reacted_mol[rxnID]]; + if (insert_atoms(global_mega_glove,i)) + ; else { // create aborted + reaction_count_total[rxnID]--; + continue; + } + } + for (int j = 0; j < max_natoms+1; j++) update_mega_glove[j][update_num_mega] = global_mega_glove[j][i]; update_num_mega++; @@ -2702,6 +2777,9 @@ void FixBondReact::update_everything() delete [] iskip; // mark to-delete atoms + nlocal = atom->nlocal; + mark = new int[nlocal]; + for (int i = 0; i < nlocal; i++) mark[i] = 0; for (int i = 0; i < update_num_mega; i++) { rxnID = update_mega_glove[0][i]; onemol = atom->molecules[unreacted_mol[rxnID]]; @@ -2731,25 +2809,6 @@ void FixBondReact::update_everything() } } - //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,"Bond/react: 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++) { @@ -3205,6 +3264,279 @@ void FixBondReact::update_everything() } } +/* ---------------------------------------------------------------------- +insert created atoms +------------------------------------------------------------------------- */ + +int FixBondReact::insert_atoms(tagint **my_mega_glove, int iupdate) +{ + // inserting atoms based off fix_deposit->pre_exchange + int flag; + imageint *imageflags; + double **coords,lamda[3],rotmat[3][3],vnew[3]; + double *newcoord; + double **v = atom->v; + double t,delx,dely,delz,rsq; + + memory->create(coords,twomol->natoms,3,"bond/react:coords"); + memory->create(imageflags,twomol->natoms,"bond/react:imageflags"); + + double *sublo,*subhi; + if (domain->triclinic == 0) { + sublo = domain->sublo; + subhi = domain->subhi; + } else { + sublo = domain->sublo_lamda; + subhi = domain->subhi_lamda; + } + + // find current max atom and molecule IDs + tagint *tag = atom->tag; + double **x = atom->x; + tagint *molecule = atom->molecule; + int nlocal = atom->nlocal; + + tagint maxtag_all,maxmol_all; + tagint max = 0; + for (int i = 0; i < nlocal; i++) max = MAX(max,tag[i]); + MPI_Allreduce(&max,&maxtag_all,1,MPI_LMP_TAGINT,MPI_MAX,world); + + max = 0; + for (int i = 0; i < nlocal; i++) max = MAX(max,molecule[i]); + MPI_Allreduce(&max,&maxmol_all,1,MPI_LMP_TAGINT,MPI_MAX,world); + + int dimension = domain->dimension; + + // only proc that owns reacting atom (use ibonding), + // fits post-reaction template to reaction site, for creating atoms + int n2superpose = 0; + for (int j = 0; j < twomol->natoms; j++) { + if (modify_create_fragid[rxnID] >= 0) + if (!twomol->fragmentmask[modify_create_fragid[rxnID]][j]) continue; + if (!create_atoms[j][rxnID] && !delete_atoms[equivalences[j][1][rxnID]][rxnID]) + n2superpose++; + } + + int ifit = atom->map(my_mega_glove[ibonding[rxnID]+1][iupdate]); // use this local ID to find fitting proc + Superpose3D superposer(n2superpose); + int fitroot = 0; + if (ifit >= 0 && ifit < atom->nlocal) { + fitroot = me; + + // get 'temperatere' averaged over site, used for created atoms' vels + t = get_temperature(my_mega_glove,1,iupdate); + + double **xfrozen; // coordinates for the "frozen" target molecule + double **xmobile; // coordinates for the "mobile" molecule + memory->create(xfrozen,n2superpose,3,"bond/react:xfrozen"); + memory->create(xmobile,n2superpose,3,"bond/react:xmobile"); + tagint iatom; + tagint iref = -1; // choose first atom as reference + int fit_incr = 0; + for (int j = 0; j < twomol->natoms; j++) { + if (modify_create_fragid[rxnID] >= 0) + if (!twomol->fragmentmask[modify_create_fragid[rxnID]][j]) continue; + int ipre = equivalences[j][1][rxnID]-1; // equiv pre-reaction template index + if (!create_atoms[j][rxnID] && !delete_atoms[ipre][rxnID]) { + if (atom->map(my_mega_glove[ipre+1][iupdate]) < 0) { + printf("WARNING: eligible atoms skipped for created-atoms fit on %d\n",me); + continue; + } + iatom = atom->map(my_mega_glove[ipre+1][iupdate]); + if (iref == -1) iref = iatom; + iatom = domain->closest_image(iref,iatom); + for (int k = 0; k < 3; k++) { + xfrozen[fit_incr][k] = x[iatom][k]; + xmobile[fit_incr][k] = twomol->x[j][k]; + } + fit_incr++; + } + } + double rmsd = superposer.Superpose(xfrozen, xmobile); + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + rotmat[i][j] = superposer.R[i][j]; + memory->destroy(xfrozen); + memory->destroy(xmobile); + } + MPI_Allreduce(MPI_IN_PLACE,&fitroot,1,MPI_INT,MPI_SUM,world); + MPI_Bcast(&t,1,MPI_DOUBLE,fitroot,world); + + // get coordinates and image flags + for (int m = 0; m < twomol->natoms; m++) { + if (create_atoms[m][rxnID] == 1) { + // apply optimal rotation/translation for created atom coords + // also map coords back into simulation box + if (fitroot == me) { + MathExtra::matvec(rotmat,twomol->x[m],coords[m]); + for (int i = 0; i < 3; i++) coords[m][i] += superposer.T[i]; + imageflags[m] = atom->image[ifit]; + domain->remap(coords[m],imageflags[m]); + } + MPI_Bcast(&imageflags[m],1,MPI_LMP_IMAGEINT,fitroot,world); + MPI_Bcast(coords[m],3,MPI_DOUBLE,fitroot,world); + } + } + + // check distance between any existing atom and inserted atom + // if less than near, abort + if (overlapsq[rxnID] > 0) { + int abortflag = 0; + for (int m = 0; m < twomol->natoms; m++) { + if (create_atoms[m][rxnID] == 1) { + for (int i = 0; i < nlocal; i++) { + delx = coords[m][0] - x[i][0]; + dely = coords[m][1] - x[i][1]; + delz = coords[m][2] - x[i][2]; + domain->minimum_image(delx,dely,delz); + rsq = delx*delx + dely*dely + delz*delz; + if (rsq < overlapsq[rxnID]) { + abortflag = 1; + break; + } + } + if (abortflag) break; + } + } + MPI_Allreduce(MPI_IN_PLACE,&abortflag,1,MPI_INT,MPI_MAX,world); + if (abortflag) { + memory->destroy(coords); + memory->destroy(imageflags); + return 0; + } + } + + // clear ghost count and any ghost bonus data internal to AtomVec + // same logic as beginning of Comm::exchange() + // do it now b/c inserting atoms will overwrite ghost atoms + atom->nghost = 0; + atom->avec->clear_bonus(); + + // check if new atoms are in my sub-box or above it if I am highest proc + // if so, add atom to my list via create_atom() + // initialize additional info about the atoms + // set group mask to "all" plus fix group + int preID; // new equivalences index + int add_count = 0; + for (int m = 0; m < twomol->natoms; m++) { + if (create_atoms[m][rxnID] == 1) { + // increase atom count + add_count++; + preID = onemol->natoms+add_count; + + if (domain->triclinic) { + domain->x2lamda(coords[m],lamda); + newcoord = lamda; + } else newcoord = coords[m]; + + flag = 0; + if (newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] && + newcoord[1] >= sublo[1] && newcoord[1] < subhi[1] && + newcoord[2] >= sublo[2] && newcoord[2] < subhi[2]) flag = 1; + else if (dimension == 3 && newcoord[2] >= domain->boxhi[2]) { + if (comm->layout != Comm::LAYOUT_TILED) { + if (comm->myloc[2] == comm->procgrid[2]-1 && + newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] && + newcoord[1] >= sublo[1] && newcoord[1] < subhi[1]) flag = 1; + } else { + if (comm->mysplit[2][1] == 1.0 && + newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] && + newcoord[1] >= sublo[1] && newcoord[1] < subhi[1]) flag = 1; + } + } else if (dimension == 2 && newcoord[1] >= domain->boxhi[1]) { + if (comm->layout != Comm::LAYOUT_TILED) { + if (comm->myloc[1] == comm->procgrid[1]-1 && + newcoord[0] >= sublo[0] && newcoord[0] < subhi[0]) flag = 1; + } else { + if (comm->mysplit[1][1] == 1.0 && + newcoord[0] >= sublo[0] && newcoord[0] < subhi[0]) flag = 1; + } + } + + int root = 0; + if (flag) { + root = me; + + atom->avec->create_atom(twomol->type[m],coords[m]); + int n = atom->nlocal - 1; + atom->tag[n] = maxtag_all + add_count; + + // locally update mega_glove + my_mega_glove[preID][iupdate] = atom->tag[n]; + + if (atom->molecule_flag) { + if (twomol->moleculeflag) { + atom->molecule[n] = maxmol_all + twomol->molecule[m]; + } else { + atom->molecule[n] = maxmol_all + 1; + } + } + + atom->mask[n] = 1 | groupbit; + atom->image[n] = imageflags[m]; + + // guess a somewhat reasonable initial velocity based on reaction site + // further control is possible using bond_react_MASTER_group + // compute |velocity| corresponding to a given temperature t, using specific atom's mass + double vtnorm = sqrt(t / (force->mvv2e / (dimension * force->boltz)) / atom->mass[twomol->type[m]]); + v[n][0] = random[rxnID]->uniform(); + v[n][1] = random[rxnID]->uniform(); + v[n][2] = random[rxnID]->uniform(); + double vnorm = sqrt(v[n][0]*v[n][0] + v[n][1]*v[n][1] + v[n][2]*v[n][2]); + v[n][0] = v[n][0]/vnorm*vtnorm; + v[n][1] = v[n][1]/vnorm*vtnorm; + v[n][2] = v[n][2]/vnorm*vtnorm; + modify->create_attribute(n); + + // initialize group statuses + // why aren't these more global... + int flag; + int index1 = atom->find_custom("limit_tags",flag); + int *i_limit_tags = atom->ivector[index1]; + + int *i_statted_tags; + if (stabilization_flag == 1) { + int index2 = atom->find_custom(statted_id,flag); + i_statted_tags = atom->ivector[index2]; + } + + int index3 = atom->find_custom("react_tags",flag); + int *i_react_tags = atom->ivector[index3]; + + i_limit_tags[n] = update->ntimestep + 1; + if (stabilization_flag == 1) i_statted_tags[n] = 0; + i_react_tags[n] = rxnID; + } + // globally update mega_glove and equivalences + MPI_Allreduce(MPI_IN_PLACE,&root,1,MPI_INT,MPI_SUM,world); + MPI_Bcast(&my_mega_glove[preID][iupdate],1,MPI_LMP_TAGINT,root,world); + equivalences[m][0][rxnID] = m+1; + equivalences[m][1][rxnID] = preID; + reverse_equiv[preID-1][0][rxnID] = preID; + reverse_equiv[preID-1][1][rxnID] = m+1; + } + } + + // reset global natoms + // if global map exists, reset it now instead of waiting for comm + // since other pre-exchange fixes may use it + // invoke map_init() b/c atom count has grown + atom->natoms += add_count; + if (atom->natoms < 0) + error->all(FLERR,"Too many total atoms"); + maxtag_all += add_count; + if (maxtag_all >= MAXTAGINT) + error->all(FLERR,"New atom IDs exceed maximum allowed ID"); + if (atom->map_style != Atom::MAP_NONE) { + atom->map_init(); + atom->map_set(); + } + // atom creation successful + memory->destroy(coords); + memory->destroy(imageflags); + return 1; +} + /* ---------------------------------------------------------------------- read superimpose file ------------------------------------------------------------------------- */ @@ -3222,6 +3554,7 @@ void FixBondReact::read(int myrxn) // skip blank lines or lines that start with "#" // stop when read an unrecognized line + ncreate = 0; while (1) { readline(line); @@ -3240,6 +3573,7 @@ void FixBondReact::read(int myrxn) "equal number of atoms in reaction templates"); } else if (strstr(line,"deleteIDs")) sscanf(line,"%d",&ndelete); + else if (strstr(line,"createIDs")) sscanf(line,"%d",&ncreate); else if (strstr(line,"chiralIDs")) sscanf(line,"%d",&nchiral); else if (strstr(line,"constraints")) { sscanf(line,"%d",&nconstraints[myrxn]); @@ -3276,6 +3610,8 @@ void FixBondReact::read(int myrxn) Equivalences(line, myrxn); } else if (strcmp(keyword,"DeleteIDs") == 0) { DeleteAtoms(line, myrxn); + } else if (strcmp(keyword,"CreateIDs") == 0) { + CreateAtoms(line, myrxn); } else if (strcmp(keyword,"ChiralIDs") == 0) { ChiralCenters(line, myrxn); } else if (strcmp(keyword,"Constraints") == 0) { @@ -3312,7 +3648,7 @@ void FixBondReact::Equivalences(char *line, int myrxn) for (int i = 0; i < nequivalent; i++) { readline(line); sscanf(line,"%d %d",&tmp1,&tmp2); - if (tmp1 > onemol->natoms || tmp2 > onemol->natoms) + if (tmp1 > onemol->natoms || tmp2 > twomol->natoms) error->one(FLERR,"Bond/react: Invalid template atom ID in map file"); //equivalences is-> clmn 1: post-reacted, clmn 2: pre-reacted equivalences[tmp2-1][0][myrxn] = tmp2; @@ -3335,6 +3671,17 @@ void FixBondReact::DeleteAtoms(char *line, int myrxn) } } +void FixBondReact::CreateAtoms(char *line, int myrxn) +{ + create_atoms_flag[myrxn] = 1; + int tmp; + for (int i = 0; i < ncreate; i++) { + readline(line); + sscanf(line,"%d",&tmp); + create_atoms[tmp-1][myrxn] = 1; + } +} + void FixBondReact::CustomCharges(int ifragment, int myrxn) { for (int i = 0; i < onemol->natoms; i++) diff --git a/src/USER-REACTION/fix_bond_react.h b/src/USER-REACTION/fix_bond_react.h index b5922de327..87a5945d45 100644 --- a/src/USER-REACTION/fix_bond_react.h +++ b/src/USER-REACTION/fix_bond_react.h @@ -67,6 +67,9 @@ class FixBondReact : public Fix { int custom_exclude_flag; int *stabilize_steps_flag; int *custom_charges_fragid; + int *create_atoms_flag; + int *modify_create_fragid; + double *overlapsq; int *molecule_keyword; int maxnconstraints; int *nconstraints; @@ -84,10 +87,10 @@ class FixBondReact : public Fix { int max_natoms; // max natoms in a molecule template tagint *partner,*finalpartner; double **distsq,*probability; - int *ncreate; - int maxcreate; - int allncreate; - tagint ***created; + int *nattempt; + int maxattempt; + int allnattempt; + tagint ***attempt; class Molecule *onemol; // pre-reacted molecule template class Molecule *twomol; // post-reacted molecule template @@ -117,7 +120,7 @@ class FixBondReact : public Fix { int *ibonding,*jbonding; int *closeneigh; // indicates if bonding atoms of a rxn are 1-2, 1-3, or 1-4 neighbors - int nedge,nequivalent,ndelete,nchiral; // # edge, equivalent atoms in mapping file + int nedge,nequivalent,ndelete,ncreate,nchiral; // # edge, equivalent atoms in mapping file int attempted_rxn; // there was an attempt! int *local_rxn_count; int *ghostly_rxn_count; @@ -133,6 +136,7 @@ class FixBondReact : public Fix { int **landlocked_atoms; // all atoms at least three bonds away from edge atoms int **custom_charges; // atoms whose charge should be updated int **delete_atoms; // atoms in pre-reacted templates to delete + int **create_atoms; // atoms in post-reacted templates to create int ***chiral_atoms; // pre-react chiral atoms. 1) flag 2) orientation 3-4) ordered atom types int **nxspecial,**onemol_nxspecial,**twomol_nxspecial; // full number of 1-4 neighbors @@ -156,6 +160,7 @@ class FixBondReact : public Fix { void EdgeIDs(char *, int); void Equivalences(char *, int); void DeleteAtoms(char *, int); + void CreateAtoms(char *,int); void CustomCharges(int, int); void ChiralCenters(char *, int); void ReadConstraints(char *, int); @@ -169,7 +174,7 @@ class FixBondReact : public Fix { void ring_check(); int check_constraints(); void get_IDcoords(int, int, double *); - double get_temperature(); + double get_temperature(tagint **, int, int); int get_chirality(double[12]); // get handedness given an ordered set of coordinates void open(char *); @@ -185,6 +190,7 @@ class FixBondReact : public Fix { void glove_ghostcheck(); void ghost_glovecast(); void update_everything(); + int insert_atoms(tagint **, int); void unlimit_bond(); void limit_bond(int); void dedup_mega_gloves(int); //dedup global mega_glove