Stoichiometry with substances on both sides

[rcolpo]

Hi,

I would like to calculate the stoichiometry of reactions like the one bellow:

reac, prod = balance_stoichiometry({'H4P2O7', 'HPO3', 'H2O'}, {'H4P2O7', 'HPO3'})

Answer would be:

H4P2O7 + 4 HPO3 + H2O <-> 2 H4P2O7 + 2 HPO3

OR

2 HPO3 + H2O = H4P2O7

Any idea on how this would be performed? At the moment I receive an error: "Stoichiometry with substances on both sides".

[bjodah]

It should be possible to devise an algorithm which comes up with the second solution (which looks preferable to me in the sense that it is arguably more "canoncial"). There is however no such algorithm implemented at the moment.

[rcolpo]

Thanks bjodah. Any idea on how this could be done?

[bjodah]

I would start considering small examples using "pen and paper" first. And solve a bunch of those to develop test cases and get an intuition.

C + CO + CO2 -> C + CO + CO2  # should fail
C + CO + CO2 -> C + CO        # suggested solution:  C +      CO2 ->     2CO
C + CO + CO2 -> C +      CO2  # suggested solution:      2CO      -> C +      CO2
C + CO + CO2 ->     CO + CO2  # suggested solution:  C +      CO2 ->     2CO
C + CO       -> C + CO + CO2  # suggested solution:      2CO      -> C +      CO2
C +      CO2 -> C + CO + CO2  # suggested solution:  C +      CO2 ->     2CO
    CO + CO2 -> C + CO + CO2  # suggested solution:      2CO      -> C +      CO2

here's a start. Do those solutions look like the canonical ones to you?

I think I can throw together a naïve implementation using brute force.

[rcolpo]

Hi Björn.

Yes, this is exactly the kind of solution I'm looking for. I'm not having much success in implementing a code to replicate these results. Thank you.

[bjodah]

@rcolpo I released chempy-0.7.6 with this feature. Make sure you pass underdetermined=None as well as allow_duplicates=True.

[rcolpo]

You are the best. Thank you so much. This will be very useful.

[bjodah]

Thanks! (please don't forget to cite the article mentioned in the README if you use it in an academic setting)

[rcolpo]

I will most certainly do it. Thank you so much.

[rcolpo]

Hi bjodah,

I receive an error (Failed to remove duplicate keys: ['Mn1']) with the following command:

balance_stoichiometry({'H2O2', 'Mn1', 'H1'}, {'Mn1', 'H2O1'}, allow_duplicates=True, underdetermined=None)

Solution should be: 2 H2O1 = H2O2 + 2 H1

Do you know what the problem is? Thank you.

[bjodah]

I see, here's one possible remedy: https://github.com/bjodah/chempy/pull/126

Please take a look at the tests and see if you agree on my choice of expected results.

[rcolpo]

Thank you so much. This really helped!

[rcolpo]

Hi bjodah,

Do you think it's also possible to check if a given reaction is already stoichiometrically balanced? (or whether a stoichiometric solution is correct?) I'm asking this because every reaction has multiple possible stoichiometric solutions, and I would like to test if a suggested solution is valid before trying to find another solution (of course, with repeated metabolites in both sides). Thank you.

[bjodah]

Sure, ReactionSystem checks for imbalance by default:

>>> ReactionSystem([Reaction({'H2O': 1}, {'H2': 1, 'O2': 1})],
...     substances=[Substance.from_formula(f) for f in 'H2O H2 O2'.split()])
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-4-82729cf4de88> in <module>()
----> 1 ReactionSystem([Reaction({'H2O': 1}, {'H2': 1, 'O2': 1})], substances=[Substance.from_formula(f) for f in 'H2O H2 O2'.split()])

~/vc/chempy/chempy/reactionsystem.py in __init__(self, rxns, substances, name, checks, substance_factory, missing_substances_from_keys, sort_substances)
    102 
    103         for check in (self.default_checks if checks is None else checks):
--> 104             getattr(self, 'check_'+check)(throw=True)
    105 
    106         if sort_substances:

~/vc/chempy/chempy/reactionsystem.py in check_balance(self, strict, throw)
    308                     if throw:
    309                         raise ValueError("Composition violation (%s: %s) in %s" %
--> 310                                          (k, net, rxn.string(with_param=False)))
    311                     else:
    312                         return False

ValueError: Composition violation (8: 1) in H2O -> H2 + O2

(Any improvements to the documentation or the README is of course most welcome)