Numerically solving chemical equilibrium equations

Question

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?

Answer 1

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.

Answer 2

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