/
_equilibrium.py
93 lines (80 loc) · 2.47 KB
/
_equilibrium.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
# -*- coding: utf-8 -*-
from __future__ import (absolute_import, division, print_function)
try:
import numpy as np
except ImportError:
np = None
from .chemistry import equilibrium_quotient
def equilibrium_residual(rc, c0, stoich, K, activity_product=None):
"""
Parameters
---------
rc: float
Reaction coordinate
c0: array_like of reals
concentrations
stoich: tuple
per specie stoichiometry coefficient
K: float
equilibrium constant
activity_product: callable
callback for calculating the activity product taking
concentration as single parameter.
"""
if not hasattr(stoich, 'ndim') or stoich.ndim == 1:
c = c0 + stoich*rc
else:
c = c0 + np.dot(stoich, rc)
Q = equilibrium_quotient(c, stoich)
if activity_product is not None:
Q *= activity_product(c)
return K - Q
def _get_rc_interval(stoich, c0):
""" get reaction coordinate interval """
limits = c0/stoich
if np.any(limits < 0):
upper = -np.max(limits[np.argwhere(limits < 0)])
else:
upper = 0
if np.any(limits > 0):
lower = -np.min(limits[np.argwhere(limits > 0)])
else:
lower = 0
if lower is 0 and upper is 0:
raise ValueError("0-interval")
else:
return lower, upper
def _solve_equilibrium_coord(c0, stoich, K, activity_product=None):
from scipy.optimize import brentq
mask, = np.nonzero(stoich)
stoich_m = stoich[mask]
c0_m = c0[mask]
lower, upper = _get_rc_interval(stoich_m, c0_m)
# span = upper - lower
return brentq(
equilibrium_residual,
lower, # + delta_frac*span,
upper, # - delta_frac*span,
(c0_m, stoich_m, K, activity_product)
)
def solve_equilibrium(c0, stoich, K, activity_product=None):
"""
Solve equilibrium concentrations by using scipy.optimize.brentq
Parameters
----------
c0: array_like
Initial guess of equilibrium concentrations
stoich: tuple
per specie stoichiometry coefficient (law of mass action)
K: float
equilibrium constant
activity_product: callable
see ``equilibrium_residual``
delta_frac: float
to avoid division by zero the span of searched values for
the reactions coordinate (rc) is shrunk by 2*delta_frac*span(rc)
"""
stoich = np.array(stoich)
c0 = np.array(c0)
rc = _solve_equilibrium_coord(c0, stoich, K, activity_product)
return c0 + rc*stoich