Save labeled adjacency lists of all reactions

[rwest]

Motivation or Problem

For third-party tools wishing to calculate reactions, eg. via automated TST calculations, it's helpful to have adjacency lists of reactant and product structures, with the atom order invariant, and with reaction recipe labels on the atoms.

Description of Changes

With this PR, if you put generateLabeledReactions=True in the options section of the input file, then it will save a YAML file with the adjacency list for every reaction. If you also have saveEdgeSpecies=True then it'll save edge reactions too.

Testing

Added the generateLabeledReactions=True option to the examples/rmg/catalysis/ch4_o2/input.py example input file and ran it.

Work in Progress

• [ ] Add unit tests?
• [ ] Add documentation

please add to this list as appropriate

[mazeau]

when trying to run the meoh synthesis model, I get this error:

/home/mazeau.e/.conda/envs/rmg_env/lib/python3.7/site-packages/scipy/optimize/optimize.py:1978: LinAlgWarning: Ill-conditioned matrix (rcond=9.79183e-18): result may not be accurate.
  fu = func(x, *args)
/home/mazeau.e/.conda/envs/rmg_env/lib/python3.7/site-packages/scipy/optimize/optimize.py:1978: LinAlgWarning: Ill-conditioned matrix (rcond=5.9437e-17): result may not be accurate.
  fu = func(x, *args)
Traceback (most recent call last):
  File "/home/mazeau.e/Code/RMG-Py/rmg.py", line 118, in <module>
    main()
  File "/home/mazeau.e/Code/RMG-Py/rmg.py", line 100, in main
    rmg_profiler.runctx(command, global_vars, local_vars)
  File "/home/mazeau.e/.conda/envs/rmg_env/lib/python3.7/cProfile.py", line 100, in runctx
    exec(cmd, globals, locals)
  File "<string>", line 1, in <module>
  File "/home/mazeau.e/Code/RMG-Py/rmgpy/rmg/main.py", line 737, in execute
    trimolecular_react=self.trimolecular_react)
  File "/home/mazeau.e/Code/RMG-Py/rmgpy/rmg/model.py", line 595, in enlarge
    procnum=procnum)
  File "/home/mazeau.e/Code/RMG-Py/rmgpy/rmg/react.py", line 172, in react_all
    return react(spc_fam_tuples, procnum), [fam_tuple[0] for fam_tuple in spc_fam_tuples]
  File "/home/mazeau.e/Code/RMG-Py/rmgpy/rmg/react.py", line 66, in react
    reactions = list(map(_react_species_star, spc_fam_tuples))
  File "/home/mazeau.e/Code/RMG-Py/rmgpy/rmg/react.py", line 79, in _react_species_star
    return react_species(*args)
  File "/home/mazeau.e/Code/RMG-Py/rmgpy/rmg/react.py", line 94, in react_species
    reactions = get_db('kinetics').generate_reactions_from_families(species_tuple, only_families=only_families)
  File "/home/mazeau.e/Code/RMG-Py/rmgpy/data/kinetics/database.py", line 534, in generate_reactions_from_families
    prod_resonance=resonance))
  File "/home/mazeau.e/Code/RMG-Py/rmgpy/data/kinetics/database.py", line 561, in react_molecules
    prod_resonance=prod_resonance))
  File "/home/mazeau.e/Code/RMG-Py/rmgpy/data/kinetics/family.py", line 1727, in generate_reactions
    self._generate_reactions(reactants, products=products, forward=True, prod_resonance=prod_resonance))
  File "/home/mazeau.e/Code/RMG-Py/rmgpy/data/kinetics/family.py", line 2122, in _generate_reactions
    rxn = self._create_reaction(reactant_structures, product_structures, forward)
  File "/home/mazeau.e/Code/RMG-Py/rmgpy/data/kinetics/family.py", line 1648, in _create_reaction
    if same_species_lists(reactants, products):
TypeError: Argument 'list2' has incorrect type (expected list, got rmgpy.molecule.molecule.Molecule)

[mazeau]

I think the change to line 1626 is the cause for this error because same_species_lists expects lists and can't deal with the product structures

[mazeau]

After the rebase, I still see the same error:

Initialization complete. Starting model generation.

conditions choosen for reactor 0 were: T = 483.0 K, 
Generating initial reactions...
For reaction generation 1 process is used.
Error: Problem family: Surface_Adsorption_Dissociative
Error: Problem reactants: (Molecule(smiles="[Pt]"), Molecule(smiles="[H][H]"))
Error: <Molecule "[Pt]">
1 *3 X u0 p0 c0

Error: <Molecule "[H][H]">
1 *1 H u0 p0 c0 {2,S}
2 *2 H u0 p0 c0 {1,S}

Traceback (most recent call last):
  File "/Users/emilymazeau/Code/RMG-Py/rmg.py", line 118, in <module>
    main()
  File "/Users/emilymazeau/Code/RMG-Py/rmg.py", line 112, in main
    rmg.execute(**kwargs)
  File "/Users/emilymazeau/Code/RMG-Py/rmgpy/rmg/main.py", line 737, in execute
    trimolecular_react=self.trimolecular_react)
  File "/Users/emilymazeau/Code/RMG-Py/rmgpy/rmg/model.py", line 599, in enlarge
    procnum=procnum)
  File "/Users/emilymazeau/Code/RMG-Py/rmgpy/rmg/react.py", line 172, in react_all
    return react(spc_fam_tuples, procnum), [fam_tuple[0] for fam_tuple in spc_fam_tuples]
  File "/Users/emilymazeau/Code/RMG-Py/rmgpy/rmg/react.py", line 66, in react
    reactions = list(map(_react_species_star, spc_fam_tuples))
  File "/Users/emilymazeau/Code/RMG-Py/rmgpy/rmg/react.py", line 79, in _react_species_star
    return react_species(*args)
  File "/Users/emilymazeau/Code/RMG-Py/rmgpy/rmg/react.py", line 94, in react_species
    reactions = get_db('kinetics').generate_reactions_from_families(species_tuple, only_families=only_families)
  File "/Users/emilymazeau/Code/RMG-Py/rmgpy/data/kinetics/database.py", line 534, in generate_reactions_from_families
    prod_resonance=resonance))
  File "/Users/emilymazeau/Code/RMG-Py/rmgpy/data/kinetics/database.py", line 561, in react_molecules
    prod_resonance=prod_resonance))
  File "/Users/emilymazeau/Code/RMG-Py/rmgpy/data/kinetics/family.py", line 1727, in generate_reactions
    self._generate_reactions(reactants, products=products, forward=True, prod_resonance=prod_resonance))
  File "/Users/emilymazeau/Code/RMG-Py/rmgpy/data/kinetics/family.py", line 2122, in _generate_reactions
    rxn = self._create_reaction(reactant_structures, product_structures, forward)
  File "/Users/emilymazeau/Code/RMG-Py/rmgpy/data/kinetics/family.py", line 1648, in _create_reaction
    if same_species_lists(reactants, products):
TypeError: Argument 'list2' has incorrect type (expected list, got rmgpy.molecule.molecule.Molecule)

[mazeau]

THis is the unitt test error that I get:

======================================================================
ERROR: ensure rxns with multiple transition states are kept as separate reactions
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/Users/emilymazeau/Code/RMG-Py/rmgpy/data/kinetics/kineticsTest.py", line 421, in test_separate_transition_states_generated_regardless_of_reactant_order
    reaction_list = family._generate_reactions([mol_a, mol_b], products=[mol_c, mol_d])
  File "/Users/emilymazeau/Code/RMG-Py/rmgpy/data/kinetics/family.py", line 1925, in _generate_reactions
    save_adjacency_lists = get_input('generate_labeled_reactions')
  File "/Users/emilymazeau/Code/RMG-Py/rmgpy/rmg/input.py", line 1240, in get_input
    raise Exception('Could not get variable with name: {}'.format(name))
Exception: Could not get variable with name: generate_labeled_reactions

[rwest]

So it now saves the files if you set

options(
    generateLabeledReactions=True,
)

and doesn't crash if you don't. The unit tests don't pass yet, but this will at least create files to send to Sandia

[rwest]

I added @yunsiechung as a reviewer, because some of the ideas in here may help her with her reaction atom-mapping question (for predicting ∆∆Gsolv of transition states). This PR probably needs another update.

[yunsiechung]

@rwest Thank you Prof. West for tagging me here! I'm happy to review this PR and it will help me understanding how atom can be mapped for reactions in RMG. But it looks like this branch has some conflicts. Is this active PR?

[rwest]

Looks like it was last rebased 10 months ago.

[github-actions[bot]]

This pull request is being automatically marked as stale because it has not received any interaction in the last 90 days. Please leave a comment if this is still a relevant pull request, otherwise it will automatically be closed in 30 days.

[rwest]

@mjohnson541 this is how we generated labeled adjacency lists of transition states in order to give to tools like pynta, back in the day. Is this still a helpful way to do it? Should this PR be brought back from the dead, rebased, and merged? Or is it totally redundant?

[mjohnson541]

Nominally merging such a PR would allow easy movement from an output RMG job to Pynta or similar tools. However, there are a couple important caveats. 1) Pynta is agnostic to RMG's families so it needs both the labeled reactants and products...I'm also running into some data analysis cases where its valuable to have the reactants and products be in the same order (implying a complete atom map) 2) Pynta is aware of "site reuse" in a way that RMG doesn't seem to be, for Pynta the distinction between a new empty site and a site that is occupied in the reactants is very important, but if I recall correctly as a result of this issue I was unable to use RMG to generate useful Pynta inputs from the reaction generator in earlier attempts...