A general computer approach for calculating rate constants from near-equilibrium kinetic studies

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.