-
Notifications
You must be signed in to change notification settings - Fork 0
/
reaction.py
177 lines (157 loc) · 6.12 KB
/
reaction.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
import rbnmol
import AChemKit
import AChemKit.utils
import AChemKit.utils.bag
from AChemKit.utils.bag import OrderedFrozenBag
def decompose(*mols):
assert len(mols) > 0
if len(mols) == 1:
newmols = mols[0].decomposition()
else:
newmols = []
for mol in mols:
for outmol in decompose(mol):
newmols.append(outmol)
assert len(newmols) > 0
newmols = tuple(newmols)
if newmols != mols:
return decompose(*newmols)
else:
return mols
def validate(mol):
if mol.composing is None:
#top level, can drill through it checking bonding sites
assert 1 not in mol.rbn.bonding.values(), product
drill = mol
while drill.composition is not None:
drill = drill.composition[0]
assert drill.rbn.bonding[0] == 0, drill
drill = mol
while drill.composition is not None:
drill = drill.composition[-1]
assert drill.rbn.bonding[max(drill.rbn.bonding)] == 0, drill
else:
assert mol.composing.composition is not None
assert mol in mol.composing.composition
if mol.composition is None:
return
else:
assert len(mol.composition) > 1
for test in mol.composition:
validate(test)
def react(reactants):
#dont do complicated flipping-ness
#reactions are either A+B or B+A. Assume this has
#been done before this function.
#This determines the outcome of one reaction between A and B
#An alternative would be to determine all outcomes
#and then pick one.
if isinstance(reactants[0], list) or isinstance(reactants[0], tuple):
reactants = reactants[0]
A00, B00 = reactants
#go down to the smallest components in each one
while A00.composition is not None:
A00 = A00.composition[-1]
while B00.composition is not None:
B00 = B00.composition[0]
reacting = False
while reacting is False:
#if these can react, react
if A00.bonding_criterion(B00):
reacting = True
else:
if A00.composing is None and B00.composing is None:
#nothing has reacted
#nothing is composing larger things
#end reaction
return OrderedFrozenBag(reactants)
elif A00.composing is not None and B00.composing is None:
#only one has a larger structure
#go to the larger structure
A00 = A00.composing
elif A00.composing is None and B00.composing is not None:
B00 = B00.composing
#only one has a larger structure
#go to the larger structure
elif A00.composing is not None and B00.composing is not None:
#both of them have a larger structure
#which one do we use?
#go for the smaller one
if A00.composing.rbn.n < B00.composing.rbn.n:
#As composing structure is smaller
#use it
A00 = A00.composing
elif A00.composing.rbn.n > B00.composing.rbn.n:
#Bs composing structure is smaller
#use it
B00 = B00.composing
elif A00.composing.rbn.n == B00.composing.rbn.n:
#both the composing stuctures are the same size
#use both of them
A00 = A00.composing
B00 = B00.composing
assert reacting is True
#fill in bonding sites in only these bRBNs
A01 = A00.fill_this_bonding("r")
B10 = B00.fill_this_bonding("l")
A01_top, upstepsA = A01.top_steps()
B10_top, upstepsB = B10.top_steps()
products = ((A01_top,), (B10_top,))
#check nothing breaks
products = (decompose(A01_top), decompose(B10_top))
A01_top = products[0][-1]
B10_top = products[-1][0]
A01 = A01_top.drill("r", upstepsA)
B10 = B10_top.drill("l", upstepsB)
#see if bonding can continue
if A01.bonding_criterion(B10):
#if this is the top-most, then create a new composite and return it
#have to do complicated extensions stuff if it is not
#but this function will deal with that...
A01 = A01.fill_all_bonding("r")
B10 = B10.fill_all_bonding("l")
C00 = A01.extend(B10)
C00_top, upsteps = C00.top_steps()
products = products[0][:-1]+(C00_top,)+products[1][1:]
else:
#need to revert the bonding process
A00_2 = A01.empty_all_bonding("r")
B00_2 = B10.empty_all_bonding("l")
#get the top so it is easier to make products
#upsteps is a consequence of this, but not needed here
A00_top, upsteps = A00_2.top_steps()
B00_top, upsteps = B00_2.top_steps()
products = products[0][:-1]+(A00_top, B00_top)+products[1][1:]
#check nothing breaks again
products = decompose(*products)
#collapse single compositions
newproducts = []
for mol in products:
newproducts.append(mol.collapse())
products = newproducts
if __debug__:
for product in products:
validate(product)
return OrderedFrozenBag(products)
def all_reactions(*reactants):
a,b = reactants
rates = {}
for ordering in ((a,b), (b,a)):
reactants = OrderedFrozenBag(ordering)
products = react(ordering)
#sanity check the total atoms are the same
if __debug__:
reactantsatoms = ""
for mol in reactants:
reactantsatoms += mol.atomstr()
productsatoms = ""
for mol in products:
productsatoms += mol.atomstr()
assert reactantsatoms == productsatoms, [reactantsatoms, productsatoms, tuple(reactants), tuple(products)]
for mol in reactants:
assert mol.composing is None, mol
for mol in products:
assert mol.composing is None, mol
if products != reactants:
rates[reactants, products] = 1.0
return rates