Two coupled equilibria: protolysis of ammonia in water

In this notebook we will look at how ChemPy can be used to formulate a system of (non-linear) equations from conservation laws and equilibrium equations. We will look att ammonia since it is a fairly well-known subtance.

In [1]:
from __future__ import division, print_function
from operator import mul
from functools import reduce
from itertools import product
import chempy
from chempy.chemistry import Species, Equilibrium
from chempy.equilibria import EqSystem, NumSysLog, NumSysLin
import numpy as np
import sympy as sp
sp.init_printing()
import matplotlib.pyplot as plt
%matplotlib inline
sp.__version__, chempy.__version__
Out[1]:
('1.4', '0.8.0.dev0+git')

First, let's define our substances. ChemPy can parse chemical formulae into chempy.Substance instances with information on charge, composition, molar mass (and pretty-printing):

In [2]:
substance_names = ['H+', 'OH-', 'NH4+', 'NH3', 'H2O']
subst = {n: Species.from_formula(n) for n in substance_names}
assert [subst[n].charge for n in substance_names] == [1, -1, 1, 0, 0], "Charges of substances"
print(u'Composition of %s: %s' % (subst['NH3'].unicode_name, subst['NH3'].composition))
Composition of NH₃: {7: 1, 1: 3}

Let's define som e initial concentrations and governing equilibria:

In [3]:
init_conc = {'H+': 1e-7, 'OH-': 1e-7, 'NH4+': 1e-7, 'NH3': 1.0, 'H2O': 55.5}
x0 = [init_conc[k] for k in substance_names]
H2O_c = init_conc['H2O']
w_autop = Equilibrium({'H2O': 1}, {'H+': 1, 'OH-': 1}, 10**-14/H2O_c)
NH4p_pr = Equilibrium({'NH4+': 1}, {'H+': 1, 'NH3': 1}, 10**-9.26)
equilibria = w_autop, NH4p_pr
[(k, init_conc[k]) for k in substance_names]
Out[3]:
[('H+', 1e-07), ('OH-', 1e-07), ('NH4+', 1e-07), ('NH3', 1.0), ('H2O', 55.5)]
In [4]:
eqsys = EqSystem(equilibria, subst)
eqsys
Out[4]:
H2OH+ + OH- 1.8018⋅10-16
NH4+H+ + NH3 5.4954⋅10-10

We can solve the non-linear system of equations by calling the root method (the underlying representation uses pyneqsys:

In [5]:
x, sol, sane = eqsys.root(init_conc)
x, sol['success'], sane
Out[5]:
(array([2.34917141e-12, 5.54957436e+01, 9.95743506e-01, 4.25659361e-03,
        4.25649361e-03]), True, True)

This system is quite easy to solve for, if we have convergence problems we can try to solve a transformed system. As an example we will use NumSysLog:

In [6]:
logx, logsol, sane = eqsys.root(init_conc, NumSys=(NumSysLog,))
logx, logsol['success'], sane
Out[6]:
(array([2.34917141e-12, 5.54957436e+01, 9.95743506e-01, 4.25659361e-03,
        4.25649361e-03]), True, True)

In thise case they give essentially the same answer:

In [7]:
x - logx
Out[7]:
array([0., 0., 0., 0., 0.])

We can create symbolic representations of these systems of equations:

In [8]:
ny = len(substance_names)
y = sp.symarray('y', ny)
i = sp.symarray('i', ny)
K = Kw, Ka = sp.symbols('K_w K_a')
w_autop.param = Kw
NH4p_pr.param = Ka
ss = sp.symarray('s', ny)
ms = sp.symarray('m', ny)
In [9]:
numsys_log = NumSysLog(eqsys, backend=sp)
f = numsys_log.f(y, list(i)+list(K))
f
Out[9]:
$\displaystyle \left[ y_{0} - y_{1} + y_{4} - \log{\left(K_{w} \right)}, \ y_{0} + y_{2} - y_{3} - \log{\left(K_{a} \right)}, \ - i_{0} - i_{3} + i_{4} + e^{y_{0}} + e^{y_{3}} - e^{y_{4}}, \ - i_{0} - 2 i_{1} - 3 i_{2} - 4 i_{3} - i_{4} + e^{y_{0}} + 2 e^{y_{1}} + 3 e^{y_{2}} + 4 e^{y_{3}} + e^{y_{4}}, \ - i_{2} - i_{3} + e^{y_{2}} + e^{y_{3}}, \ - i_{1} - i_{4} + e^{y_{1}} + e^{y_{4}}\right]$
In [10]:
numsys_lin = NumSysLin(eqsys, backend=sp)
numsys_lin.f(y, i)
Out[10]:
$\displaystyle \left[ - i_{0} - i_{3} + i_{4} + y_{0} + y_{3} - y_{4}, \ - i_{0} - 2 i_{1} - 3 i_{2} - 4 i_{3} - i_{4} + y_{0} + 2 y_{1} + 3 y_{2} + 4 y_{3} + y_{4}, \ - i_{2} - i_{3} + y_{2} + y_{3}, \ - i_{1} - i_{4} + y_{1} + y_{4}\right]$
In [11]:
A, ks = eqsys.stoichs_constants(False, backend=sp)
[reduce(mul, [b**e for b, e in zip(y, row)]) for row in A]
Out[11]:
$\displaystyle \left[ \frac{y_{0} y_{4}}{y_{1}}, \ \frac{y_{0} y_{2}}{y_{3}}\right]$
In [12]:
from pyneqsys.symbolic import SymbolicSys
subs = list(zip(i, x0)) + [(Kw, 10**-14), (Ka, 10**-9.26)]
numf = [_.subs(subs) for _ in f]
neqs = SymbolicSys(list(y), numf)
neqs.solve([0, 0, 0, 0, 0], solver='scipy')
Out[12]:
(array([-3.62343921e+01, -1.33889952e-07, -2.33889967e-07, -1.49124544e+01,
         3.99820071e+00]),
    cov_x: array([[ 2.42967116e+00,  1.41467839e+00, -6.03670487e-01,
          1.82600034e+00, -1.49931131e-02],
        [ 1.41467839e+00,  1.39999983e+00, -6.00000425e-01,
          8.14677631e-01, -1.46788922e-02],
        [-6.03670487e-01, -6.00000425e-01,  4.00000323e-01,
         -2.03670497e-01,  3.66972815e-03],
        [ 1.82600034e+00,  8.14677631e-01, -2.03670497e-01,
          2.62232916e+00, -1.13233789e-02],
        [-1.49931131e-02, -1.46788922e-02,  3.66972815e-03,
         -1.13233789e-02,  3.14227056e-04]])
     fjac: array([[-9.44020661e+01,  0.00000000e+00, -5.77317876e-01,
          5.77317876e-01,  0.00000000e+00,  5.77317876e-01],
        [-1.73195322e+00, -2.82848605e+00,  3.49800834e-01,
          7.10836992e-01,  3.53545942e-01, -3.49800834e-01],
        [-1.72136041e+00, -1.06724305e+00, -1.37764628e+00,
          7.45786834e-01, -2.02347779e-01,  6.26954537e-01],
        [-1.05929885e-02, -3.47059671e-01,  1.00797381e+00,
          9.29207280e-01,  5.32982258e-01, -2.31503941e-01],
        [-5.78281895e-07,  3.53544844e-01, -2.73887397e-01,
         -6.47032732e-01, -6.17527635e-01, -9.49498512e-01]])
      fun: array([ 0.00000000e+00,  0.00000000e+00, -1.42108547e-14,  1.42108547e-14,
         1.33226763e-15,  1.42108547e-14])
     ipvt: array([5, 3, 2, 1, 4], dtype=int32)
  message: 'The relative error between two consecutive iterates is at most 0.000000'
     nfev: 16
     njev: 13
      qtf: array([-3.65093557e-12, -2.24978587e-12, -1.35957930e-12,  6.04379901e-13,
        -1.31254525e-12])
   status: 2
  success: True
        x: array([-3.62343921e+01, -1.33889952e-07, -2.33889967e-07, -1.49124544e+01,
         3.99820071e+00]))
In [13]:
j = sp.Matrix(1, len(f), lambda _, q: f[q]).jacobian(y)
init_conc_j = {'H+': 1e-10, 'OH-': 1e-7, 'NH4+': 1e-7, 'NH3': 1.0, 'H2O': 55.5}
xj = eqsys.as_per_substance_array(init_conc_j)
jarr = np.array(j.subs(dict(zip(y, xj))).subs({Kw: 1e-14, Ka: 10**-9.26}).subs(
            dict(zip(i, xj))))
jarr = np.asarray(jarr, dtype=np.float64)
np.log10(np.linalg.cond(jarr))
Out[13]:
$\displaystyle 24.404551922433335$
In [14]:
j.simplify()
j
Out[14]:
$\displaystyle \left[\begin{matrix}1 & -1 & 0 & 0 & 1\\1 & 0 & 1 & -1 & 0\\e^{y_{0}} & 0 & 0 & e^{y_{3}} & - e^{y_{4}}\\e^{y_{0}} & 2 e^{y_{1}} & 3 e^{y_{2}} & 4 e^{y_{3}} & e^{y_{4}}\\0 & 0 & e^{y_{2}} & e^{y_{3}} & 0\\0 & e^{y_{1}} & 0 & 0 & e^{y_{4}}\end{matrix}\right]$
In [15]:
eqsys.composition_balance_vectors()
Out[15]:
$\displaystyle \left( \left[ \left[ 1, \ 0, \ 0, \ 1, \ -1\right], \ \left[ 1, \ 2, \ 3, \ 4, \ 1\right], \ \left[ 0, \ 0, \ 1, \ 1, \ 0\right], \ \left[ 0, \ 1, \ 0, \ 0, \ 1\right]\right], \ \left[ 0, \ 1, \ 7, \ 8\right]\right)$
In [16]:
numsys_rref_log = NumSysLog(eqsys, True, True, backend=sp)
numsys_rref_log.f(y, list(i)+list(K))
Out[16]:
$\displaystyle \left[ y_{0} + y_{2} - y_{3} - \log{\left(K_{a} \right)}, \ y_{1} + y_{2} - y_{3} - y_{4} - \log{\left(\frac{K_{a}}{K_{w}} \right)}, \ - i_{0} - i_{3} + i_{4} + e^{y_{0}} + e^{y_{3}} - e^{y_{4}}, \ - i_{1} - i_{4} + e^{y_{1}} + e^{y_{4}}, \ - i_{2} - i_{3} + e^{y_{2}} + e^{y_{3}}\right]$
In [17]:
np.set_printoptions(4, linewidth=120)
scaling = 1e8
for rxn in eqsys.rxns:
    rxn.param = rxn.param.subs({Kw: 1e-14, Ka: 10**-9.26})
In [18]:
x, res, sane = eqsys.root(init_conc, rref_equil=True, rref_preserv=True)
x, res['success'], sane
Out[18]:
(array([1.7739e-11, 5.5469e+01, 9.6873e-01, 3.1270e-02, 3.1270e-02]),
 True,
 True)
In [19]:
x, res, sane = eqsys.root(init_conc, x0=eqsys.as_per_substance_array(
        {'H+': 1e-11, 'OH-': 1e-3, 'NH4+': 1e-3, 'NH3': 1.0, 'H2O': 55.5}))
res['success'], sane
Out[19]:
(True, True)
In [20]:
x, res, sane = eqsys.root(init_conc, x0=eqsys.as_per_substance_array(
        {'H+': 1.7e-11, 'OH-': 3e-2, 'NH4+': 3e-2, 'NH3': 0.97, 'H2O': 55.5}))
res['success'], sane
Out[20]:
(True, True)
In [21]:
init_conc
Out[21]:
{'H+': 1e-07, 'OH-': 1e-07, 'NH4+': 1e-07, 'NH3': 1.0, 'H2O': 55.5}
In [22]:
nc=60
Hp_0 = np.logspace(-3, 0, nc)
def plot_rref(**kwargs):
    fig, axes = plt.subplots(2, 2, figsize=(16, 6), subplot_kw=dict(xscale='log', yscale='log'))
    return [eqsys.roots(init_conc, Hp_0, 'H+', plot_kwargs={'ax': axes.flat[i]}, rref_equil=e,
                        rref_preserv=p, **kwargs) for i, (e, p) in enumerate(product(*[[False, True]]*2))]
In [23]:
res_lin = plot_rref(method='lm')
[all(_[2]) for _ in res_lin]
No handles with labels found to put in legend.
No handles with labels found to put in legend.
No handles with labels found to put in legend.
Out[23]:
[True, True, True, True]
In [24]:
for col_id in range(len(substance_names)):
    for i in range(1, 4):
        plt.subplot(1, 3, i, xscale='log')
        plt.gca().set_yscale('symlog', linthreshy=1e-14)
        plt.plot(Hp_0, res_lin[0][0][:, col_id] - res_lin[i][0][:, col_id])
plt.tight_layout()
/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:3: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
  This is separate from the ipykernel package so we can avoid doing imports until
In [25]:
eqsys.plot_errors(res_lin[0][0], init_conc, Hp_0, 'H+')
In [26]:
init_conc, eqsys.ns
Out[26]:
({'H+': 1e-07, 'OH-': 1e-07, 'NH4+': 1e-07, 'NH3': 1.0, 'H2O': 55.5}, 5)
In [27]:
res_log = plot_rref(NumSys=NumSysLog)
No handles with labels found to put in legend.
No handles with labels found to put in legend.
No handles with labels found to put in legend.
In [28]:
eqsys.plot_errors(res_log[0][0], init_conc, Hp_0, 'H+')
In [29]:
res_log_lin = plot_rref(NumSys=(NumSysLog, NumSysLin))
eqsys.plot_errors(res_log_lin[0][0], init_conc, Hp_0, 'H+')
No handles with labels found to put in legend.
No handles with labels found to put in legend.
No handles with labels found to put in legend.
In [30]:
from chempy.equilibria import NumSysSquare
res_log_sq = plot_rref(NumSys=(NumSysLog, NumSysSquare))
eqsys.plot_errors(res_log_sq[0][0], init_conc, Hp_0, 'H+')
No handles with labels found to put in legend.
No handles with labels found to put in legend.
No handles with labels found to put in legend.
In [31]:
res_sq = plot_rref(NumSys=(NumSysSquare,), method='lm')
eqsys.plot_errors(res_sq[0][0], init_conc, Hp_0, 'H+')
No handles with labels found to put in legend.
No handles with labels found to put in legend.
No handles with labels found to put in legend.
In [32]:
x, res, sane = eqsys.root(x0, NumSys=NumSysLog, rref_equil=True, rref_preserv=True)
x, res['success'], sane
/drone/src/github.com/bjodah/chempy/chempy/equilibria.py:260: UserWarning: Root finding indicated as failed by solver.
  warnings.warn("Root finding indicated as failed by solver.")
Out[32]:
(array([1.8480e-16, 1.0072e+00, 4.5543e-17, 1.5315e-23, 5.4500e+01]),
 False,
 True)
In [33]:
x, res, sane = eqsys.root(x, NumSys=NumSysLin, rref_equil=True, rref_preserv=True)
x, res['success'], sane
Out[33]:
(array([1.8480e-16, 1.0072e+00, 4.5543e-17, 1.5315e-23, 5.4500e+01]),
 True,
 True)