search

bjodah / chempy

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

Root finding failed #161

Open zplizzi opened 5 years ago

zplizzi commented 5 years ago

When trying to find the equilibria of the following system, I get "root finding failed" and "too much of at least one component" warnings.


eqsys = EqSystem.from_string("""

CO2 + H2O = H2CO3; 1.7e-3
H2CO3 = H+ + HCO3-; 10**-6.3
HCO3- = H+ + CO3-2; 10**-10.3

HPO4-2 + H2O = H2PO4- + OH-; 10**-6.79
Na2HPO4 = 2 Na+ + HPO4-2; 10**5
KH2PO4 = K+ + H2PO4-; 10**5

HCl = H+ + Cl-; 10**5.9
CH3CO2H = CH3CO2- + H+; 10**-4.756

H2O = H+ + OH-; 10**-14
""")

arr, info, sane = eqsys.root(defaultdict(float, \
    {'H2O': 55.4,
     'CO2': 3, 
     "KH2PO4": 0.0422 * 1,
     "Na2HPO4": 0.0578 * 1,
     "HCl": 0,
     "CH3CO2H": 0.0,
     "CH3CO2-": 0
    }))
conc = dict(zip(eqsys.substances, arr))
print(conc)

Even with somewhat simplified versions of this system, I still get the errors at times.

bjodah commented 5 years ago

Hi, it was a while since I was solving equilibria, but I remember that the solvers provided by scipy had trouble with convergence quite often. I tried applying various transformations with some success, but I think using another solver is probably a better solution long term, I started looking into using IPOPT (cyipopt provides Python bindings), but I got busy with other work.

For trying variable transformations, I have an example notebook (ammonical cupric solution), if you want to try IPOPT there will be some development that needs to be done, I won't have time to do this myself in the foreseeable future, but if you're interested in trying that avenue, let me know and I'll try to help out.

Either way, you'll probably want to indicate that some of these species belong to another phase (it's solubility products for the potassium and sodium salts I presume?). I seem to remember that the default parser respected the (s) suffix?

zplizzi commented 5 years ago

No worries, I probably don't have the time/skill to dig into improving the solver. I'll just work with what you have here - this library is awesome overall!

bjodah commented 5 years ago

I'm happy to hear that you are finding the library useful! And I appreciate that you're opening issues as well, it will make it easier to prioritize what to improve when I find some spare time.

bjodah commented 3 years ago

Just saw that openturns wraps Ipopt in their python wrapper, leaving a link here for future reference: http://openturns.github.io/openturns/latest/user_manual/_generated/openturns.Ipopt.html#ipopt