USING GDOC WITH A CHEMICAL KINETICS MODEL

The Problem

In 2004, Taylor et al.[1] observed the cyclohexadienyl radical reaction with oxygen in non-polar solvents using a liquid-phase spectrometer. The spectrometer measured the amount of UV light exiting a solution of 1,4-cyclohexadiene and a peroxide that was exposed to the pulse of an excimer laser. This pulse initiated a series of reactions which produced cyclohexadienyl radicals. Cyclohexadienyl radicals react with oxygen to ultimately form benzene and hydroperoxyl radical. This technique, known as laser-flash photolysis, both produced and measured cyclohexadienyl radical concentration over a period of 5 millionths of a second or 5 microseconds.

The reactions proceeded under a variety of different conditions. Both the concentration of oxygen and the temperature were varied. Photodetectors collected the thousands of data points. Numerical models based on the hypothetical chemistry (see Table 1) simulated the data within error. However, optimizing the kinetic parameters in the numerical model proved disastrous. No reasonable fits were found that modeled the data with confidence. Taylor et al. were not sure if the model was incorrect or if the procedure they were using to fit the data was incorrect. They needed a method to determine the best set of parameters, also known as the globally optimum set, that would fit the data to their model.


Table: Proposed kinetic model. % latex2html id marker 521 $ k_{\ref{tr:oneFour2chd}}$ and $ k_5$ were fixed; % latex2html id marker 525 $ k_{\ref{tr:chd2pchdoo}}$, % latex2html id marker 527 $ k_{\ref{tr:chd2ochdoo}}$, and % latex2html id marker 529 $ k_{\ref{tr:ochdoo2benzene}}$ were adjusted to fit the data within the stated bounds.
# Reaction $ \boldsymbol{k_{298}}$ ( $ \mathbf{M^{-1}\boldsymbol{\mu} s^{-1}}$ or $ \mathbf{\boldsymbol{\mu} s^{-1}}$)
1 $ \mathrm{(CH_3)_3CO}+ 1,4$-$ \mathrm{C_6H_8}\rightarrow c$-$ \mathrm{C_6H_7}+ \mathrm{(CH_3)_3COH}$ 53
2 $ c$-$ \mathrm{C_6H_7}+ \mathrm{O_2}\rightleftharpoons p$-$ \mathrm{C_6H_7OO}$ [1, 1200]
3 $ c$-$ \mathrm{C_6H_7}+ \mathrm{O_2}\rightleftharpoons o$-$ \mathrm{C_6H_7OO}$ [1, 1200]
4 $ o$-$ \mathrm{C_6H_7OO}\rightarrow \mathrm{C_6H_6}+ \mathrm{HO_2}$ [0.001, 100]
5 $ 2 c$-$ \mathrm{C_6H_7}\rightarrow$   Products 1200



Mathematical Description

In formulating the optimization problem, first we need convenient variables to represent the concentrations and rate constants in the kinetic model. Let the vector $ \boldsymbol{x}$ be the state variables that correspond to the species concentrations. In brackets are the physically reasonable ranges for each of the state variables.


$\displaystyle x_1$ $\displaystyle =$ $\displaystyle C_{c\mbox{-}\mathrm{C_6H_7}} \in [0, 1.4 \times 10^{-4}]$  
$\displaystyle x_2$ $\displaystyle =$ $\displaystyle C_{\mathrm{(CH_3)_3CO}} \in [0, 1.4 \times 10^{-4}]$  
$\displaystyle x_3$ $\displaystyle =$ $\displaystyle C_{1,4\mbox{-}\mathrm{C_6H_8}} \in [0, 0.4 \times 10^{-1}]$  
$\displaystyle x_4$ $\displaystyle =$ $\displaystyle C_{p\mbox{-}\mathrm{C_6H_7OO}} \in [0, 1.4 \times 10^{-4}]$  
$\displaystyle x_5$ $\displaystyle =$ $\displaystyle C_{o\mbox{-}\mathrm{C_6H_7OO}} \in [0, 1.4 \times 10^{-4}]$  

We define the vector $ \boldsymbol{q}$ to represent the known constants in the system.


$\displaystyle q_1$ $\displaystyle =$ $\displaystyle k_1 = 53$  
$\displaystyle q_2$ $\displaystyle =$ $\displaystyle k_5 = 1200$  
$\displaystyle q_3$ $\displaystyle =$ $\displaystyle C_{\mathrm{O_2}} = 0.019$  
$\displaystyle q_4$ $\displaystyle =$ $\displaystyle K_1 = 46\mathrm{e}^{(6500/T-18)}$  
$\displaystyle q_5$ $\displaystyle =$ $\displaystyle K_2 = 2 q_4$  

Also, let the vector $ \boldsymbol{p}$ represent the adjustable or unknown values in the system. In the brackets are the physically reasonable bounds for the variables.


$\displaystyle p_1$ $\displaystyle =$ $\displaystyle k_2 \in [10.0, 1200.0]$  
$\displaystyle P_2$ $\displaystyle =$ $\displaystyle k_3 \in [10.0, 1200.0]$  
$\displaystyle p_3$ $\displaystyle =$ $\displaystyle k_4 \in [0.001, 40.0]$  

The data vector $ \boldsymbol{d}$ consists of the average value of the absorbance at each time-point in the vector $ \boldsymbol{t}$ and has an error estimate based on the replicate values of absorbance $ \boldsymbol{\sigma}$.

The model is formed from the rules for elementary rates in kinetics.


$\displaystyle \frac{\mathrm{d}x_1}{\mathrm{d}t}$ $\displaystyle =$ $\displaystyle q_1 x_2 x_3 - q_3 (p_1+p_2) x_1 + p_1/q_4 x_4$  
  $\displaystyle +$ $\displaystyle p_2/q_5 x_5 - q_2 x_1^2$  
$\displaystyle \frac{\mathrm{d}x_2}{\mathrm{d}t}$ $\displaystyle =$ $\displaystyle -q_1 x_2 x_3$  
$\displaystyle \frac{\mathrm{d}x_3}{\mathrm{d}t}$ $\displaystyle =$ $\displaystyle -q_1 x_2 x_3$  
$\displaystyle \frac{\mathrm{d}x_4}{\mathrm{d}t}$ $\displaystyle =$ $\displaystyle p_1 x_1 q_3 - p_1/q_4 x_4$  
$\displaystyle \frac{\mathrm{d}x_5}{\mathrm{d}t}$ $\displaystyle =$ $\displaystyle p_2 q_3 x_1 - (p_3+p_2/q_5) x_5$  

Initial conditions for the model are based on the experimental values for the reactants.


$\displaystyle {x_1}^0$ $\displaystyle =$ 0  
$\displaystyle {x_2}^0$ $\displaystyle =$ $\displaystyle 1.4 \times 10^{-4}$  
$\displaystyle {x_3}^0$ $\displaystyle =$ $\displaystyle 4.0 \times 10^{-1}$  
$\displaystyle {x_4}^0$ $\displaystyle =$ 0  
$\displaystyle {x_5}^0$ $\displaystyle =$ 0  

Finally, using the known values of the absorbance for cyclohexadienyl radical and a value for the background absorbance we can create an objective function for comparing the data at each point with the state variables in the model. Let i be the index corresponding to each data point.

$\displaystyle {\left(\frac{2100 x_1 +200 x_4 + x_5 - d_i}{\sigma_i}\right)}^2$    


Input and Data Files

With the problem formulated, we can now write an input file for the GDOC program. A data file must also be made which contains the time, absorbance, and the standard deviation of the absorbance. In order to scale the problem correctly, the standard deviation has been divided by 100. In order to get the actual $ \chi^2$ values from the output it is necessary to multiply the global optimum found by 100$ ^2$.

An input file and data set for an experiment at T = 298 K is given here:

Input file:chi298.oai (1468 bytes)
Data file: chi298.inp (14096 bytes)

An input file and data set for an experiment at T = 323 K is given here:

Input file:chi323.oai (1468 bytes)
Data file: chi323.inp (14096 bytes)

Running GDOC

Information on how to obtain GDOC can be found on the GDOC homepage under the heading How can I get it? . In order to run GDOC you will need to have gcc version 3.2.3 or higher.

When you have obtained a copy of the distribution, unpack it and cd into the GDOC-1.0 directory. 

gunzip GDOC-1.0.tar.gz
tar xvf GDOC-1.0.tar
cd GDOC-1.0

Follow the steps in the README file carefully to install the GDOC distribution. Once the distribution has been built, change into the examples directory. Copy the chi298.oai and chi323.oai files into the examples/input directory. Copy chi298.inp and chi323.inp into the examples/data directory. Replace the makefile in the examples directory with this Makefile, which includes the make commands for the chi298 and chi323 examples. If you wish to keep the old makefile rename it to makefile.old.

Once the files are copied, type 

make chi298
at the prompt. A suggested instruction for running the chi298 data set will be provided. Copy that instruction and enter it at the prompt. You should see the GDOC output sent to the screen listing the objective function value and the branch-and-bound nodes.

After a short time (approximately 20-30 s on a fast computer) the global optimum will be displayed along with the ln scaled parameters. To convert to the actual values for the rate constants simply take the exponential of the values. Within the tolerance specified the number should match the value found with the SLSQP solver.

The chi323 example is run the same way. Type 

make chi323
at the prompt and follow the instructions. However, the time required to obtain the global optimum with SLSQP and no domain reduction could be several hours. The authors recommend downloading the free IPOPT solver and using the domain reduction option set to true to speed things up. This will require re-compiling the solver and linking the IPOPT library to GDOC. The output from GDOC in verbose mode is shown in the log files below:

Output for the fit to experiment at T = 298 K. The results are at the bottom of the file:

Output file:chi298.log (108668 bytes)

Output for the fit to experiment at T = 323 K. A file containing just the results for 323 K is listed, as well as a gzipped file containing the full output file:

Results file:chi323r.log (315 bytes)
Output file:chi323.log.gz (1174139 bytes)

Values for the global optima and results using other capabilities of GDOC will soon be published in Reference 2.

References

1
James W. Taylor, Gerhard Ehlker, Hans-Heinrich Carstensen, Leah Ruslen, Robert W. Field, and William H. Green.
Direct measurement of the fast, reversible reaction of cyclohexadienyl radicals with oxygen in nonpolar solvents.
Journal of Physical Chemistry A, 108(35):7193-7203, 2004. http://dx.doi.org/10.1021/jp0379547

2
Adam B. Singer, James W. Taylor, Paul I. Barton, and William H. Green.
Global dynamic optimization for parameter estimation in chemical kinetics.
Submitted to the Journal of Physical Chemistry A


GDOC Homepage | GDOC Kinetics Example | GDOC User's Manual | License Terms and Conditions |
Web site designed and maintained by James W. Taylor, jameswtaylor@gmail.com.