bjodah / chempy

⚗ A package useful for chemistry written in Python
BSD 2-Clause "Simplified" License
536 stars 79 forks source link

Issue with Balancing Reactions (from examples) #209

Open jacaughran-uga opened 2 years ago

jacaughran-uga commented 2 years ago

I discovered chempy this morning and was working through some of the examples so I apologize if this is something I've done wrong and isn't really an issue.

I tried to use the example code from Balancing Reactions to balance a couple of redox reactions from the end-of-chapter problems in the textbook we use. The example code worked flawlessly but when I built the reactions from the textbook I got odd coefficients: [-12, -12] in one instance and [24, 12] in the other. These should have been [-1, -1] or [1, 1] for the first and [2, 1] for the second.

Here is my code for each instance:

# balance the reaction Cr2O7-2 + I- --> Cr+3 + IO3- in acidic solution
# the answer should be Cr2O7-2 + 8 H+ + I- --> 2 Cr+3 + 4 H2O + IO3-
import chempy as chem
ox_r, ox_p = chem.balance_stoichiometry({'I-', 'H2O'}, {'IO3-', 'H+', 'e-'})
red_r, red_p = chem.balance_stoichiometry({'Cr2O7-2', 'H+', 'e-'}, {'Cr+3', 'H2O'})

print(ox_r, ox_p)
print(red_r, red_p)

ox = chem.Equilibrium(ox_r, ox_p, K1)
red = chem.Equilibrium(red_r, red_p, K2)

coeff = chem.Equilibrium.eliminate([ox, red], 'e-')
print(f'coeff: {coeff}')

redox = ox*coeff[o] + red*coeff[1]
print(redox)
# balance the reaction MnO4- + br- --> MnO2 + BrO3-
# the answer should be  2 MnO4- + Br- + 2 H+ -->  2 MnO2 + BrO3- + H2O
import chempy as chem
ox_r, ox_p = chem.balance_stoichiometry({'MnO4-', 'H+', 'e-'}, {'MnO2', 'H2O'})
red_r, red_p = chem.balance_stoichiometry({'Br-', 'H2O'}, {'BrO3-', 'H+', 'e-'})

print(ox_r, ox_p)
print(red_r, red_p)

ox = chem.Equilibrium(ox_r, ox_p, K1)
red = chem.Equilibrium(red_r, red_p, K2)

coeff = chem.Equilibrium.eliminate([ox, red], 'e-')
print(f'coeff: {coeff}')

redox = ox*coeff[0] + red*coeff[1]
print(redox)
jeremyagray commented 2 years ago

It looks like this may be a bug in how eliminate is calculating its elimination coefficients. Let me take a look later today.

I was able to replicate the behavior you are seeing and I was able to generate the correct coefficients with a small fix, but I need to verify with some more tests.

jeremyagray commented 2 years ago

So it looks like this is a bug in how the coefficients by which the equations should be multiplied are calculated. The code was using the result of sympy.primefactors() and division to get the coefficients instead of counting the factors from sympy.factorint(). The old code worked for non-repeating, relatively prime factors but was failing on others and the only test was on the former. I've managed to patch that and rework the elimination tests and add your examples. I've pushed the changes here. If you can test it, let me know how it works.

The sign of the coefficients still needs some work I believe. Everything works, signs included, except the manganese-bromate redox reaction. Right now the elimination coefficients are both negative and it seems like they should be positive but I need to chase more of the stoichiometry logic in other functions to be certain. I would also like to add some documentation on the signs of these coefficients (at least in the code; it may be elsewhere) and some more tests cases.

jeremyagray commented 2 years ago

I've pushed the latest fixes to the above branch. The permanganate-bromate reaction elimination coefficients should both be positive but were both returning negative. The important thing was that both coefficients have the same sign; the negative would only reverse the reaction. This was happening arbitrarily based on whether the first reaction was eliminating a reactant and the second a product, or vice-versa. I've added a check for all negative coefficients, and that seems to solve all the current tests.

At this point, we should probably add more tests to increase confidence in the solution and get some more eyes on the changes. Currently, the elimination method will only handle two reactions and not any number of reactions in general. So it works fine for redox half reactions but would need to be generalized further before using for something like elimination in a sequence of Hess's Law reaction steps. There's a commented example of a failing test that demonstrates the problem.