Talanta, 1971, Vol. 18, pp. 1137 to 1155. PcroamooPress. Printedin NorthernIreland
A GENERAL COMPUTER APPROACH FOR CALCULATING RATE CONSTANTS FROM NEAR-EQUILIBRIUM KINETIC STUDIES
v. s. SHARMA and Department
D. L. LEUSSING
of Chemistry, The Ohio State University,
Columbus, Ohio 43210, U.S.A.
(Received 15 October 1970. Accepted 30 November 1970) Smnmary-A general computer approach for estimating rate constants from relaxation times is described. The programme CORNEK is essentially a least-squares retkement programme applied to non-linear systems. It uses directly the differential forms of the first derivatives of mass-balance and rate equations, thus avoiding the time-consuming derivations of near-equilibrium rate equations. The programme has been tested for binary systems such as Cu-histamine, Cu-serine, and the ternary system Cu-histamine-wine.
of relaxation techniques has resulted in the elucidation of reaction mechanisms of many fast reactions which were previously inaccessible by using classical techniques .l However, at present there is no general method of evaluating rate constants from the expressions relating relaxation time to equilibrium concentrations, equilibrium constants and rate constants. Also, the derivation of such equations from mass balance relationships involves lengthy, though not complicated, algebra.a As a result, earlier workers made various simplifying assumptions with or without justification, or tested their kinetic data for only a small number of mechanisms out of a large number of possibilities. The available methods of calculating rate constants are also clumsy when applied to more complicated systems where more than one type of metal ion or ligand exists in the system. In this paper we present a general approach for calculating rate constants from relaxation times. The computer programme CORNEK is essentially a least-squares refinement programme for nonlinear systems and is almost completely general. The unique feature of the method discussed here is that in order to make a best fit, the calculated value of the dependent variable, relaxation time, is not required. The programme has been tested for simple and mixed ligand systems. 34 It uses directly the differential form of the fist derivatives of mass balance and rate equations in order to calculate rate constants from relaxation times, thus avoiding the time-consuming derivations of near-equilibrium rate equations. THE DEVELOPMENT
THEORY
theory of relaxation spectral and the general least-squares5 approach as applied to parameter estimation have been described in detail elsewhere and therefore only the portions relevant to the present problem will be discussed here. For the sake of clarity the present discussion is in terms of a specific system, a bivalent metal ion with a co-ordination number of four (such as Cu2+) and a bidentate ligand with two ionizing groups (such as protonated glycine). The extension to more general cases will follow essentially the same line and will be apparent towards the end of the discussion. The computer programme CORNEK given in the appendix is, however, a general one and can be used for any system. The
1137
1138
V. S. SIURMAand D. L. LEUSSINQ
For aqueous solutions containing a bivalent metal ion, M, with co-ordination number four and a bidentate ligand, HsA, with two dissociable protons, the following equilibria must be considered in order to characterize the dist~bution of proton or metal complex species (charges on the species have been omitted): H-l-A”,HA
(1)
H+HA+H,A
(2)
M-i-A+MA
(3)
MA+A+MA,
(4)
In addition to these expressions, as many more may be written as logically appropriate for describing various reaction paths. For example, we may like to consider the possibility of the monoprotonated form of the ligand reacting with the metal ion. Therefore : M+HA+MA+H
(9
MA+RA+MA%+H
(6)
Similarly, any other reaction path or species can be added to this scheme. If the relaxation effect is followed by means of a coloured pH indicator, In, the appropriate equilibrium should be taken into consideration: H+In+HIn
(7)
The mass balance equations for the system under consideration are: Total metal (Z’,) = [M] + [MA] + [MA,]
(8)
Total ligand (TJ = [A] + [HA] + [HaA] + [MA] + 2[MA,]
(9
Total proton (Z’u) = [H] + [HA] + 2[H,A] + [HIn] - [OH] Total indicator (Tr,,) = [HIn] + [In]
(10) (11)
Since the protolytic reactions reach equilibrium much faster than the metal complex reactions under consideration, they can be assumed to be at equilibrium at all times. For species MA, MA, we can write the following rate expressions: -4MAl/dr=
-kJM][A]
+ 2 [MA] + k,[MA][A] - g
-kiEMIIHAl+
s
2 MN~W
a
[MAJ
k6 [MA~[Hl -I- k,IMA’l[HAl - F
5
(12)
6
and -O%l/d~
= -~,[MAI[AI+
g4 IMA~ - k,[MAl[HAl
+ % [MAJ[HJ
6
(13)
Computer
1139
approach for calculating rate constants
On differentiating (12) and (13), and substituting (S[A] . [HA]/[A] + 6[H] . [HA]/[H]) for 8[HA], and (QH] .2[H,A]/[H] + 6[A] . [H,A]/[A]) for 6[H,A] and simplifying, the differential forms of the resulting expressions are obtained: -dB[MA] dt =
--k,M + MMAI - k,
+ (-k,[A]
- k&HA]) 6[M] + (-k,
MEW + k EMAIFW 61Al 6 LAl [Al ) + 2 [MA] + k, ‘“;;BA]
I”;F
5
-
W&I) W-U+ (2 + ; [HI + k&41)O,fAl
g
5
6
+ (-
2
- ; 4
-dB[MA,] dt -
- =
2 M&I)
k,[MA] + k, r”;y)
WU - (WI
[H]) 6[MA,] a
(14)
fi[A] _ (kl [M;;y
+ k.,JHAl)4MAl
+ re + $ PI) 4MA,l
(15)
8
6
Subscripts on rate and equilibrium constants refer to the equation number, k,‘s are the rate constants for the forward reactions and KJ’s are the stepwise formation constants; parentheses represent equilibrium concentrations. For a system containing IZcomplex species, II such near-equilibrium rate equations are obtained. Similarly, the differential form of the mass balance equations (8-l 1) is: 6T, = 0 = 6[M] + 6[MA] + 6[MAJ 6T, = 0 =
DWI +;flAl + [Al dLAl 2WA + [HAI WI + ) ) ( WI + ~[MAI+ 24MA,l (17) ~[H,A] + [HAI + [HI +
ST, = 0 =
(16)
gl + WW 4Hl
[HI +
2’HzA;G
do,, = 0 = ([HInl/[In]+ 1)4hl+
IHA’) 6[A] + ($I$)
(D-W/EHIP[W
6[In]
(18)
(1%
where K, = ionic product of water. Again substitution for 6[HA], 6[H,A], G[HIn] has been made in terms of 6[H], 6[A], and 6[In], respectively, as explained earlier. Equations equivalent to (l-19), for a ternary system, are given in Appendix I. On designation of the coefficients of S[MA], 6[MAJ, 8[M], 6[A], 6[H], and 6[In] as ail, aiz, a,,, ai4, ais, and aia, respectively,
1140
v. s. SEARMAand D. L. LEUSSING
where i = 1, 2, 3, 4,5, and 6 for equations (1419),
respectively, equations (14-19)
can be written in the following matrix notation. AX=Y all
al2
al3
a14
a15
al6
a21
a22
a23
a24
a25
a28
a31
a32
a33
a34
a35
a36
a41
a42
a43
a44
a45
a46
4A]
a51
a52
a53
a54
as5
a56
4Hl
0.0
aI31
a62
a63
a64
a65
a66
Wnl
0.0
11
o-0 =
(20)
0.0
By the process of elimination involving elementary row operations, from the matrix equation (20) we can easily obtain: A’X = Y ai
a12'
0
0
0
0
a21'
a22'
0
0
0
0
a31'
a32'
a33'
0
0
0
WI
a41'
a42'
a43'
a44'
0
0
6[A]
a51'
a52'
am
a54'
a55'
0
d[Hl
a61'
a62'
a63'
a64r
'65'
a66'
6[In]
0.0
=
(21)
0.0
o-0 1
I
04
I
The first two rows now correspond to the near-equilibrium rate equations of Hammes and Steinfield. The corresponding characteristic determinant is :
all’ - ;z
a12’ a22f - 3,
a2i’
=
0
or
]A” - 111 = 0
(22)
where iz = I/T, r = relaxation time, I = the identity matrix, and
and the characteristic polynomial is :
(all’ -
jl)(a22) - A) - ali . a211 = 0
In general for a square matrix of order n corresponding
to equation (22), the char-
acteristic polynomial is given by: il” + c,&‘+--l
+ c,,_~P-~ + - * - + c,A + (-1)”
]A]” = 0
(23) where c, is the k, coefficient of the polynomial. The last term is a constant equal to (-1)” times the value of the determinant of matrix A”. Estimation of rate constants by the generalized method of least squares as applied to non-linear systems In equation (23) ct)s are related to a,,‘s of the secular determinant and are functions
Computer approach for calculatingrate constants
1141
of equilibrium concentrations of all the species, stability constants and the rate constants. In near-equilibrium kinetic studies equilibrium concentrations of various species and the corresponding formation constants are supposed to be obtained from prior knowledge of solution equilibria. A reasonable initial approximation of rate constants to be refined can also be made, and any one of the observed relaxation times can be used in order to calculate the function R where: R = A” + c,+~A”--~+ c,,_&+~ + . . * + c,l + (-1)“. The least-squares problem can now be stated mathematically
IA”I.
(24)
as:
NPT S =
izl
R,2
=
minimum (NPT = number of observations).
(25)
The general problem of least-squares adjustment of data has been presented thoroughly by Deming and the reader is referred to Ref. 5 for details of the treatment. The function R in equation (24), in Deming’s notation corresponds to the condition equation: F(x, y, a, b, c) = 0 where x, y are the independent and dependent variables, respectively, and a, b, c are the parameters the least-squares estimate of which is needed. The least-squares solution of (25) is:
2
RF. R:1 Ak, + 2
RP;R:l+...+Z:
I
i
Li
R:l L,’ Ri” Ak, = 2\ Rt1.L R, i
I: RikyiR:L Ak, + 1 RikyiRiks Ak, + . . . + 1 Rpi,Rp
Ak, = 2 y
2 Rik2,RF’ Ak, + c ‘Pi$:’
Ak, = z y
Ak, _/_. . . + C Rp~iRik’
(26) (27)
(28)
where Rikj = 6Ri/6k,, (j = 1,2,. . . , n, and n = number of rate constants to be refined or estimated), R, = function R at ith observation, Li = (~R,/c%$)~. All NPT
summations refer to 2. The partial derivatives GRJSk,, etc, can be obtained by i=l numerical differentiation using the method of symmetric finite difference. The equations (26)-(28) are known as normal equations, which can be solved for Ak,, Ak,, . . . , Ak,,, etc. There are numerous ways in which these equations can be solved. In the programme CORNEK, the matrix inversion subroutine (MTXIN) has been used for this purpose. The new improved set of the rate constants is thus obtained: krf = kj” - Ak,,
(j = 1,2, . . . , n).
(29)
The superscript here refers to iteration number. The procedure is repeated until the function S converges to a steady value. The iterative adjustments are necessary except for the simple, linear case where the first iteration yields the unique solutions. The simplest iterative procedure 6*7*8 to apply is to take the final estimates in equation (29) as initial estimates and repeat the calculations. THE PROGRAMME CORNEK CONSTANTS FROM NEAR
(CALCULATIONS OF RATE EQUILIBRIUM KINETICS)
The programme is in five parts, the main programme and four subroutines. The main programme deals with the input data, the calculations of derivatives, aR,/ak,, etc, by numerical differentiation, the setting up of normal equations and
V. S. SHARMAand D. L. LEUSSING
1142
output of results. The statistical parameters such as variance of the function R, standard deviation in each of the constants estimated, covariance matrix, and correlation matrix are also calculated in the main programme. The subroutine MTORIG simply defines the elements of the coeflicient matrix of equation (20). The subroutine MTXOP performs three matrix operations : 1. It operates on the coefficient matrix of equation (20) and reduces it to the coe~cient matrix of equation (21). 2. It obtains the values of the determinant of the coefficient matrix of equation (22) by reducing the matrix to a lower, triangular form. The determinant is then given by the product of the elements of the main diagonal:
3. In the subroutine EIGEN it is again called for obtaining the numerical value of the polyno~al given by equation (23). The subroutine MTXIN is a standard matrix inversion subroutine and is used to solve the normal equations. After the final refinement of rate constants, the subroutine EIGEN calculates the predicted values of the relaxation times for comparison with the experimentally obtained values. The programme EIGEN is also completely general and is based on Muller’ss method of quadratic interpolation. The programme has provision for applying relative (% of Iz observed) weights to the dependent variable, which is L in this case. In the absence of the knowledge of weight, unit weights are applied to all observations. The programme can also be used for calculating only the roots (or relaxation times) without refining the constants if so desired. When the shifts in the constants have been calculated, the computer checks their values and reduces any which are excessively large. At no time is a shift allowed to be greater than the corresponding rate constant in that iteration. For the purpose of numerical differentiation to obtain the derivatives in equations (26-28), a small increment (4-l %) is given to the jth parameter, the other parameters remaining unchanged, and the value of aR,/ak, is calculated for all points. The same increment (O-1%) is then given to each of the other parameters in turn and the remaining derivatives are thus obtained. For each reaction scheme a separate MTORIG subroutine is needed. However, if only new reaction paths are to be considered without changing the number of rate equations, no change is needed in the subroutine, provided in the original derivation of the rate equations these paths have been taken into account. It, therefore, saves considerable time and computations, if all possible reaction paths are taken into account in the derivation of rate equations. The computer can then be asked to ignore any one or more of these paths by simply writing zeros for the corresponding rate constants in the input cards. The programme calculates variance in the function R in two different ways:
Var(R) = zy.
Computer approach for calculating rate constants
1143
Equation (30) is valid only at the minimum, whereas equation (30’) is a general one. Near the minimum the two equations are equivalent. In a good least-squares fit, variance in the function [equation (25)] being minimized continuously decreases and converges to a constant value, shifts in the constants [Ak, in equations (26-291, become progressively smaller, the gradient of the residual surface with respect to the parameters being estimated becomes very small, and finally on back-calculation, the calculated values of the dependent variable should be well within the acceptable limits of the technique. However, divergence or inability to converge is sometimes observed. This behaviour is quite often due to an inappropriate reaction mechanism chosen as the model, insufficient accuracy in the calculations of the derivatives, bad initial approximations of the parameters to be refined, and in some cases, depending on the nature of the “residual surface,” inappropriate size of shifts in the constants. The programme CORNEK contains several print statements converted into comment statements (see listing of the programme). These can readily be used to ascertain that the differences in the quantities [E(l), . . . , E(5)], which calculate the derivatives, are significant but as small as possible. The divergence due to bad initial approximations and inappropriate “step-size” can be taken care of by increasing the number of refinement cycles. In most cases, if the model is a correct one, convergence to a stationary minimum is attained in spite of initial wobbling of the function being minimized. However, if the initial approximations of the parameters are too far away from the values at the minimum, convergence in the error function is unlikely, no matter what size of shifts is chosen. The programme has provision for applying fractional shifts if it is felt that in a particular run divergence or wobbling in variance is due to the shifts being too big. At the end of each cycle the computer checks the variance and if it is higher than in the previous cycle, the shifts are reduced in the steps of SHIFT x [Fractional Shift Factor]. This may be done fifty times in a cycle and if, even then, no step-size gives lower variance the appropriate message is printed out and the process of refinement is terminated. The calculated values of relaxation times are then obtained by using the constants estimated in the last refinement cycle. If reduction in step-size improves variance, the programme continues the refinement with reduced step-size in the usual manner. A good model with reasonable initial guesses rarely requires this treatment. In some runs it may be observed that the value of variance at the final stationary point of convergence is slightly higher than in some earlier cycle. This usually arises from instability in the process of numerical differentiation. However, the corresponding error in the parameters will rarely exceed the estimated standard deviation in the parameters, provided all other requirements of a good fit are satisfied. Our experience with the programme is that in most cases all the conditions of a good fit are met only when the variances in R calculated by both equations (30) and (30’) agree closely (i.e. to the first or second decimal place). CONCLUSIONS
The least-squares refinement approach provides a powerful and useful tool for calculating the rate constants from relaxation times. In the CORNEK-II version of the present programme, provision is being made to refine the stability constants also. Deming’s general treatment does not require that the calculated value of the dependent variable, relaxation time in the present case, should be obtained in order to make the “best fit.” This is very helpful in relaxation studies where experimentally only one or
1144
v. s. sHAaM.4and D. L. Lmussr~o
two relaxation times are observed, whereas the calculations give a spectrum of relaxation times. The programme has been tested to calculate the rate constants for binary systems such as Cu2+-se&e, Cu2+-histamine, as well as for ternary systems involving five relaxations, such as Cu-histamine-serine. In all cases reasonable values of rate constants were obtained and good agreement between the observed and predicted relaxation times was achieved for all observations. At present the limiting factor in the near-equilibrium kinetic investigations is the accuracy with which relaxation times are measured. It will increase the reliability of the calculated values of rate constants considerably if ten relaxation spectra are taken for each observation so that the variance can be calculated and hence the weight of each observation. The values of the rate constants calculated also depend on the values of formation constants chosen for calculating the equilibrium concentrations of various species. The practice of using the formation constants determined at some ionic strength and temperature other than that maintained in the relaxation studies can also lead to meaningless results. In order to get reliable results from the programme it is essential that the relaxation times cover at least a tenfold variation in the concentrations of reactants and as many observations as possible. AcknowJec&ment-The authors wish to express their appreciation to the National Science Foundation for support of this work. APPENDIX
I
Ternary system Cu(II)-histamine-scrine M = Copper
A = Neutral histamine B = Serine anion The enthalpy change for the formation of copper chelates under consideration is appreciable (-20 kl/mole per M-N bond), and therefore, for a 10” temperaturejump relaxations could easily be obtained without using an indicator. The equilibrium expressions characterlzin~ the distribution of proton and metal complex species are: H+A+HA (31) H+I-IA+H,A H+B+HB H +HB trH,B
(32) (33) (34)
A-l-M+MA
(35)
A+MA+MA,
(36)
M+B+MB
(37)
B+MB+MB,
(38)
M+A+B+MAB.
(39)
In addition to equations (31-39), the following likely reaction paths were also considered: M+HA+MA+H
(40)
MA+HA+MA,fH
(40
MA+B+MAB
(42)
MB+A+MAB MB,+A+MAB+B MB,+HA+MAB+B+H MB+HA+MAB+H.
(43) (44) (45) (46)
Computer
approach
for calculating rate constants
1145
The mass balance equations for the system are: TH =
[Ml + [MAI + M&l + [MB1-I-MB,1 + [MAB]
(47)
TA
=
[Al + [MAI + 2Mhl+
(4)
TB
=
PI + [MB1 + 2WB,l + [H,Bl f [HB] + [MABI
TH = 2[H,Al+
[Hal
[HAI + 2[H,Bl+
+ [HA] + [MAB]
[HB] + [HI - [OH].
(49) (50)
The differential form of the mass balance equations is: 67-x = 0 = d[Ml + S[MA] + S[MAd + d[MB*] + d[MAB] 6TA=O=
(51)
IA1 + ‘H;;; + IHA1)d[Al + d[MA] + 2d[MA,l + d[MAB]
+ dTB
=
o =
(PI + ‘H;; + IHBI) WI
+
d[MBl+ ti[MB,]
2cH,A;H; rHA1)dM
+ d[MAB]
2’HsB;H; IHB1)d[HI
(53)
[HAI -I- 4[H,Bl f [HI + K,IIH 1 dWl. [HI
(54)
+ dT,=O=
‘WA;;
IHAl)5iAl + (Z1H.B;;
(52)
[HBl)dpl
(
+
4[H.A1+ (
From equations (3146) one can obtain five near-equilibrium rate expressions for dd[MA]/dr, dd[MAd/dt, da[MB]/dr, dd[MBd/dt, and da[MAB]/dt. These, with four mass balance equations [differential forms (51)-(54)], give nine equations, which in matrix notation acquire the form of equation (20), i.e.: z11.
....
.a,,
......
...... zs1. .....
aID
dIMA
dd[MA]/dr
dIMA,,
dd[MAA/dt
d[MBl
dd[MB]/dr
d[MBzl
dd[MB&dr
d[MAB]
dd[ MAB]/dr
d[Ml
o-0
d[Al
0.0
d[Bl
0.0
d[Hl
0.0
APPENDIX
II
Input data andformats 1. NJ (Number of jobs. For trying more than one [e.g. n] reaction model, put n.) FORMAT (12) 2. FSF (If fra&&al shifts are desired, put the value of “fractional shift factor,” e.g. for halfshifts write 0.5; otherwise, leave blanks.) FORMAT (F10.4) 3. NPT (number of observations); NC (number of constants required to define all the equilibrium expressions, e.g. 16 for the ternary system described in Appendix I); NCV (number of constants being refined or estimated); NL (number of types of ligands present in the system, e.g. two for a ternary
1146
V. S. SHARMAand D. L. LEUSSING
system and one for a binary system); TNEQ (total number of equations necessary to define the subroutine MTORIG; i.e. the number of mass balance equations plus the number of rate equations, e.g. nine for the ternary system); NREQ (number of rate equations, equal to the number of complex species); NCD (number of cycles of refinement required, usually 5-10); NDPA (number of dissociable protons on ligand A, e.g. two for serine, three for histidine and one for acetic acid); NDPB (number of dissociable protons on l&and B); WTF (weighting factor, unity in the absence of knowledge of weights of the relaxation times; if, on the basis of experimental data, it is concluded that in a particular set of experiments the relaxation times are accurate within &x%, then the weighting factor is approximately equal to x/200). FORMAT (2 x ,9 (12,4 x), F10.5) 4. PH @H); M (free metal concentration); A (concentration of free ligand A); B (concentration of free ligand B); LAMBDA (inverse of relaxation time). Number of such cards = NPT FORMAT (F10.4, 3E10.5, F10.4, E10.5) 5. SC (logarithm of stability constant); RC (rate constants for the forward reaction. In the case of rate constants accurately known, put the exact value. In the case of reaction paths to be ignored put 0.0 and put approximate values of the constants to be refined.); EQEXP (equilibrium expression corresponding to the equilibrium constant on the card, e.g. M + A = MA. This information is only for the purpose of print-out and does not enter into calculations.) Number of such cards = NC FORMAT (F10.4, D10.4, 5 x, lOA4) 6. CV (equation number of the constant to be varied; for example, to refine the forward rate constant for the reaction M + A = MA, the number needed in this card is 35 on the basis of the reaction scheme given in Appendix I); ICV (increment in the corresponding constant for the purpose of numerical differentiation; this number should usually be less than 1% of the expected value of the corresponding constant.) Number of such cards = NCV FORMAT (13, 7 x, E10.4) 7. SCHEME (This card is only for the purpose of record. On these cards can be written the left-hand side of rate equations, which in a way records the types of relaxations dealt with, e.g., DDMA/DT, DDMA*/DT, etc.) Number of such cards = NREQ FORMAT (2OA4) 8. If NJ is greater than one, repeat steps 2-7 NJ times. Preparation of subroutine ‘MTORIG’ The manner in which the present programme is written makes the consistency of notation very imnortant. However. there should be no difficultv if the followine nrocedure is strictlv adhered to. ‘1. Write down the stepwise equilibrium expre6sions characterrAg the formation df proton and metal complex species [e.g. equations (l-7)]. Give each equation a number (1,2, . . .). These equations should define all the known and postulated reaction paths. The corresponding equilibrium constant and the forward rate constant will have this number as the subscript. The tirst n equations (n = NDPA) should be for the proton-ligand association equilibria, starting with the lowest protonligand complex (e.g. H + A + HA). If the system has two types of ligands, then the next m equations (m - NDPB) should be for the proton association equilibria of the second ligand. The first and second ligands are always referred to as A and B, respectively. 2. On the basis of the reactions written in step 1, derive N (number of metal complex species) independent near-equilibrium rate equations [e.g. equations (14) and (15)]. 3. Write mass balance equations and obtain their differential forms [e.g. equations (8-11) and (16)-(19)]. This gives all the equations needed for the subroutine MTORIG, which consists of equations that define the elements of matrix A [equation (20)]. The programme requires that the first n (n = NREQ) column vectors [matrix A, equation (20)] should be the coefficients of differentials of metal complex species (6[MA], 6[MA,], etc). In the listing of the programme the subroutine MTORIG gives equations for the copper(H)-histamine system. This subroutine will apply to any system containing a metal ion with co-ordination number of four and a ligand with two dissociable protons. It is assumed that the only metal complex species are MA and MA*. These equations can be checked against equations (14-18). Since no indicator was used, equation (19) and the terms [HIn], ([HIn]/[In])GDn] from equation (18) have been deleted. Variables used in the main programme Variables used for input data have been explained earlier and are not listed here. Reserves the value of RC(I) at the start of each cycle. RCORIGQ
Computer LMORIG(I) DLAM(I) H(I)
ii3 E(2), E(3) E(4). E(5) WT(I) DRRC(I,J) DRLAMQ ROU) SSQR SNUTS SNUTl(NCC) Z(C) RCBEK(Cl Bl(I,J) . ’
Y(I)
S,SSQRo =-DO ROOT(N)
approach
for calculating rate constants
1147
Reserves the value of LAMBDA(I). Increment in LAMBDA(I) for the purpose of numerical differentiation. Hydrogen ion concentration. Elements of the coefficient matrix corresponding to equation (20). Value of the function R (equation 24) with current values of the parameters. Index used to identifv the constants to be refined. Counter for E(I). ’ Two values of the function R for the purpose of numerical differentiation of the function with resnect to constants IRC(C)l to be estimated. Two values of the function R for the l&u-pose of numerical differentiation with respect to LAMBDA(I). Weight of the ith observation (= 1 in the absence of knowledge of weights). Derivative of the function R at ith observation with respect tojth constant. Derivative of the function R at ith observation with respect to dependent variable [LAMBDA(I)]. Reserves the value of E(1). Sum-square of the function R over all the observations. Variance in function R calculated from equation (30). Reserves the value of SNUTS in each cycle for comparison with the last cycle. Number of iterations while applying fraction shifts. Shift in variable constant C, obtained by solving normal equations. Reserves the value of Cth constant. Left-hand side elements of normal equations. Later on, Bl(I,J) represents the elements of the inverse of the coefficient matrix (variance-covariance matrix) of normal equations. Still later (statement 177) it represents the elements of the correlation matrix. Right-hand side elements of normal equations. Variance in function R calculated from equation (29). Standard deviation in ith constant. Calculated values of LAMBDA(I).
1148
V. S. SHARMAand D. L. LEUESIN~
C
PROGRAMME
000 1
CORNEK
DOUBLE PRECISION A1~20,20)rB1l2012011D(20~20~~Yl20 1RCf2O~,RCDRIG12OlrDRRC(20,2O~,RO~2O~~SSORO~S~BB~2O
l,DRC(lO),E r20)*WT(20)
RODT~10),A11~20,20lrDRLAM~2O)tLMORIG(ZO) lSTOl10~. 1L~201,R1,R2,R3,RM,SNUT1120)1RCBEK~20) INTEGER CV,C,TNEQ DIMENSION H~20),PH(20),A(20),BotH2A(ZO),HA120),H2B~2O~~ 1HB~2Ol,LA~BDAl2D~,SCl2D~,CV~lO~,ICV(10~,SCLG~2O~~ 1SCHEMEl20,201,REACTl20,20~,EQEXP~20~201~0LAM~20~~
,SIGMAllO),
0002 0003 0004 0005
(5). t
0006 0007 0008 0009 0010 0011
lH3AT20),H4AT20),H3B(201,H48(20) REAL Hl20~,1CV,lAMBOA,1Nl2O~ COMMON M,A,B,H,RC,SC,HA,H2A,HB,H2B,Al,IN, lH3A,H4A,H3B,H4B,I 1 FORMATlIZ) 101 FDRMAT(2X,9112,4X),FlO.5~
0012 0013 0014 0015 0016 0017
102 103 104
0018 0019 0020
FORMATlFlO.4,3E10.5,FlD.4,ElO.5l FORMATI F10.4,010.4,5X,10A41 FORMAT(2OA41
0021 0022 0023 0024 0025 0026
105 FORMATI 13,7X,E10.4) 106 FflRMAT~13,2X,~10.4,5X,020.11,5X,10A4,5X,010.41 107 FORCATIlX,ZOA4) 108 FDRHAT(2X,4D2D.lD) 109 110
FORMATT lH1.4DlD.4) FORMATl1H1,10010.4,/,lOOlO.4) FORMAT(///lX,’ VARIANCE
111
EXT
UTS:
‘020.10,5X,’
VARIANCE
EXT
NUTS:
1020.10,5X,'NIT =' 13) 112 F0RMATl2X,13,3X~F6.0,5X,5(~17.10,2X~) 113 FORMAT~2X,E10.4,2X,FlO.4,5X,OlO.4~5X,lOA4,5X,I2~ 114
FORMAT(////lX,’ FORMATl70X.Dl0.4) FORMATlb(2X.ElO.4)) FORMAT (1Hlr FORMATI lH1,
115 116
121 124 125
FORWATllHl,‘NPT WTF’) 1’ FORMATllHO,‘THE
126 127
’
VARIANCE-COVARIANCE CDRRELATION MATRIX: NC RELAXATIONS
NCV
NL
128
1 129 131 132 133
B
MATRIX: ‘/) TNEQ
CONSIOE’REO
FORMATl2X,‘EQIL.‘,llX,‘LOG~SC)‘,4X,’RATE 1’/2X,‘CONSTlSCl FORMATllX,’ IRC)
CRC) STO.OEV.
‘/I NREQ ARE
NC0
NOPA
0l-j 66 I=l,NREQ REAOt5,104)
CV(I7,ICVtI)
(SCHEHE(I~IC)~
It=11201
NDPB’
FOLLOWING:‘/) CONSTANT’,2X,‘EDUILIBRIUM EXPRESS ION’//) EQUILIBRIA
SHIFTS’ 1 FORMAT{ 1X. ’ OBSD LAMBDA CALCULATED VALUES OF LAMBDA’ FORMAT(5(2X,D20.10)) FORMATl/2X,‘NO STEP SIZE GAVE LOWER VARIANCE THAN THE LAST FORMAT(F10.41 READt5,l) NJ NJO-0 READ(5,133) FSF READl5,lOl) NPT,NC~NCV~NL,TNEQvNREQ1NCD1NDPAvNOPB~WTF WRITEl6,1251 HRITEl6,lOl) NPT,NC,NCV,NL,TNEO,NREQ,NCD,NDPA,NDPB,WTF DC i I=l,NPT R~~o15,1021 PH~I~~M~I~~A~I~~B~I~1LAMRDA~II1INo 00 5 111, NC REAOf5,1031 SC~I)~RC~I)~(EQEXP~I~IC)~IC~~IC=~~~O~ 00 6 I=l,NCV READ(5.105) NRITEt6,126)
0027 0028 0029 0030 0031
‘I
’
’
/I CYCLE’)
0032 0033 0034 0035 0036 0037 0038 0039 0040 0041 0042 0043 0044 0045 0046 0047 0048 0049 0050 0051 0052 0053 0054 0055 0056 0057 0058 0059 0060 0061
1149
Computer approach for calculating rate constants
66
WRITE16,107)
0062
~SCHEME~l~ICl~IC=1,20~
0063
WRITEt6.114) NREQPl=NRED+l DO 12 I=l,NC SCLGII)=SCtI) 12
7
sctIT=lo.**sclI~ RCORIGII)=RClIT WRITEf6,127T DO 7 I=l,NC WRITEl6rll3~ SC~I~,SCLG~I~,RC~Il,(EQEXPfI,IC~,IC=l,lO~,I WRITE(6.114) DO 13 I=lrNPT LNORIGIIl=LAMRDAtI) DLAMII)=00.0001*LAMBDAlIT CALCULATES
CONCENTRATIONS
HII)=lO.**(-PHIIT) HA(I)=SC(l)*H(I)*A(I) IF(NOPA.GE.2) IF(NDPA.GE.3) IFlNDPA.GE.4) IFINL.LT.2) Ll=NDPA+l HBlI)=SCTL1)*H~Il*B~Il LZ=NOPA+Z
13
31 130
1000
15
c C
49 150
t
16
s 177
PROTON
COMPLEX
SPECIES
0079 0080 0081
H2A~Il~SCl2l*H~I~*HA~f~ H3AfIl=SC(3l*H~I~*H2AII) H4AlIT=SC~4)*HlIT*H3A(I) GO TO 13
IF(NOPG.GE.2) L3=NDPA+3 IFlNOPE,GE.3l L4=NDPA+4 IF(NOPB.GE.4)
48
OF
0064 0065 0066 0067 0068 0069 0070 DO71 0072 0073 0074 0075 0076 0077 0078
0082 0083 0084 0085 0086 0087
H2EfI)=SC(L2l*HI
I)*HB(I
I
H3B(IJ=SCfL3T+Hf
I )*HtB(
I)
0088 0089 0090
CONTINUE NCC=O.O CONTINUE NCCMl=NCC-1 NIT=0 CONTINUE
0091 0092 0093 0094 0095 0096 0097 0098 0099
DO 19 DO 19 ID=1
0100 0101 0102
H4RII)=SCfL4)*HlI)*H381I)
J=lsNCV I=lrNPT
CONTINUE CALL MTORIG CALL MTXOPlAlrTNECI,YREOPl) DO 15 II=l,NREQ Al(IIrII)~AlfIIrII)-LAHBOAO IFiNREO.EO.l)
GO
0103 0104 0105
TO
0106 OlOi 0108
150
CALL HTXOP(Al*NRE0,2) DO 49 IV=lrNREQ HRITE(6v131) (Al~IV~IVST,IVS~lrNREP) E(ID)*l. DO 16 N=l,NREO EtID)=ElIDT*Al(N,N) BEGINS
NUMERICAL
0109 0110 0111
DIFFERENTIATION
fiOtI)=E(l) 15=10+1 GO TO ~177,177~17,18,61~119)( C=CVI Jl RC(Cl=RCtC)+ICV(JT
ID
0112 0113 0114 0115 0116 0117 0118 0119 0120 0121 0122
V. S. SHARMAand D. L. LEUSSING
1150
GO 17 18
C
TO
1000
0123
RCTCl=RClC)-Z.*ICVLJ) GO TO 1000 DRRClI,J)=~E~2l-El3)l/~2.~ICV(J)~ WRITE(6.1311
40
ETl),E(2)1
RClCl=RCORIGICI IFLJ.LT.NCV) GO TO LAHRDALI)=LACRDATI)+DLAMO
0124 0125 0126
E(3)
0127 0128 0129 0130
19
GO
61 119 c 19
TO 1000 LAMROAII)=LAHRDAlI1-Z.,DLPMII) GO TO 1000 DRLAM(I)=(E(4)-El5))/I2.*DLAM(I)) WRITE16,131) Ef4),ET5) LAHBDA(I)=LHDRIGII)
0131 0132 0133 0134 0135 0136 0137 0138 0139
c SSQRO=O.O SSQR=O.O DO 27 I=l,NPT WT(I~=(l./(LAHRDATIT IF(WTF.EQ.1.0) 27
190 C C
62
0140 0141 0142
*WTFII+*2 WTLI)=l.D
L(IL=DRLAMLI)**2/WTlIl DO 190 I=l,NPT SSQR=SSQR+RDTI)*RO(I)
0143 0144 ii45
SSQRO=SSQRO+RDlI,**2/LO DO 62 I=l.NPT WRITET6.1311 ~DRRC~I,Jl,J=l,NCV~,L(I),RD~I~,SSQR SNUTS=SSQRO/ INPT-NCV) SNUTl (#CC I-SNUTS
C APPLIES FRACTIONAL IS GREATER THAN
c
THE
SHIFTS LAST
IF THE CYCLE
VARIANCE
E
662
33
34
IF(NCC.EQ.0) GO TO IFTFSF) 662134,662 1FISNUT1~NCCI.LT.SNlJT1~NCCH1l~ NIT=NIT+l DO 33 I=l,NCV C=CVL I) DRCTI)=DRC(I)*FSF RCLC~=RCORIGlCl+DRClI) RCDRIGLCi-RCICT
34 GO
IFLNIT.GT.50) CONTINUE IFlNIT.GT.50)
RCLCl.RCBEKlC)
IFlNIT.GT.50) GO TO 31 CONTINUE
GD
WRITET6~1321 TO
98
c SETS
UP
NORMAL
EQUATIONS
:
20
21 C 51
DO 20 I=l,NCV DO 20 J=l,NCV Rl~I,J)=O.O 00 20 KK=l,NPT B1~I,J)=B1~I~J)+DRRC~KK11)*DRRC~KK~J~/LIKK~ I=1 ,NCV DO 2: Y(I)=0 DO 21 KK=l,NPT YTI)=YTIl+RDLKK~*DRRC(KK11)/L(KK) no 51 I=l,NCV MRITET6~116) ~Bl~I,Jl~J=l~N~V~~YII~
TO
34
IN THE
CURRENT
CYCLE
0146 0141 0148 0149 0150 0151 0152 0153 0154 0155 0156 0157 0158 0159 0160 0161 0162 0163 0164 0165 0166 0167 0168 0169 0110 0171 0172 0173 0174 0175 0176 0171 0178 Dl79 0180 0181 0182 0183
Computer approach for cabdating
CALL
c c
DO
22
1151
0184
MTXINlBl,NCVL
CALCULATES
c
rate constants
22
0185 0186
SHIFTS
0187 0188
I=l,NCV
0189 0190 0191
DRC( I)=O.rJ DO 22 J=l, NCV DRC~I)=DRC~l)+RllI~J~oy(J~
0192 0193 0194
C
122
123
23
120 c C C
S=SSCJRO DO 122 I=l,NCV S=S-YIlI+ORClI) S=S/lNPT-NCVT HRITEl6,lll~ 00 123 I=l.NCV SIGMAll~=R1~I,I)*SNUTS STD(I)=SIGMAlI)*+.5 WRlTE~6,1281 DO 23 N=l,NCV C=CVLN) WRITE(6.106) CONTINUE lF(NCC.EO.NCD)
0195 0196 0197
SrSNUTS,NIT
0198 0199 0200 0201 0202 0203 0204 0205 0206
C,RClC),STD(Nl,iEDEXPO,NIl,NI=l,lO~,DRCLN~ GD
TO
98
0207 020R 0209
DO 120 I=l,NCV DO 120 J=l.NCV BBLl,J)=B?(I,J) APPLIES GREATER
‘UNRESTRICTED THAN THE
SHIFTS’ CORRESPONDING
AND ND SHIFT RATE CDNSTANT
IS
ALLOWED
TO
0210 0211 0212
RE
0213
C DO 46 I=l,NCV c=cvI I) RCLC)=RCORIGlC) RCLCl=RC(C)-DRCLI) IF(RCiCT.LT.O.0) 1RCDRIGlC) RCBEKICT=RCDRIGLC) RCORIG(C)=RCLC) 46 CONTINUE 47
0214 0215 0216 0217 0218
RCIC)=RCORIGLCL-DRClI)/(DARS(ORCrI)))+.S*
0219 0220 0221 0222 0223 0224
c NCC=NCC+l GO TI-I 48 C C C C
AFTER FINAL CYCLE OF CDHPLJTED AND PRINTED 90
CONTINUE WRITE~6,1141 WRITEl6,1291 DD 29 I=l,NPT CALL CALL CALL
29
REFINEMENT
ItLAMBDA~I~,LROOTlN~,N=l,NREB)
C C
PRINTS
VARIANCE
OF
LAMBDA
ARE
0225 0226 0227 0228 0229 0230 0231
60
TO
2222
CflVARIANCE
0234 0235 0236 0237 023R 0239 0240 0241 0242
c IF(NCD.EU.fl.0) WRITEl6.114) WRITEI 6,121)
VALUES
0232 0233
MTORIG MTXOPLAl,TNEO,YREOPl) ElGENlAl,NRFD,RDOT)
WRITEL6.112) CDNTINIJE
CALCULATED
MATRIX
0243 0244
11,52
V. S. SHARMAand D. L. LEIJSSING
0245
C 9
c c c 117
118 2222
PRINTS
CORRELATIDN
DO DO
I=l,NCV J=l,NCV
0246 0247 0248 0249 0250 0251 0252
BlLI,Jl=B8LI~Jl/ WRITEL6,1141 WRlTEL6,124) DO 118 I=l,NCV HRITEl6.1161 CONT INVE
0253 0254 0255 0256 0257 0258
00 9 I=lrNCV wt11~I6,1161
117 117
(Rl(IvJ)eJ=ltNCV) MATRIX
LBlII,JIvJ=l,NCVl
NJD=NJD+l IFINJD.LT.NJl STOP
GO
TO
8
0262 0263 0264
END ct SUBROUTINE MTORIG DOUBLE PRECISION REAL M120I,INl2DI
0259 0260 0261
AlL20,2OI,RC(2OI
DXMENSION A1201,8i20I,H120l,SC01HB01H2B01H3BL20 1HAi20l.H2Al2OI.H3A~2Ol~H4AL2OI
0265 0266 0267 0268 0269 0270 0271 0272 0273 0274 0275 0276 0277 0278 0279 0280 0281 0282
AlL3,3)=1.0 AlI3.4)=0.0 Aii3;5)=0.0 AlL4,1l=l.O A1(4,21=2.0 A1(413I=O.O Al~4~4~~~H2A~I~+MAlI)+A~I)I/A~I~ A1L4,5l=l2.*H2A~Il+HAo)/H(I) AiL5,1)=0.0
c c
0283 0284 0285 0286 0287 0288 0289 0290 0291
0292 0293
0294 0295 0296 0297
AlL5,21=0.0 A1(5,31=0.0 A1(5,4)=(2.*H2A(II+HA(I))/A(I)
0298 0299 0300
F1(5,5)=(4.9HZA(I)+HA(I)+H(I))/H(I) RETURN END
0301 0302 0303 0304 0305
Computer approach for calculating rate constants
SUBROUTINE DOUBLE
MTXOPtOvIP PRECISION
00 13 N=l,IP IR=IP-N+l IFIIR.LT.IR)
RRIO
GO
TO
,101 ~20,201
1153
0306 0307 0308 0309 0310
13
0311 0312 0313
M=IR-1 00 13
I=l,M RR=DiI.IR)/D~IR~IRl IFIRR)i0,13,10 10 00 12 J=l,IR 12
OII,J)=D1I,JkRR*D~IR,J)
0314 0315 0316
13
CONTINUE RETURN END
0317 0318 0319
SUBROUTINE
0320 0321 0322
c C DOUBLE DO 764
HTXINlA760,N760) PRECISION K760=1.N760
0323 0324 0325
6760120,201
765
DO 765 J760=l,N760 A760fK760rJ760)=A760IK7~O,J76O~/CM760
0326 0327 0328
762
DO 764 1760=lrN760 IF( 1760-K760)762,764v762 CM760=A760( 1760,K760)
0329 0330 0331
763
A760lI760,K760)=0.0 DO 763 J760=1,N760 A760~1760~J7601~A760(I760~J7601~CM760*A760~K76DiJ7601
0332 0333 0334
CONTINUE RETURN END
0335 0336 0337 0338 0339 0340
764
C C SUBROUTINE
111 112 113
EIGENIAZINR,RRTT
DOUBLE PRECISION A2l20~2DT~R~b)~RT~HH~G~D~LM~XN,TEMP,RRT~lOT~ lFPRT,FF(3),A3(20,20T FORHATl50X,E17.10~ FORHAT~2X,S(D17.10,2XTT FORMAT1 2X. 1 NONCONVERGENCE FOR THIS POINT’) DO 1 I=lqNR
1 RRT( 570 581
I)=O.O ERR=.01 DRT= .99 K=O K=K+l NIT=0 IF(RRT(KT.EO.O.0) RT=DRT*RRTfKl n=1
571
572
513
GO TO 580 RT=-1.0 l-B-l= -1.0 H=3 GO TO 580 FF( l)=FPRT RT=RRT IK )/DRT n=2 GC TO 580 Ff-(Z)=FPRT Hh-RRTIKI-RT
0341 0342 0343 0344 0345 0346 0347 0348 0349 0350 0351 0352
GO
TO
571
0353 0354 0355 0356 0357 0358 0359 0360 0361 0362 0363 0364 0365 0366
v. s. s-
1154
and D. L. -0
0367 0368 0369 0370 0371 0372 0373 0374 0375 0376 0377 0378 0379 0380
RT=RRTLK) H=5 GO TO 580 FFL l)=FPRT RT=l .O H=4 GO TO 580 FFLZ)=FPRT
5573
5574
RT=O.O H=5 CD TO FFl31=FPRT LH=-.5 GO TO
5575
5576
580
574
FFll)=FFIZ) FFL2)=FF(3) FF(3)=FPRT
574
DEL=l.+LM G=FF(1~*LM*~2-FFl2)*OEL**2tFFo*(LH+DEL~ D.G~+2-4.~FF~3l~LM~DEL~~FF~ll*LM-FFo+DELtFF~3~)
0381 0382 0383 0384 0385 0386
D=G+G/DABS(G)*DSDRTLDABSo) LH=-Z.*FF I3 )+DEL/D HH=LM+HH
0387 0388 0389
XN=RT+HH IFIDABSIHH/XN).LT.ERR) RT=XN
GO
XXX=DABS tHH/XN) NIT=NITtl lF(NIT.GT.100~ M=6 GO TO 580 575
580
2
3
4
600
590 591
RRTlK)=XN IFtK.LT.NR) GO TO 591 CONTINUE DO 2 I=l,NR DO 2 J=l,MR A3iI,J)=AZ(I+J)
GO
GO
TO
50
TO
0390 0391 0392
575
0393 0394 0395
591
0396 0397 0398 0399 0400 0401 0402 0403 0404
581
DO 3 I=l,NR A3lI,I~=A2lI,I)-RT CALL MTXOP(A3,NR,Z)
0405 0406 0407
FPRT-I .O DO 4 I=l,NR FPRT=FPRT*A3(1,11
0408 0409 0410
J=l
0411 0412 0413
IFLK.ECd) GO .TO 590 J=J+l TEMP=RT-RRT L J-I I FPRTsFPRT/TEMP 1FLJ.NE.K) GO TO 600 GO ‘TO (572,573,5573,5574$5575155761( CONTINUE RETlJRM END
0414
M
0415 0416 0417 0418 0419 0420
Computer Note added in proof: changes are made.
approach for calculating rate constants
For some computers line line line line
(card) (card) (card) (card)
the programme
0095; 0155; 0206; 0240;
change change change change
1155
will operate better if the following
0.0 to 1 0 to 1 EQ to GE 0.0 to 0.
Zusammeufassang--Ein allgemeines Maschinenprogram zur Ermittlung von Geschwindigkeitskonstanten aus Relaxationszeiten wird beschrieben. Das Programm CORNEK ist im wesentlichen ein auf nichtlineare Systeme angewandtes Verfahren der kleinsten Quadrate. Es verwendet direkt die differentiellen Formen der ersten Ableitungen der Geschwindigkeitsund der Massengleichgewichts-Gleichungen; damit werden die zcitraubenden Ableitungen von gleichgewichtsnahen Geschwindigkeitsgleichungen vermieden. Das Programm wurde an bin&en Systemen wie Cu-Histamin und Cu-Serin sowie dem terrdiren System Cu-Histamin-Serin tiberprtift. Rt%uru&Gn d&rit un ac& general par calculatrice a l’estimation des constantes de vitesse a partir de temps de relaxation. Le programme CORNEK est essentiellement un programme d’dtude plus fine par moindres car&s applique aux systemes non lineaires. 11 utilise directement les formes differentielles des premieres d&iv&es des equations d’dquilibre de masse et de vitesse, evitant ainsi les derivations des equations de vitesse au voisinage de l’equilibre qui absorbent du temps. Le programme a et6 essay6 pour des systemes binaires tels que Cu-histamine, Cu-serine, et le systeme ternaire Cu-histamine-&me. REFERENCES 1. M. Eigen and R. G. Wilkins, Mechanism of Inorganic Reactions. American Chemical Society, Washington, D.C., 1965. 2. 0. G. Hammes and J. I. Steinfield, J. Am. Chem. Sot., 1962, 84,4639. 3. V. S. Sharma and D. L. Leussing, Znorg. Chem., in press. 4. Zdem., Chem. Commun., 1970,1218. 5. W. E. Deming, Statistical Adjustment of Data. Dover Publications, New York, 1964. 6. H. 0. Hartley, Technometrics, 1961, 3,269. I. K. Levenberg, Quart. Appl. Math., 1944,2,164. 8. D. W. Marquardt, J. SIAM Sot., 1963, 11,431. 9. S. D. Conte, Elementary Numerical Analysis. McGraw-Hill, New York, 1965.