Numerically solving chemical equilibrium equations

Here are chemical equilibrium equations:

$$\begin{cases} \mathbf{N}^\text{T}\mathbf{X}+\mathbf{C}=\mathbf{Y}\\ \mathbf{N}\ln\mathbf{Y}=\ln\mathbf{K} \end{cases}$$

Here $$\mathbf{C}=\begin{pmatrix}c_1\\c_2\\\vdots\\c_n\end{pmatrix}$$ are initial concentrations (mol/L) of n types of substances;

$$\mathbf{N}=\begin{pmatrix} \nu_{11}&\nu_{12}&\cdots&\nu_{1n}\\ \nu_{21}&\nu_{22}&\cdots&\nu_{2n}\\ \vdots&\vdots&\ddots&\vdots\\ \nu_{m1}&\nu_{m2}&\cdots&\nu_{mn}\\ \end{pmatrix}$$ are coefficients of m chemical reactions (negative for reagents, positive for products, and zero if the substance is not related to the reaction);

$$\mathbf{X}=\begin{pmatrix}x_1\\x_2\\\vdots\\x_m\end{pmatrix}$$ are amounts of reactions (concentration changes for each reaction);

$$\mathbf{Y}=\begin{pmatrix}y_1\\y_2\\\vdots\\y_n\end{pmatrix}$$ are final (equilibrium state) concentrations;

$$\mathbf{K}=\begin{pmatrix}k_1\\k_2\\\vdots\\k_m\end{pmatrix}$$ are equilibrium constants for each reaction.

I need to solve X (m variables) and Y (n variables) by given C, N and K, either analytic or numeric solutions.

For Example

0.1 mol/L ammonia solution, we have 2 equilibria:

(A) $$\ce{H2O <=> H+ + OH-}$$ (concentration of $$\ce{H2O}$$ is not considered)

(B) $$\ce{NH3 + H2O <=> NH4+ + OH-}$$ (concentration of $$\ce{H2O}$$ is not considered)

Equilibrium constants:

• $$[\ce{H+}][\ce{OH-}] = 10^{-14}$$
• $$[\ce{NH4+}][\ce{OH-}]/[\ce{NH3}] = 1.77 \times 10^{-5}$$

So we have 2 reactions (A and B) and 4 substances ($$\ce{H+}$$, $$\ce{OH-}$$, $$\ce{NH3}$$, and $$\ce{NH4+}$$), and known values are:

$$\mathbf{C}=\begin{pmatrix}c_\ce{H+}\\c_\ce{OH-}\\c_\ce{NH3}\\c_\ce{NH4+}\end{pmatrix}=\begin{pmatrix}0\\0\\0.1\\0\end{pmatrix}$$

$$\mathbf{N}=\begin{pmatrix}\nu_{\text{A},\ce{H+}}&\nu_{\text{A},\ce{OH-}}&0&0\\0&\nu_{\text{B},\ce{OH-}}&\nu_{\text{B},\ce{NH3}}&\nu_{\text{B},\ce{NH4}}\end{pmatrix}=\begin{pmatrix}1&1&0&0\\0&1&-1&1\end{pmatrix}$$

$$\mathbf{K}=\begin{pmatrix}k_A\\k_B\end{pmatrix}=\begin{pmatrix}10^{-14}\\1.77\times10^{-5}\end{pmatrix}$$

The solutions are:

$$\mathbf{X}=\begin{pmatrix}x_A\\x_B\end{pmatrix}=\begin{pmatrix}7.57\times10^{-12}\\1.32\times10^{-3}\end{pmatrix}$$

$$\mathbf{Y}=\begin{pmatrix}y_\ce{H+}\\y_\ce{OH-}\\y_\ce{NH3}\\y_\ce{NH4+}\end{pmatrix}=\begin{pmatrix}7.57\times10^{-12}\\1.32\times10^{-3}\\0.0987\\1.32\times10^{-3}\end{pmatrix}$$

Here I used the approximation of $$y_\ce{OH-}=y_\ce{NH4+}$$, so just solved a quadratic equation.

But how can I solve a more complicated system like:

\begin{align} \ce{H2O &<=> H+ + OH-} \\ \ce{NH3 + H2O &<=> NH4+ + OH-} \\ \ce{CO2 + H2O &<=> H+ + HCO3-} \\ \ce{HCO3- &<=> H+ + CO3^2-} \end{align}

which includes 4 reactions and 7 substances?

• I'm trying to get numeric solution via SciPy but failed. Here is my example: stackoverflow.com/questions/68302523 Jul 11, 2021 at 14:08
• Please, consider taking a look at FAQ: How can I format math/chemistry expressions on Chemistry Stack Exchange?. Clearly, you already know a good part of it, so specifically I suggest section 3 which describes the mhchem package. Jul 11, 2021 at 14:16
• You may find it useful to make your variables be $\ln [\ce{A}]$ and $\ln [\ce{B}]$ rather than $[\ce{A}]$ and $[\ce{B}]$, which should be more robust numerically. Jul 13, 2021 at 18:32
• Smith and Missen cover this in their textbook. First of all, get rid of your set of equations. You can generate the minimum necessary set using only the species, dont get trapped into thinking like an experimentalist and thinking you need reaction paths. Jul 13, 2021 at 18:47
• @a-cyclohexane-molecule Use $\ln$ works better than original. But also need some luck to try initial values. e.g. I changed equation $\mathbf{N}\ln(\mathbf{N}^\text{T}\mathbf{X}+\mathbf{C})=\ln\mathbf{K}$ to $\mathbf{N}\ln(\mathbf{N}^\text{T}e^\mathbf{X}+\mathbf{C})=\ln\mathbf{K}$, then with inital value X=(-1,-1) and X=(-10,-10) fsolve can't converge, while X=(-3,-3) can get correct result. Jul 16, 2021 at 6:24

First, just so you know, you don't need to use a set of reactions that make sense... See Smith & Missen Chemical Reaction Equilibrium Analysis: Theory and Algorithms. It shows how, given a set of species, you can make an infinite number of set of chemical reactions that all satisfy elemental abundance etc...

Part 1: Systems of equations - Enter the Jacobian

Each one of these expressions with the evil equilibrium constant, is a function, a nonlinear function generally, in the case of reactions. For a function, we can perform a taylor expansion to linearize it, and make it workable. We then can use linear algebra to solve the set of linear equations, here is an example of two variables... your problem has 4, but it is all general, just add the extra terms similar to the first two variables, and you will have 4 equations rather than 2

This is the familiar $$Ax=b$$ where $$A$$ is the Jacobian, and is a matrix of derivatives, x is what you want... in reactions these are usually your extents of reactions, one per reaction For example, above, $$f(\xi_1, \xi_2, \xi_3, \xi_4) = K - \frac{[H^+][OH^-]}{[H_2O]} = 0$$, but put concentrations in terms of $$\xi's$$ using initial amounts, changes (extents of reactions) possibly via an ice table. The other three equations will be more important, or, at least 2 of them will.

Anyways, you need the 4 equations, one for each reaction, and the Jacobian matrix (matrix of derivatives.). I prefer using Sympy to calculate the derivatives analytically for me, and populate the matrix. Back in the day though, this was done by hand (yikes). Then we have the matrix A which is just numbers, and the vector b, again, just a vector of numbers, and just need to solve Ax=b. Solving Ax=b is trivial using software, you just need to provide a guess for what each extent of reaction is in order to evaluate the b vector (functions) and A matrix (derivatives of the 4 functions), and off it goes (a bad guess will give a bad answer, unfortunately). Solving this once will give you a new vector of extents of reactions... you repeat this, recalculate the b vector (functions) with the new extents of reactions, and the new Jacobian, solve for the new extents of reactions, and keep going until the extents of reactions stop changing.

You can also solve this using blackbox solvers in programs like Julia, Python, Matlab etc. You need to be able to provide them with the functions though, in terms of initial concentrations and extents of reactions.

To solve for x, you just need to solve $$x = A^{-1}b$$. normally you don't take the inverse, but this is so small, you can just tell the software to take the inverse of A and multiply it by b. There are many ways to solve it though, depending on the software. feel free to comment for more info you are after a Julia, python or matlab method.

I will add some symbolic python code in awhile, I am at my end though... this is not my first version :)

Also, as has been pointed out, there are certainly the possibility of numerical issues. However, this is why the method for creating systems of equations of Smith & Missen is so important, if concentrations get small, you can create a new set of reactions automatically to deal with this. Smith has done this, I am not sure if he ever published that code, I think he at least published the method.

Also, a common method is to solve these using LaGrange multipliers as well. There are many ways to solve systems of equations.

Finally, deriving everything from the chemical potential and avoiding equilibrium constants (instead using reference chemical potentials of the species) really adds a lot of robustness. Proffs teaching first year chem love using equilibrium constants because they don't need to solve anything other than toy problems. The pros don't mess around with kiddie equilibrium constants. Chemical potential or bust.

• Your Jacobian approach looks like Newton's method. However, I failed to get the result just as a-cyclohexane-molecule mentioned. See my another post stackoverflow.com/questions/68302523 Jul 16, 2021 at 6:04
• SymPy can't give solution like $(a\ln x_1 + b\ln x_2)=c$ because the result is not analytic. If I replaces with power, e.g. $x_1^a\cdot x_2^b=c$ , if the system is too complicated (Such as $\ce{H3PO4}$ + $\ce{H2CO3}$ + $\ce{NH3}$ + ...), we may encounter the degree larger than 4, which is also unable to be solved. Jul 16, 2021 at 6:13
• I would suggest sending Professor Smith an email. I dont think he has anymore copies, but he may have a pdf version he could give you Jul 19, 2021 at 0:49
• Numerical i stabilities are a real issue, but there are methods for dealing with them. The textbook covers all of this and more. Jul 19, 2021 at 0:50
• And yes i was suggesting a Newton method, for simplicity, but there are better methods . Jul 19, 2021 at 0:51

Inspired by Cesareo's answer in MathSE, I tried SciPy's basinhopping and some tips:

1. I use "Nelder-Mead" as minimization method, although very slow, because the function may be non-differentiable. I haven't tried "Powell"; "CG" sometimes doesn't work and BFGSs don't work at all.
2. Set large and gradient "bounds" to avoid negative logarithm, and this also help basinhopping .
3. Set very large iterate count (for both basinhopping and minimize) if system is very complicated, such as McIlvaine buffer.

Here is my code with SciPy:

import numpy as np
from scipy.optimize import basinhopping

'''
# H2O <-> H+ + OH-
# NH3 + H2O <-> NH4+ + OH-
# [H+, OH-, NH3, NH4+]
# 0.1M NH3
C = np.array([0, 0, 0.1, 0])
N = np.array([
[1, 1, 0, 0],
[0, 1, -1, 1]
])
K = np.array([1e-14, 1.77e-5])
X0 = np.array([0, 0])
'''

# H2O <-> H+ + OH-
# CO2 + H2O <-> H+ + HCO3-
# HCO3- <-> H+ + CO32-
# [H+, OH-, CO2, HCO3-, CO32-]
# 0.001M CO2 + 0.01M HCO3-
C = np.array([0, 0, 0.001, 0.01, 0])
N = np.array([
[1, 1, 0, 0, 0],
[1, 0, -1, 1, 0],
[1, 0, 0, -1, 1],
])
K = np.array([1e-14, 4.3e-7, 4.8e-11])
X0 = np.array([0, 0, 0])

def concentr(X):
return np.matmul(N.transpose(), X) + C

def chem_eq(X):
C1 = concentr(X)
bounds = 0
for c1 in C1:
if c1 <= 0:
bounds += (1 - c1)*1e+14
if bounds > 0:
return bounds
return np.linalg.norm(np.matmul(N, np.log(C1)) - np.log(K))

def main():
print('Inequilibrium:', chem_eq(X0))
print('Concentrations:', concentr(X0))
# niter should be so large for H3PO4
result = basinhopping(chem_eq, X0, niter = 1000, \
minimizer_kwargs = {'method': 'Nelder-Mead', 'options': {'xatol': 1e-14, 'maxiter': 1000}})
print(result)
X = result.x
print('Inequilibrium:', chem_eq(X))
print('Concentrations:', concentr(X))

if __name__ == '__main__':
main()


Result for 0.1M $$\ce{NH3}$$:

• $$[\ce{H+}] = 7.56662656\times10^{-12}$$ (pH = 11.12)
• $$[\ce{OH-}] = 1.32159291\times10^{-3}$$
• $$[\ce{NH3}] = 9.86784071\times10^{-2}$$
• $$[\ce{NH4+}] = 1.32159290\times10^{-3}$$

Result for 0.1M $$\ce{NH4HCO3}$$:

• $$[\ce{H+}] = 4.35771835\times10^{-8}$$ (pH = 7.36)
• $$[\ce{OH-}] = 2.29477887\times10^{-7}$$
• $$[\ce{CO2}] = 1.01117642\times10^{-3}$$
• $$[\ce{HCO3-}] = 9.97783306\times10^{-3}$$
• $$[\ce{CO32-}] = 1.09905218\times10^{-5}$$

btw there is an algorithm called HALTAFALL for equilibrium calculations, which may be much more effective than basinhopping. I found a Java implementation in this project