Incremental isoconversional method

This a simple implementation of the incremental isoconversional method described in our paper (see ref. 1 in the list of related publications). It is a 32-bit command-line program written in Fortran 90. The program performs orthogonal distance fitting of isoconversional data to a regression model defined as $$\beta=\int_{T_{i}}^{T_{i+1}}\frac{dT}{A\cdot\exp(E/RT)},$$ where $\beta$ is the (constant) heating rate, $T_i$ and $T_{i+1}$ are the temperatures corresponding to conversion levels $\alpha_i$ and $\alpha_{i+1}$ respectively. There are two parameters estimated — the activation energy (E) and the preexponential factor combined with the conversion function as $$A=\frac{F(\alpha_{i+1})-F(\alpha_i)}{A_k}.$$ Here, $A_k$ it the "real" preexponential factor of the Arrhenius equation $k(T)=A_k\exp(-E/RT)$ and $F$ is the antiderivate to the conversion function $f(\alpha)$.

Alternatively, the data can also be fitted to a model based on non-arrhenian temperature function (also known as the Berthelot—Hood function): $$k(T)=A\exp(DT)$$ In this case, the regression model is defined as $$T_{i+1} = \frac{1}{D}\ln\left[AD\beta+\exp(DT_i)\right]$$ where $A$ and $D$ are the non-arrhenian kinetic parameters.

Input file

The program reads the data from an input file named input.txt. The input file should be formatted as follows:

   	0.0250	0.0500	0.0750	...
   4.	158.33	202.30	225.72	...
   6.	163.93	207.09	230.41	...
   8.	168.68	211.92	235.09	...
   12.	173.58	216.85	240.25	...
   16.	177.41	220.58	243.62	...
   ...	...	...	...	...

The first line contains the conversion levels at which the isoconversional temperatures were evaluated. The values are not used for any calculation; they are only supplied for easier navigation in the resulting kinetic parameters. The first column contains the heating rates (in °C/min) at which the conversion was measured. The rest of the entries are the isoconversional temperatures (in °C) evaluated at corresponding conversion level and heating rate.

Parameter initialization

The optimization routine is an iterative process based on the Newton method. The program requires starting values of the kinetic parameters. This is done automatically using linear regression. Applying the trapezoidal rule to the regression model above we obtain: $$\ln\left(\frac{\beta}{T_{i+1}-T_i}\right)=\ln\frac{1}{A}-\frac{E}{R}\cdot\frac{2}{T_i+T_{i+1}},$$ which can be handled as a simple linear regression model.

Output

The program automatically calculates the parameters for all conversion levels supplied in the input file. The on-screen output should read:

 1   0.6403E-14   145481.41   0.8386912350E+04  0.61E+04
 2   0.7320E-15   144535.43   0.2766101674E+03  0.97E+03
 3   0.4262E-15   144501.70   0.1590469265E+01  0.50E+02
 4   0.2256E-14   137802.00   0.5759198034E+00  0.51E+01
 5   0.1504E-14   139410.34   0.5405388440E+00  0.16E+00
 6   0.1479E-14   139476.90   0.5404884149E+00  0.22E-03
 7   0.1479E-14   139477.00   0.5404884147E+00  0.18E-07

 ***********************************************************
 Conversion =  0.02500 --  0.05000
 A =       0.1479E-14  0.16E-14
 E =         139477.0    4312.1
 ***********************************************************

The table of iterations contains several columns, from left to right: iteration counter, A (min), E (J/mol), total sum of squares (the objective function), and its current gradient ($L_2$-norm). The last three lines contain the conversion range in which the parameters were calculated and, finally, the optimized values of parameters with their standard errors. In addition to the on-screen output, a simplified output file named output.txt is created. It should have the following structure:

  Conversion range        A         S.E.    E (kJ/mol)  S.E.
0.02500 --  0.05000   0.1479E-14   0.16E-14  139.48     4.31
0.05000 --  0.07500   0.7654E-16   0.89E-16  157.99     4.85
0.07500 --  0.10000   0.3823E-17   0.49E-17  175.11     5.54
....... --  .......   ..........   ........  ......     ....