Author: Björn Dahlgren.
Let us consider the reaction:
$$ Fe^{3+} + SCN^- \rightarrow FeSCN^{2+} $$the product is strongly coloured and we have experimental data (from a stopped-flow appartus) of the absorbance as function of time after mixing for several replicates. The experiment was performed at 7 different temperatures and for one temperature, 7 different ionic strengths. For each set of conditions the experiment was re-run 7 times (replicates). In this notebook, we will determine the activation enthalpy and entropy through regressions analysis, we will also look at the ionic strength dependence.
import bz2, codecs, collections, functools, itertools, json
import numpy as np
import matplotlib.pyplot as plt
import chempy
import chempy.equilibria
from chempy.electrolytes import ionic_strength
from chempy.kinetics.arrhenius import fit_arrhenius_equation
from chempy.printing import number_to_scientific_latex, as_per_substance_html_table
from chempy.properties.water_density_tanaka_2001 import water_density
from chempy.units import rescale, to_unitless, default_units as u
from chempy.util.regression import least_squares, irls, avg_params, plot_fit, plot_least_squares_fit, plot_avg_params
from chempy._solution import QuantityDict
%matplotlib inline
print(chempy.__version__)
Experimental conditions, the two solutions which were mixed in 1:1 volume ratio in a stopped flow apparatus:
sol1 = QuantityDict(u.mM, {'SCN-': 3*u.mM, 'K+': 3*u.mM, 'Na+': 33*u.mM, 'H+': 50*u.mM, 'ClO4-': (33+50)*u.mM})
sol2 = QuantityDict(u.mM, {'Fe+3': 6*u.mM, 'H+': 50*u.mM, 'ClO4-': (3*6+50)*u.mM})
sol = (sol1 + sol2)/2 # 1:1 volume ratio at mixing
sol.quantity_name = 'concentration'
Ibase = ionic_strength(rescale(sol/water_density(293*u.K, units=u), u.molal))
print(Ibase)
sol
ionic_strength_keys = 'abcd'
ionic_strengths = dict(zip(ionic_strength_keys, [0, 20, 40, 60]*u.molal*1e-3))
temperature_keys = '16.5 18.5 20.5 22.5 24.5'.split()
T0C = 273.15*u.K
temperatures = {k: T0C + float(k)*u.K for k in temperature_keys}
nrep = 7
indices = lambda k: (ionic_strength_keys.index(k[0]), temperature_keys.index(k[1]))
We will read the data from a preprocessed file:
transform = np.array([[1e-3, 0], [0, 1e-4]]) # converts 1st col: ms -> s and 2nd col to absorbance
_reader = codecs.getreader("utf-8")
_dat = {tuple(k): np.dot(np.array(v), transform) for k, v in json.load(_reader(bz2.BZ2File('specdata.json.bz2')))}
data = collections.defaultdict(list)
for (tI, tT, tR), v in _dat.items():
k = (tI, tT) # tokens for ionic strength and temperatures
data[k].append(v)
assert len(data) == len(ionic_strengths)*len(temperatures) and all(len(serie) == nrep for serie in data.values())
Let's plot the data:
def mk_subplots(nrows=1, subplots_adjust=True, **kwargs):
fig, axes = plt.subplots(nrows, len(ionic_strengths), figsize=(15,6), **kwargs)
if subplots_adjust:
plt.subplots_adjust(hspace=0.001, wspace=0.001)
return axes
def _set_axes_titles_to_ionic_strength(axes, xlim=None, xlabel=None):
for tI, ax in zip(ionic_strength_keys, axes):
ax.set_title(r'$I\ =\ %s$' % number_to_scientific_latex(Ibase + ionic_strengths[tI], fmt=3))
if xlabel is not None:
ax.set_xlabel(xlabel)
if xlim is not None:
ax.set_xlim(xlim)
def plot_series(series):
axes = mk_subplots(sharey=True, sharex=True)
colors = 'rgbmk'
for key in itertools.product(ionic_strengths, temperatures):
idx_I, idx_T = indices(key)
for serie in series[key]:
axes[idx_I].plot(serie[:, 0], serie[:, 1], c=colors[idx_T], alpha=0.15)
_set_axes_titles_to_ionic_strength(axes, xlim=[0, 3.4], xlabel='Time / s')
axes[0].set_ylabel('Absorbance')
for c, tT in zip(reversed(colors), reversed(temperature_keys)):
axes[0].plot([], [], c=c, label=tT + ' °C')
axes[0].legend(loc='best')
plot_series(data)
We see that one data series is off: 16.5 ℃ and 0.0862 molal. Let's ignore that for now and perform the fitting, let's start with a pseudo-first order assumption (poor but simple):
def fit_pseudo1(serie, ax=None):
plateau = np.mean(serie[2*serie.shape[0]//3:, 1])
y = np.log(np.clip(plateau - serie[:, 1], 1e-6, 1))
x = serie[:, 0]
# irls: Iteratively reweighted least squares
res = irls(x, y, irls.gaussian)
if ax is not None:
plot_least_squares_fit(x, y, res)
return res
beta, vcv, info = fit_pseudo1(data['a', '16.5'][0], ax=True)
def fit_all(series, fit_cb=fit_pseudo1, plot=False):
if plot:
axes = mk_subplots(nrows=len(temperatures), sharex=True, sharey=True)#, subplots_adjust=False)
_set_axes_titles_to_ionic_strength(axes[0, :])
avg = {}
for key in itertools.product(ionic_strengths, temperatures):
idx_I, idx_T = indices(key)
opt_params, cov_params = [], []
for serie in series[key]:
beta, vcv, nfo = fit_cb(serie)
opt_params.append(beta)
cov_params.append(vcv)
ax = axes[idx_T, idx_I] if plot else None
avg[key] = avg_params(opt_params, cov_params)
plot_avg_params(opt_params, cov_params, avg[key], ax=ax, flip=True, nsigma=3)
for tk, ax in zip(temperature_keys, axes[:, 0]):
ax.set_ylabel('T = %s C' % tk)
return avg
result_pseudo1 = fit_all(data, plot=True)
def pseudo_to_k2(v):
unit = 1/sol['Fe+3']/u.second
k_val = -v[0][1]*unit
k_err = v[1][1]*unit
return k_val, k_err
k_pseudo1 = {k: pseudo_to_k2(v) for k, v in result_pseudo1.items()}
k_pseudo1['a', '16.5']
k2_unit = 1/u.M/u.s
axes = mk_subplots(sharex=True, sharey=True)
for idxI, (tI, I) in enumerate(ionic_strengths.items()):
series = np.empty((len(temperatures), 3))
for idxT, (tT, T) in enumerate(temperatures.items()):
kval, kerr = [to_unitless(v, k2_unit) for v in k_pseudo1[tI, tT]]
lnk_err = (np.log(kval + kerr) - np.log(kval - kerr))/2
series[idxT, :] = to_unitless(1/T, 1/u.K), np.log(kval), 1/lnk_err**2
x, y, w = series.T
res = b, vcv, r2 = least_squares(x, y, w)
plot_least_squares_fit(x, y, res, w**-0.5, plot_cb=functools.partial(
plot_fit, ax=axes[idxI], kw_data=dict(ls='None', marker='.'), nsigma=3))
axes[idxI].get_lines()[-1].set_label(r'$E_{\rm a} = %.5g\ kJ/mol$' % (-b[1]*8.314511e-3))
axes[idxI].legend()
_set_axes_titles_to_ionic_strength(axes)