Degeneracy & Multiple TS

[goldmanm]

This PR solves two related problems in RMG: reaction degeneracy values that take into account resonance structures and allows multiple transition states to be stored in RMG's output file. This PR will solve underterminable kinetics problems (#494) , muptiple TS (#323, #142), and proper degeneracy values (#876, #833, #619, #552) . These two problems both involve updates to the degeneracy-finding algorithm, which is why they are in one PR. The major modification for this pull request is the order in which the generate reactions steps take place. The reactions of all resonance structures are now generated before algorithms like degeneracy or finding reverse kinetics are run. In addition, the degeneracy algorithm was redone to take into account which atoms have actually been switched (this correct the excessive degneracy from certain classes # ) and the templates that are matched (to allow multiple transition states). In addition, the adding of new reactions to CoreEdgeReactionModel also checks the family and template, and marks duplicate reactions when they exist (allowing for multiple TS and removing underterminatee kinetics problems)

1. redefines the degeneracy variable as a private float. When degeneracy is updated, the kinetic rate (if set) is too.
2. modularizes parts of generateReactions - finding rxn.reverse and degeneracy are separate methods
3. runs reactions betwween all resonance structures of species before finding degeneracy
4. utilizes reaction templates to find and keep multiple transition states (as separate reactions)
5. uses atomIDs to identify when a reaction path has already been added to degeneracy
6. adding proper cleanup to many test methods (including solvation database cleanup)
7. some code cleanup and other modularization of methods.

This PR is ready for code review. Testing this could involve merging it and running a job of interest, though I think RMG-Tests is already doing that. If it is merged, PR #323 can be removed.

[mjohnson541]

Before this branch it was virtually impossible to pick up the propane low T ketohydroperoxide feedback loop without a library due to a multiple TS issue with a key chain branching reaction. I can confirm that this pull request enables RMG to now pick up this key pathway at appropriate temperatures and feasible tolerances.

[mjohnson541]

Why do we want to be able to generate entirely random atom ID's rather than just increment them from a set value?

[rwest]

...and then if they're not unique (not self.atomIDValid()) we make them sequentially increasing from a random start? (Suggesting sequential is fine and it needn't be random?). Just using sequential sounds simpler to me, unless Matt and I are both missing something.

[goldmanm]

Good question. Generating IDs from a starting value (like previously going 1 to N) is almost certainly going to cause clashes of atom IDs after two labeled atoms are merged in generateReactions. Generating IDs over the entire range of integer values is much more likely to reduce the clashing IDs when merging two molecules (which is what random IDs seeks to accomplish). When a clash occurs, we would need to reassign IDs properly for all resonance structures (possibly remaking resonance structures which is quite time intensive)

One alternative to completely random numbers could be to label starting from -2^-15 till 2^15 and then loop back over after labeling ~65,000 atoms (not resetting after each call). If you know that this would reduce conflicts and/or computational time, let me know.

The atomIDs cannot be set post-merging since this can result in incorrect degeneracy values when trying to generate reactions when multiple resonance structures are present. Lets say you have two reactants reacting together, and each has multiple resonance structures. In generateReactions, all combinations of resonance structures will be merged and the possible reactions output. Then findDegeneracies would use the atomIDs for all the sets of products to find the degeneracy. If different resonance structures are labeled differently (which can happen when relabeling after the merging of structures), the degeneracy could be higher than expected.

[goldmanm]

I can make a class variable that tracks the entire range of integers to avoid overlap, and then all the atom IDs can be added sequentially. That seems more intuitive than random numbers and probably will save computational time by not checking atomIDValid when setting the IDs. Do you think that would be preferable? Any other ideas on how to avoid conflicting IDs when merging without random Numbers?

[rwest]

Now that I understand what the IDs are doing, I like @goldmanm's last suggestion of a class variable that gradually increments.

[mliu49]

I only skimmed through the changes and made a few minor comments. I'm very excited to try it out.

Have you tested to see if pdep is affected at all? I ran into some issues with it when I tried changing generateReactions and deflate/inflate to work with species objects.

[mliu49]

Should this line be indented less one level?

[goldmanm]

put an else statement there in squashed commit.

[mliu49]

I think it's generally undesirable to have mutable default arguments. Was there a reason for this change?

[goldmanm]

This was an old commit from awhile ago (which I fixed), I am not sure how it got in here. I will remove this commit from the PR.

[mliu49]

Upper limit should be 2**15?

[goldmanm]

I put it lower to ensure it could still count up, but this will be replaced with a more intuitive method.

[mliu49]

Unfinished comment?

[goldmanm]

fixed comment in squashed commit.

[mliu49]

Should this say Species or Molecule Objects?

[goldmanm]

Yup! Fixed in squashed commit.

[mliu49]

I think you can check identical only if isomorphic is true or check isomorphic only if identical is false. Either way should reduce the number of calls made.

• If two reactions are identical, then they are isomorphic.
• If two reactions are not isomorphic, then they are not identical.

[goldmanm]

Great idea! This would be an awesome speedup. Added this in a squash commit (for the second case under the assumption that not isomorphic is more common than is identical)

[mjohnson541]

Wouldn't it be much faster and shorter to iterate through the lists and pop out entries that are the same from both lists as you go. I believe the below code would reduce the maximum calls to comparison_method to 6 rather than the current 18.

list1Temp = []
list2Temp = []
indexList = []
k = 0
while list1 != []:
   entry = list1[k]
   for i,entry2 in enumerate(list2):
       if comparison_method(entry,entry2):
           temp = list1.pop(k)
           list1temp.append(temp)
           temp = list2.pop(i)
           list2temp.append(temp)
           indexList.append(i)
           break
   else: #loop doesn't break, something is missing, templist2 is incomplete
         for j in xrange(0,len(list2)):
             if not j in indexList:
                indexList.append(j)
         list1 = list1Temp+list1
         list2 = list2Temp+list2
         list2 = [list2[i] for i in indexList]
         return False
   k += 1
#loop didn't break, list1 is fine, resort list2 so it is unchanged
list2 = [list2Temp[i] for i in indexList]
list1 = list1Temp
return True

[mjohnson541]

Or wait you don't even need to avoid copying the entries as long as you don't deepcopy them do you?

list1Temp = list1[:]
list2Temp = list2[:]
for entry in list1Temp:
   for entry2 in list2Temp:
       if comparison_method(entry,entry2):
           list2Temp.remove(entry2)
           break
   else: #loop doesn't break, something is missing
         return False
#all loops completed
return True

[rwest]

Something must be wrong with this last example. Although I'm not 100% sure what it's trying to achieve, will the while loop ever end? You never remove anything fromlist1Temp. In any case, the alternative is not obviously faster, with lots of list manipulation, creation, popping, etc. I have some thoughts on probable efficiency gains, but before sharing them can anyone tell me if this matters? Does it get called a lot? (sorry for not being aware of the context - the discussion caught my eye).

[mjohnson541]

Oh, yes, sorry lost that in translation, I think you could actually just do "for entry in list1Temp" since you don't need k anymore. I thought this would be faster under the assumption that most of our time is spent in the isomorphism checks this avoids (maximum for 2-2 lists drops from 4 to 3, maximum for 3-3 lists drops from 18 to 6). Is that not the case?

[goldmanm]

@rwest, this method is called every time Reaction isomorphism is checked, which is done in degeneracy calculations, adding reactions to the edge, etc. I think it is relatively commonly called (though this is somewhat subjective without a profiling comparison). This method modularized what already existed in reaction.isIsomorphic, though it would be awesome to get it working faster too.

@KEHANG, Is there any listing of previous profiling to see how important the Reaction.isIsomorphic method is?

[goldmanm]

this notebook profiles the recommendations to isisomorphic by @mjohnson541 for H_abstraction reactions. It seems to be faster most of the time by about 10%, though sometimes it is slower (which could be due to other processes demanding resources). There should likely be more benefit if we are looking at 3x3 reactions. Feel free to play around with the notebook, and let me know if you'd like the suggested algorithm implemented here.

[mjohnson541]

You could also avoid defining the function by passing the function pointer instead

if checkIdentical:
    comparison_method = molecule.isIdentical
else:
    comparison_method = molecule.isIsomorphic

I'm really not sure how significant the time savings would be, but the above timings suggest we do have some noticeable overhead.

[mliu49]

I like the new approach, but I don't think the time savings would be as significant as you would expect. See the table below for a comparison of the number of function calls. Note that list.remove performs an isomorphism comparison with each item in the list until it finds the one it needs to remove. Using enumerate along with del instead of list.remove might help a little bit.

reactants	case	if/elif method	loop method
1		1 comparison call	1 comparison call 1 list.remove call
2	best	2 comparison calls	2 comparison calls 2 list.remove calls
2	worst	4 comparison calls	3 comparison calls 2 list.remove calls
3	best	3 comparison calls	3 comparison calls 3 list.remove calls
3	worst	18 comparison calls	6 comparison calls 3 list.remove calls

[goldmanm]

implemented the updated deletion algorithm, in the ipynb and still no improvement in efficiency :(

[mjohnson541]

So if I came across this comment randomly in the code without reading the function I would have thought the above comments implied that multiple TS's were simply treated as degeneracies for purposes of this function especially considering the docstring doesn't mention multiple TS's being added as duplicate reactions. Can you make that a bit more clear?

[mliu49]

This was my comment from earlier changes, before we decided to tackle multiple transition states. In this comment, different transition states refers to ones that are the effectively identical but involving a different set of atoms, which i guess you could also call degenerate transition states.

It should probably be changed now that we're also dealing with actually different transition states.

[goldmanm]

updated in the most recent commit.

[mjohnson541]

Is there a reason we can't add some boolean property to the reaction object during rate estimation (such as "fromReverse") so we already know if it's the reverse template or not and we don't have to do all of these isomorphism checks?