{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Chemical Kinetics and Numerical Integration\n", "\n", "Here we will use methods of numerical integration to solve for the abundances of the H$_3^+$ isotopologues in the ion trap experiment from last week's notebook. After using integrated rate equations and curve fitting, we came up with this result:\n", "\n", "![Deuteration Results](https://leeping.github.io/che155/assets/images/week05/deuteration_results.png)\n", "\n", "The deviations, most notable in the D$_2$H$^+$ results, are because the reverse reactions were not included in our model. It would be very difficult to derive new rate equations, so we will use numerical methods instead.\n", "\n", "## Forward Euler Method\n", "\n", "First, we will reimplement the exact same model as last time, but this time we will solve using the Forward Euler Method. First, load in the `deuteration.csv` file. It contains the same experimental data as last week, but the time field has been rounded and lined up so that all abundances for each molecule are given at the same time values. This will make comparisons with the numerical models easier down the road." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from matplotlib import pyplot as plt\n", "import pandas as pd\n", "\n", "df = pd.read_csv('deuteration.csv')\n", " \n", "df" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As a reminder, the model was defined by the equations:\n", "\n", "$$ \\frac{\\text{d}[\\text{H}_3^+]}{\\text{d}t} = -k_1[\\text{H}_3^+] $$\n", "\n", "$$ \\frac{\\text{d}[\\text{H}_2\\text{D}^+]}{\\text{d}t} = k_1[\\text{H}_3^+] - k_2[\\text{H}_2\\text{D}^+] $$\n", "\n", "$$ \\frac{\\text{d}[\\text{D}_2\\text{H}^+]}{\\text{d}t} = k_2[\\text{H}_2\\text{D}^+] - k_3[\\text{D}_2\\text{H}^+] $$\n", "\n", "$$ \\frac{\\text{d}[\\text{H}_3^+]}{\\text{d}t} = k_3[\\text{D}_2\\text{H}^+] $$\n", "\n", "We can express these in a simple form with the matrix equation:\n", "\n", "$$ \\begin{bmatrix} \\text{d}[\\text{H}_3^+]/\\text{d}t \\\\ \\text{d}[\\text{H}_2\\text{D}^+]/\\text{d}t \\\\ \\text{d}[\\text{D}_2\\text{H}^+]/\\text{d}t \\\\ \\text{d}[\\text{D}_3^+]/\\text{d}t \\end{bmatrix} = \\begin{bmatrix} -k_1 & 0 & 0 & 0 \\\\ k_1 & -k_2 & 0 & 0 \\\\ 0 & k_2 & -k_3 & 0 \\\\ 0 & 0 & k_3 & 0 \\end{bmatrix} \\begin{bmatrix}[\\text{H}_3^+] \\\\ [\\text{H}_2\\text{D}^+] \\\\ [\\text{D}_2\\text{H}^+] \\\\ [\\text{D}_3^+] \\end{bmatrix} $$\n", "\n", "Then, taking a time step $\\Delta t$, we can compute new concentrations:\n", "\n", "$$ \\begin{bmatrix}[\\text{H}_3^+] \\\\ [\\text{H}_2\\text{D}^+] \\\\ [\\text{D}_2\\text{H}^+] \\\\ [\\text{D}_3^+] \\end{bmatrix}_{\\,i+1} = \\begin{bmatrix}[\\text{H}_3^+] \\\\ [\\text{H}_2\\text{D}^+] \\\\ [\\text{D}_2\\text{H}^+] \\\\ [\\text{D}_3^+] \\end{bmatrix}_{\\,i} + \\begin{bmatrix} -k_1 & 0 & 0 & 0 \\\\ k_1 & -k_2 & 0 & 0 \\\\ 0 & k_2 & -k_3 & 0 \\\\ 0 & 0 & k_3 & 0 \\end{bmatrix} \\begin{bmatrix}[\\text{H}_3^+] \\\\ [\\text{H}_2\\text{D}^+] \\\\ [\\text{D}_2\\text{H}^+] \\\\ [\\text{D}_3^+] \\end{bmatrix}_{\\,i} \\Delta t$$\n", "\n", "As of Python 3.5, matrix multiplication (and other types of dot products) can be done with the `@` operator. When used with `numpy.ndarray` objects, the [`numpy.matmul`](https://numpy.org/doc/stable/reference/generated/numpy.matmul.html) function is called. In our case, we will create a 4x4 matrix called `J` and a 1D array with 4 elements called `n` to store the abundances. When we call `J@n`, it multiplies each row of `J` by the 4 elements in `n`, and adds them up. Here we use the results from the curve fitting to ideally give us similar results as last time. We will set the step size `dt` to 0.1 ms, and take 1500 steps." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#iniitialize rate constants\n", "hd=6.3e10\n", "k1=1.42e-9*hd\n", "k2=1.33e-9*hd\n", "k3=1.05e-9*hd\n", "\n", "#H3+ at t=0 is 932, H2D+, D2H+, and D3+ start at 0.\n", "n0 = np.array([932,0,0,0])\n", "\n", "#initialize an empty 4x4 matrix, and plug in k values at the right places\n", "J = np.zeros((4,4))\n", "J[0,0] = -k1\n", "J[1,1] = -k2\n", "J[2,2] = -k3\n", "J[1,0] = k1\n", "J[2,1] = k2\n", "J[3,2] = k3\n", "\n", "#this array n will be updated with the new concentrations at each step. Initialize it at n0\n", "n = n0\n", "dt = 1e-4\n", "steps = 1500\n", "\n", "#this array will keep track of the values of n at each step\n", "nt = np.zeros((steps+1,len(n0)))\n", "nt[0] = n0\n", "\n", "#take each steps, updating n at each one; store the results in the nt array\n", "for i in range(0,steps):\n", " n = n + J@n*dt\n", " nt[i+1] = n\n", " \n", "nt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can plot the results and compare with the experimental data." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig,ax = plt.subplots()\n", "t = np.linspace(0,150e-3,len(nt))\n", "\n", "ax.scatter(df['time'],df['H3+'],color='#000000',label=r'H$_3^+$')\n", "ax.scatter(df['time'],df['H2D+'],color='#ffbf00',label=r'H$_2$D$^+$')\n", "ax.scatter(df['time'],df['D2H+'],color='#022851',label=r'D$_2$H$^+$')\n", "ax.scatter(df['time'],df['D3+'],color='#c10230',label=r'D$_3^+$')\n", "\n", "ax.set_xlabel(\"Time (s)\")\n", "ax.set_ylabel(\"Number\")\n", "\n", "lines = ax.plot(t,nt)\n", "lines[0].set_color('#000000')\n", "lines[1].set_color('#ffbf00')\n", "lines[2].set_color('#022851')\n", "lines[3].set_color('#c10230')\n", "ax.set_yscale('log')\n", "ax.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that the step size is a critical parameter! If we increase the step size too much, we can get some bad results." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "n = n0\n", "dt = 5e-3\n", "steps = round(.15/dt)+1\n", "\n", "nt = np.zeros((steps+1,len(n0)))\n", "nt[0] = n0\n", "\n", "for i in range(0,steps):\n", " n = n + J@n*dt\n", " nt[i+1] = n\n", " \n", "fig,ax = plt.subplots()\n", "t = np.linspace(0,len(nt)*dt,len(nt))\n", "\n", "ax.scatter(df['time'],df['H3+'],color='#000000',label=r'H$_3^+$')\n", "ax.scatter(df['time'],df['H2D+'],color='#ffbf00',label=r'H$_2$D$^+$')\n", "ax.scatter(df['time'],df['D2H+'],color='#022851',label=r'D$_2$H$^+$')\n", "ax.scatter(df['time'],df['D3+'],color='#c10230',label=r'D$_3^+$')\n", "\n", "ax.set_xlabel(\"Time (s)\")\n", "ax.set_ylabel(\"Number\")\n", "\n", "lines = ax.plot(t,nt)\n", "lines[0].set_color('#000000')\n", "lines[1].set_color('#ffbf00')\n", "lines[2].set_color('#022851')\n", "lines[3].set_color('#c10230')\n", "ax.set_yscale('log')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Least Squares Fitting and Numerical Integration\n", "\n", "It is possible (though not very common) to implement least squares fitting together with the numerical integration in order to estimate the kinetics parameters. We'll walk through the process here. Last time we used `scipy.optimize.least_squares`, which required us to calculate the residuals vector between the model and the experimental data. When using integrated rate equations, this was straightforward because we could just plug in the time for each data point into the model and compute the model's prediction. With numerical integration; however, we do not have such a function!\n", "\n", "Instead, what we can do is save the model's outputs whenever the time matches the time at which an experimental data point is taken. If we choose time steps judiciously, we can make sure that we always sample the model at each time point needed. If we inspect the data frame, we can see that all of the time points are at a multiple of 0.1 ms." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "df" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Therefore, a time step `dt` of 0.1 ms (or some integer factor smaller) will ensure that the model samples each time point we need to compare with the experimental data. The code below checks to see if `i` (the current time in units of `dt`) is in the array `tvals`, which is the time array converted to units of dt, and if so it stores the current model abundances in a list for later use. Importantly, this is chosen such that all of the time comparisons are between integers so that we don't have to worry about issues with floating point comparisons.\n", "\n", "At the end of the clode block, `nm` is a 2D numpy array where each row is a time point and each column is the abundance of one of the ions." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "n = n0\n", "dt = 1e-4\n", "steps = 1500\n", "\n", "nmodel = []\n", "tvals = df['time'].to_numpy()/dt\n", "tvals = tvals.astype(int)\n", "\n", "# print(tvals)\n", "\n", "for i in range(0,steps+1):\n", " n = n + J@n*dt\n", " if i in tvals:\n", " nmodel.append(n)\n", " \n", "nm = np.array(nmodel)\n", "nm" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "tvals" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we'll plot the results. A quick side note here: we've been doing a lot of repetitive manual color changing. If you have a set of colors you want to consistently use, you can change matplotlib's default color cycling (see this [tutorial](https://matplotlib.org/tutorials/intermediate/color_cycle.html) for a quick example). Below I create a new `cycler` object that tells matplotlib to cycle between the 4 colors we have been using instead of its defaults. As the tutorial shows, you can either set the cycler on an `Axes` object like in the code below, which only affects that object, or you can apply the cycler to all subsequently created plots." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from cycler import cycler\n", "\n", "ucd_cycler = (cycler(color=['#000000','#ffbf00','#022851','#c10230','#266041','#8a532f']))\n", "\n", "fig,ax = plt.subplots()\n", "ax.set_prop_cycle(ucd_cycler)\n", "ax.plot(df['time'],nm,'o')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's turn that into a function that takes the kinetics parameters (`h30`, `k1`, `k2`, `k3`) as arguments. We also need to pass the time values at which the model should be sampled, the step size, and the number of steps." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def runmodel(params,tvals,dt,steps):\n", " h30 = params[0]\n", " k1 = params[1]\n", " k2 = params[2]\n", " k3 = params[3]\n", " n = np.asarray([h30,0,0,0])\n", " nmodel = []\n", " J = np.zeros((4,4))\n", " J[0,0] = -k1\n", " J[1,1] = -k2\n", " J[2,2] = -k3\n", " J[1,0] = k1\n", " J[2,1] = k2\n", " J[3,2] = k3\n", "\n", " for i in range(0,steps+1):\n", " n = n + J@n*dt\n", " if i in tvals:\n", " nmodel.append(n)\n", " \n", " return(np.array(nmodel))\n", " \n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Test to make sure the `runmodel` function works as intended:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "tvals = df['time'].to_numpy()/dt\n", "tvals = tvals.astype(int)\n", "hd=6.3e10\n", "k1=1.43e-9*hd\n", "k2=1.33e-9*hd\n", "k3=1.05e-9*hd\n", "h30 = 932\n", "\n", "runmodel(np.array([h30,k1,k2,k3]),tvals,1e-4,1500)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To perform the `least_squares` optimization, we need to create a function that computes the residuals of the model. This function must have the signature `f(x,*args,**kwargs)` where `x` is an array containing the parameters that will be optimized (`h30`, `k1`, `k2`, and `k3`), `*args` contains any additional arguments that are needed, and `**kwargs` can contain any other information.\n", "\n", "Like last time, we'll use `**kwargs` to pass in the experimental data. `*args` will contain the `tvals`, `dt`, and `steps` parameters that need to be passed to `runmodel.` Ance we have the results of the model, we need to compute the residuals." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def total_fit(x,*args,**kwargs):\n", " \n", " df = kwargs['df']\n", "\n", " nm = runmodel(x,*args)\n", " \n", " #a naive algorithm using for loops; slow, but flexible!\n", "# out = []\n", "# for i,model in enumerate(nm):\n", "# for j,mol in enumerate(['H3+','H2D+','D2H+','D3+']):\n", "# n = df.at[i,mol]\n", "# if np.isfinite(n):\n", "# out.append(n-model[j])\n", "# return out\n", " \n", " #taking advantage of numpy's array routines: fast, but requires more work if anything changes\n", " rh3 = df['H3+'] - nm[:,0]\n", " rh3 = rh3[~np.isnan(rh3)] #remove NaNs... isnan returns an array of booleans, so we take the logical not and use it as a slice to extract only the finite values\n", " \n", " rh2d = df['H2D+'] - nm[:,1]\n", " rh2d = rh2d[~np.isnan(rh2d)]\n", " \n", " #there are no NaNs in the experimental data for D2H+ or D3+\n", " rd2h = df['D2H+'] - nm[:,2]\n", " rd3 = df['D3+'] - nm[:,3]\n", " \n", " #concatenate and return\n", " return np.concatenate((rh3,rh2d,rd2h,rd3))\n", " " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "total_fit(np.array([932, 1e-9, 1e-9, 1e-9]), tvals, dt, steps, df=df)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can use `least_squares` to compute optimal parameters, and we can see that we get almost exactly the same results as the integrated rate equation approach. Note, however, that there is no problem with us starting out with `k1` and `k2` being equal! There is no divide by 0 error with numerical integration like there was with the integrated rate equations." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import scipy.optimize as opt\n", "import numpy.linalg as la\n", "\n", "data = {\n", " 'df' : df\n", "}\n", "\n", "tvals = df['time'].to_numpy()/dt\n", "tvals = tvals.astype(int)\n", "hd=6.3e10\n", "\n", "result = opt.least_squares(total_fit,[900,1e-9*hd,1e-9*hd,1e-9*hd],\n", " args=[tvals,1e-4,1500],kwargs=data,verbose=1)\n", "pcov = la.inv(result.jac.T @ result.jac)\n", "\n", "for i,x in enumerate(['[H3+]0','k1','k2','k3']):\n", " den = hd\n", " if i==0:\n", " den = 1.\n", " print(f'{x} = {result.x[i]/den:.2e} +/- {np.sqrt(pcov[i][i])/den:.2e}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Integration with `scipy.integrate`\n", "\n", "Our manual implementation of the numerical integration uned the Forward Euler Method, whose total error is proportional to $(\\Delta t)^{1}$. It is usually desirable to use a higher-order method to achieve either higher accuracy or obtain the same accuracy with fewer steps. The function we are going to explore is [`scipy.integrate.solve_ivp`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html), which is made to solve initial value problems." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "import scipy.integrate as spi\n", "\n", "spi.solve_ivp?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we can see from the function description, we need to provide at least 3 arguments:\n", "- `fun` is a function that computes the vector of derivatives. Its function signature needs to be `f(t,y,*args)`. `t` is the current time, `y` is the current state array (in our case, the array containing the molecule abundances), and the remainder of the arguments can contain anything else needed to compute the derivatives (e.g., rate coefficients, etc)\n", "- `t_span` is a tuple that specifies the initial and final time for the integration\n", "- `y0` is a vector containing the initial conditions - the starting abundances for the molecules.\n", "\n", "In addition to those required parameters, there are three other optional arguments that are useful for us:\n", "- `method` selects which numerical integration method will be employed. The default, `'RK45'`, is the fourth-order Runge-Kutta method, but several others are available, including some implicit solvers that are important when problems are \"stiff.\" A system of equations is stiff when the solutions are very sensitive to the step size even when the solution appears \"smooth.\" Chemical kinetics problems are frequently stiff when there are some very slow reactions combined with others that are very fast, and you want to evaluate the system over a long time compared with the rate of the fast reactions. In the current example, all of the reactions have comparable rates, so we will stick with `'RK45'`, but often the `'Adams'` or `'Radau'` methods are more appropriate for kinetics problems.\n", "- `t_eval` is a list of times at which the model returns abundances. If this is None, the model only gives the results at the final time. If we pass an array of times, the results will contain the abundances at all of the time values specified in `t_eval` which fall within `t_span`\n", "- `dense_output` causes the solver to construct functions that interpolate between time steps. This allows you to (approximately) evaluate the model at any time, not just at the time steps that were used in the model.\n", "\n", "Note that nowhere do you need to specify the step size! All of the methods employ various algorithms to automatically determine the step size needed to bring the error down to a certain desired value. Some even include adaptive step sizes that can take smaller or larger steps depending on the magnitudes of the derivatives.\n", "\n", "Let's re-implement the same model, but this time perform the integration with `solve_ivp`. First we need to write a function that computes the derivative." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# function must take t and y as its first 2 arguments. Since our derivatives don't explicitly depend on t, that variable isn't used in the body of the function.\n", "# to calculate the rates, we need the rate coefficients and abundances. The abundances are in y, so we need to pass the k values as arguments.\n", "def calc_derivative(t,y,k1,k2,k3):\n", " \n", " J = np.zeros((len(y),len(y)))\n", " J[0,0] = -k1\n", " J[1,1] = -k2\n", " J[2,2] = -k3\n", " J[1,0] = k1\n", " J[2,1] = k2\n", " J[3,2] = k3\n", " \n", " return J@y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With that, we can now use `solve_ivp` to compute the solution from 0 to 0.15 seconds. We'll use the default `RK45` integrator, and set the `dense_output` flag to allow us to generate a quasi-continuous model function. In addition, we'll pass our `df['time']` array to `t_eval` so that we have the exact model values at the experimental time points.\n", "\n", "Within the `result` object that is returned, we can access the dense solution with `result.sol`, which takes a time value as an argument. The solution values are in `result.y`, and the time points for each solution are in `result.t`. The plot that this cell creates shows both the dense output and the discrete solutions." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hd=6.3e10\n", "k1=1.0e-9*hd\n", "k2=1.0e-9*hd\n", "k3=1.0e-9*hd\n", "h30 = 932\n", "\n", "result = spi.solve_ivp(calc_derivative,(0,.15),y0=[h30,0,0,0],\n", " t_eval=df['time'],method='RK45',\n", " dense_output=True,args=(k1,k2,k3))\n", "fig,ax = plt.subplots()\n", "ax.set_prop_cycle(ucd_cycler)\n", "ax.scatter(df['time'],df['H3+'],color='#000000',label=r'H$_3^+$')\n", "ax.scatter(df['time'],df['H2D+'],color='#ffbf00',label=r'H$_2$D$^+$')\n", "ax.scatter(df['time'],df['D2H+'],color='#022851',label=r'D$_2$H$^+$')\n", "ax.scatter(df['time'],df['D3+'],color='#c10230',label=r'D$_3^+$')\n", "t = np.linspace(0,160e-3,1000)\n", "ax.plot(t,result.sol(t).T)\n", "#ax.plot(result.t,result.y.T,'o')\n", "ax.set_yscale('log')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Using object-oriented programming to make a generalized kinetic model\n", "\n", "Up until now, our kinetic model consisted of hard-coded species and reactions. When doing research, one might want to test different models quickly, which means a more general way of setting up kinetic models is desirable. Below is a simple design of a reaction network class that stores the species in the system and individual reactions." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "class ReactionNetwork:\n", " \n", " def __init__(self):\n", " # Our constructor simply initializes some attributes to empty lists\n", " # so they can be extended by calling the other methods.\n", " self.species_list = []\n", " self.reaction_list = []\n", " \n", " def add_species(self, species):\n", " # Add a species to the list (can be \"H3+\" or \"HD\" for example)\n", " self.species_list.append(species)\n", " \n", " def add_reaction(self, reaction_str, kfwd=0.0, kbak=0.0):\n", " # Add a reaction to the list; the syntax is \"H3+ + HD -> H2D+ + H2\" \n", " # The reactants and products are separated by \" -> \"\n", " # The species on either side are separated by \"+\"\n", " reaction_split = reaction_str.split(\"->\")\n", " reactant_str = reaction_split[0]\n", " reactants = []\n", " for word in reactant_str.split(\" + \"):\n", " reactants.append(word.strip())\n", " \n", " product_str = reaction_split[1]\n", " products = []\n", " for word in product_str.split(\" + \"):\n", " products.append(word.strip())\n", " # The reaction list is stored as a list of dictionaries\n", " self.reaction_list.append({'string':reaction_str,\n", " 'reactants':reactants,\n", " 'products':products,\n", " 'kfwd':kfwd,\n", " 'kbak':kbak})\n", " \n", " def calc_derivative(self, t, concs):\n", " # Initialize the derivative to an array of zeros.\n", " derivatives = np.zeros_like(concs)\n", " \n", " # This dictionary maps the species name to the index\n", " # in the species list used to get the corresponding concentration.\n", " species_map = {}\n", " for i, species in enumerate(self.species_list):\n", " species_map[species] = i\n", " \n", " # Loop through each reaction.\n", " for reaction in self.reaction_list:\n", " # Calculate the forward rate as a product of the forward rate constant\n", " # and the reactant species concentrations.\n", " rate_fwd = reaction['kfwd']\n", " for reactant in reaction['reactants']:\n", " rconc = concs[species_map[reactant]]\n", " rate_fwd *= rconc\n", " # Do the same for the backward rate.\n", " rate_bak = reaction['kbak']\n", " for product in reaction['products']:\n", " pconc = concs[species_map[product]]\n", " rate_bak *= pconc\n", " # The net reaction rate.\n", " tot_rate = rate_fwd - rate_bak\n", " # Convert net reaction rate to rate of change of\n", " # individual concentrations.\n", " for reactant in reaction['reactants']:\n", " derivatives[species_map[reactant]] -= tot_rate\n", " for product in reaction['products']:\n", " derivatives[species_map[product]] += tot_rate\n", " return derivatives\n", " \n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With the class defined as above, we can now add our species and reactions, which can now go in both directions. The `calc_derivative` method is defined with time as the first argument (even though the derivative has no explicit time dependence) because that's how `scipy.integrate.solve_ivp` needs the function to be defined." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "myNetwork = ReactionNetwork()\n", "myNetwork.add_species(\"H3+\")\n", "myNetwork.add_species(\"H2D+\")\n", "myNetwork.add_species(\"D2H+\")\n", "myNetwork.add_species(\"D3+\")\n", "myNetwork.add_species(\"HD\")\n", "myNetwork.add_species(\"H2\")\n", "myNetwork.add_reaction(\"H3+ + HD -> H2D+ + H2\", kfwd=1e-9, kbak=1e-10)\n", "myNetwork.add_reaction(\"H2D+ + HD -> D2H+ + H2\", kfwd=1e-9, kbak=1e-10)\n", "myNetwork.add_reaction(\"D2H+ + HD -> D3+ + H2\", kfwd=1e-9, kbak=1e-10)\n", "myNetwork.calc_derivative(0, np.array([932, 0, 0, 0, 6.3e10, 0.0]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now call `scipy.integrate.solve_ivp` to integrate the kinetic model using our method of calculating the time derivative. The other arguments to `solve_ivp()` are the time interval, the initial values of dependent variables (i.e. concentrations), and the time points to evaluate and store the dependent variables." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import scipy.integrate\n", "ts = np.arange(0, 0.15, 1e-4)\n", "result = scipy.integrate.solve_ivp(myNetwork.calc_derivative, (0, 0.15), np.array([932, 0, 0, 0, 6.3e10, 0.0]), t_eval=ts)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The simulation results can now be plotted:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "ax.set_prop_cycle(ucd_cycler)\n", "ax.plot(result.t, result.y[:4].T)\n", "ax.set_yscale('log')\n", "ax.set_ylim(0.01, 1000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using the new reaction network, you can now fit the model parameters to get better agreement with experiment. Write a new function below that calculates the residuals using the new model. (The cell below will be updated with the answer after the due date.)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Fill in this cell with your own code. You may reuse any of the codes from above." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following codes will optimize the model parameters from a starting guess, assuming that your function is titled `total_fit_2` and its argument is an array of parameters consisting of the $H_3^+$ concentration, the forward reaction rates, the backward reaction rates, and the $H_2$ concentration as a percentage of the $HD$ concentration, in that order." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "param0 = np.array([932, 1e-9, 1e-9, 1e-9, 1e-10, 1e-10, 1e-10, 0.01])\n", "result = opt.least_squares(total_fit_2,param0,bounds=(0.01*param0, 100*param0),x_scale=param0,verbose=1)\n", "print(result.x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can now make a plot that shows the new model with optimized parameters gives better agreement with experiment compared to before." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# The code below is an INCOMPLETE solution to Assignment 4. \n", "# You may use it as a guide to completing the assignment but keep in mind\n", "# that it does not give the right answer.\n", "\n", "def calc_residuals(params):\n", " myNetwork = ReactionNetwork()\n", " myNetwork.add_species(\"H3+\")\n", " myNetwork.add_species(\"H2D+\")\n", " myNetwork.add_species(\"D2H+\")\n", " myNetwork.add_species(\"D3+\")\n", " myNetwork.add_species(\"HD\")\n", " myNetwork.add_species(\"H2\")\n", " myNetwork.add_reaction(\"H3+ + HD -> H2D+ + H2\", kfwd=params[1], kbak=params[4])\n", " myNetwork.add_reaction(\"H2D+ + HD -> D2H+ + H2\", kfwd=params[2], kbak=params[5])\n", " myNetwork.add_reaction(\"D2H+ + HD -> D3+ + H2\", kfwd=params[3], kbak=params[6])\n", " result = scipy.integrate.solve_ivp(myNetwork.calc_derivative, (0, 0.15), np.array([932, 0, 0, 0, 6.3e10*(1-params[7]), 6.3e10*params[7]]), t_eval=df['time'])\n", " residual_h3 = result.y[0] - df['H3+']\n", " residual_h3 = np.array(residual_h3[~np.isnan(residual_h3)])\n", " residual_h2d = result.y[1] - df['H2D+']\n", " residual_h2d = np.array(residual_h2d[~np.isnan(residual_h2d)])\n", " all_residuals = np.hstack((residual_h3, residual_h2d))\n", " return all_residuals\n", "\n", "param0 = np.array([932, 1e-9, 1e-9, 1e-9, 1e-10, 1e-10, 1e-10, 0.01])\n", "\n", "calc_residuals(param0)\n", "\n", "result = opt.least_squares(calc_residuals,param0,jac='3-point',bounds=(0.01*param0, 100*param0),x_scale=param0,verbose=1)\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "result.x" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.18" } }, "nbformat": 4, "nbformat_minor": 4 }