multiligand systems

Analytica Chimica Acta, 191 (1986) 385-398 Elsevier Science Publishers B.V., Amsterdam - Printed in The Netherlands

THE CALCULATION OF EQUILIBRIUM CONCENTRATIONS IN LARGE MULTIMETAL/MULTILIGAND SYSTEMS

ALESSANDRO DE ROBERTIS, SILVIO SAMMARTANO* Zstituto di Chimica Analitica

CONCETTA DE STEFANO

and

dell’Universitd,

Via dei Verdi, 98100 Messina (Italy)

dell’Unzversitii,

viale A. Doria, 95125

CARMELO RIGANO Dipartimento

di Matematzca

Catania (Italy)

(Received 26th May 1986)

SUMMARY An algorithm for computing equilibrium concentrations by the “equilibrium constant” method is described. The main features of this algorithm are: (a) a damping procedure in conjunction with the Newton-Raphson technique that avoids divergence in dealing with very complicated (simultaneous presence of simple, mixed, protonated, polynuclear and hydroxypolynuclear species) and/or very large systems; (2) the use of devices to decrease core requirements, calculation time, and ill-conditioned problems; and (3) the calculation of errors in free and species concentrations from the uncertainties in analytical concentrations and in formation constants. Four systems are used for testing computer programs on calculation of equilibrium concentrations.

The problem of finding species concentrations in a multicomponent system is of great importance in many fields. Geochemists and marine chemists are interested in studying the real composition of sea, estuarine, lake, river, interstitial and, in general, ground waters. Biochemists, bio-inorganic chemists and analytical biochemists are interested in the speciation of fluids such as blood plasma and urine. Analytical chemists are interested in general speciation problems and in implementing titrations, etc. Environmental chemists are interested in all natural-fluid compositions. Other chemists are often concerned with the composition of multiphase/multicomponent systems in kinetic, complex formation, synthesis and other studies. Many papers have been published on the problems of computing the concentrations of all species in a multicomponent system [l-34]. Two methods have been commonly adopted, i.e., the “equilibrium constant” and “free energy minimization” methods [7, 81. With regard to numerical approaches for the solution of equations defining the system under study, the algorithms may or may not involve matrix inversion. The commonest algorithm that involves matrix inversion is the Newton-Raphson one [35, 361. The reported algorithms and computer programs can deal with very different systems, such as simple acid/base, up to systems comprising several 0003-2670/86/$03.50

o 1987’Elsevier

Science Publishers B.V.

336

components and several phases. The possibility of simulating titration curves, of scanning analytical and/or free concentrations and of plotting distribution diagrams, is considered in some programs. The various computer programs have been designed for work with large, medium, mini and personal computers and have been written in different languages (generally, Fortran, Basic and Algol). Interesting applications of programs for computing concentrations in analytical, geological and biological problems have been reported [23,37-411. Topics of current interest in this laboratory are titration procedures for the determination of formation constants and formation enthalpies of weak, polynuclear and mixed complexes by potentiometry and calorimetry [4251] and the speciation of biological fluids [52], thus it was thought worthwhile to investigate further the problem of computing equilibrium concentrations, with the aim of examining in detail all related features. In this first contribution, a general algorithm is described for calculating concentrations in multicomponent systems. The reliability of the algorithm in handling four systems is tested. THE ALGORITHM

For a system containing N components and M species, the mass-balance equations are Ck = ck + 1 &k pi II 4” i

ck = ck +

j

1 Pik Xi

1

(1) (1’)

where ck and Ck are the analytical and free concentrations of the kth component, respectively, Xi is the concentration of the ith species, pi is the formation constant of the ith species, Pik (or Pij) is the stoichiometric coefficient of the kth component in the ith species, and the indices are defined as k = 1 . . . N, i = 1 . . . M and i = 1 . . . N. If Ck, pi and pi are known, Eqn. 1 represents a set of N nonlinear equations to be solved simultaneously in order to find N ck values. The problem is then to solve the set of equations fk(c)=ck+

FPikfli

~~-&=O (2) J subject to the constraint Ck > 0. The Newton-Raphson technique allows ck values to be calculated by the iterative procedure

Cc”+ 1) = @V_G-l

e

(3)

wherec= (c,...ck... CN),G= [dfk(c)/acj] =(&j)aIlde=(el...ek...e~) (ek = (&&d - Ck); at a certain stage n, the calculated analytical concentration is taken from Eqn. 1:

387 c(n)

k.calcd

=

Ck(n) +

[ EPik i

pi

yp](n)

(4)

This technique offers the important advantage that few iterations are generally needed to reach convergence, but two main difficulties: (a) when G is singular or near-singular, inversion is impossible; (b) the iteration is often very unstable and, at some stage, the correction G-’ e may lead to divergence. These difficulties must always be kept in mind when the Newton-Raphson technique is ,used, but are particularly severe in dealing with highly nonlinear equations and when N > 5. In order to overcome difficulty (a), scaling was applied to matrix G and to VeCtOr e according to the equations g& = gkj(gkk gjj)-1’2 and e$ = ekg$ where g:j and e;z*are the elements of the scaled matrix and vector, respectively. In order to improve the stability of the iterations, it proved best to use the following procedure. If in a certain stage n, the ratio @’ = Ck / for the kth component lies outside the range p-l < Rr) < p (p is a CklLal limit chosen in the range 1 < p < lo), then the free concentration of the hth component is damped by the equation ‘hrf’da,ped

=

(5)

cfin’(Rgp

where q = IPikl,& (i.e., q is the reciprocal of the largest stoichiometric coefficient of the species containing the kth component). This procedure is applied to the component for which lnl Rk I assumes the maximum value, and is repeated until Rk for all components satisfies the condition p-’ < Rf’ < p. Then a new Newton-Raphson iteration step is performed. For the solution of the system outlined in Eqn. 2, the compact Gauss method was chosen (modified Gauss elimination method, more easily programmable on computers) [ 531. Another important factor for the efficiency of the Newton-Raphson procedure is the choice of initial estimates, here ck(‘) . This choice is not critical when the above damping procedure is used but, for speed of calculations, it is advisable to adopt some devices. According to the available knowledge of the system under study, users may choose estimates as follows: c;) = lo-‘, Ck(‘) = ck or ck(‘) = 10-log(Ck), with -log (ck) supplied by the user. When several points for a titration curve or for a distribution diagram are to be calculated, the estimates of second and third points are the output of the first and second points, respectively. For the fourth and subsequent points, the empirical extrapolation equation (valid for a continuous monotonic function) is crP = ck,p -I + hk,p x [(xp

-2 -xp

- 1 - ck,p -2)/bp -d/@k,p

-2 -ck,p

-I

-

-dl

3cp-2)12 (“p -xp

-I)

(6)

where p (24) is the index of the point and 3tp is the titrant volume for simulating a titration curve or -log (CN) for calculating distribution diagrams. For constant increments, Eqn. 6 reduces to

388

$f_:,= ck,p-1

+

kk,p

-I

-d2/@k,p

-Ck,p

-2

-Ck,p

(7)

-a)]

Other devices for reducing calculation time and core requirements must be taken into account. Unfortunately, very often the reduction of core requirements leads to increasing calculation time. An interesting example regarding this problem is the storage of stoichiometric coefficients (a very sparse twodimensional array). Some reported algorithms [24, 311 reduce drastically the space required, but increase significantly the execution time. In the preliminary versions of ES4EC (see Results and Discussion), the calculation time was reduced as much as possible. However, probably the best way is to reduce storage when N > 20 and M > 100 but to reduce calculation time when N < 20 and M < 100. In practice, in all cases in which species concentrations are to be calculated, it is necessary to give not only the free concentrations of components and the concentration of species but also the errors in these quantities arising from the uncertainties in the analytical concentrations and in the formation constants. Thus, given the error propagation law, it is necessary to calculate 2 (JCj =

f(acj/api)2

2 ‘Xi

f

=

O~i +

(axi/i3/3,)2

~ (aCj/aCr)2

u$, +

F

“‘c,

(&xi/aC,)2

(8)

u’c,

(9)

where r = 1 . ..N and Z= 1.. .M. The partial derivatives in Eqns. 8 and 9 can be obtained by deriving Eqns. 1 and 1’ with respect to Pi and Cj: aCj/aLC = (acjlafli) acj/aG

= (acj/aG)

+

C (aMa&) r +

z

E (Pir

tack/a&)

(10)

Pik XiICr) + Pij Xi/Pi

F (Pir 1

Pik

xi/ck)

(11)

Because aCj/a& = 0 [54] and aCj/aC, = 6,, (&jr = 0 for j # r but = 1 for j = r; 6, Kronecker delta), Eqns. 10 and 11 represent (M + N) sets of N linear equations A x = b (M systems, one for each i) where A = (ajk) with

ajk

=

z(plj

(10’) plk

Xl/Ck)

-k

6jk,

X

=

(Kj)

With

Kj

=

!

acj/a&, and b = (bj) with bj = -pij Xi/Pi; and B y = d (N systems, one for each r) where B =

(bjk

)

with

bjk

=

1 (pij

i

pik

(11’) xi/Ck

) +

6jk,

Y =

(yj)

with Yj =

acj/aC,, and d = (&jr). Once the partial derivatives of free concentrations of components have been obtained with respect to pi and Cr, the analogous derivatives of species concentrations can be obtained from the equations axil@~

= (6iZ Xi/PI) +

4 bij

Xi/cj(acjlaP~)l

(12)

389

axiiac,= f

bij

Xi/Cj (acj/aC,)]

(13)

The values of aci/a& obtained from Eqns. 10 or lo’, aci/aC, from Eqns. 11 or ll’, axi/apl with Eqn. 12, and axi/aC, from Eqn. 13 allow u:, and “ii to be calculated from Eqns. 8 and 9, respectively. The procedure adopted does not take into account covariance terms in Eqns. 8 and 9; this is correct for dealing with data collected from various sources and obtained by means of different techniques, but is hardly to be recommended for simulating titration curves by using a homogeneous set of formation constants obtained in the same laboratory by the same technique [55,56]. RESULTS

AND DISCUSSION

A series of computer programs all dealing with the computation of equilibrium concentrations, ES4EC (equilibria in solution, problem 4, equilibrium concentrations), was written or will be written, in different languages and with different features. [In this paper, the computer programs ES4ECl and ES4EC2 are briefly discussed in order to evaluate the performance of the proposed algorithm; more details will be given in subsequent papers but copies of the ES4ECl and ES4EC2 listings (not definitive), are available on request.] A first edition, ES4EC1, was written in Fortran, on the basis of the described algorithm. The flow diagram of the program ES4ECl is reported in Fig. 1. This edition can deal with systems containing 20 components and 80 species of any kind in one homogeneous phase; it can simulate titration curves and calculate a series of points for distribution diagrams [the free concentration of the last component (-log CN) is scanned in a certain range] ; it can calculate errors in free concentrations of components and species. A second program, ES4EC2, was written in Basic, and shows the same characteristics as ES4ECl. The ES4ECl program was run in a Honeywell Level 6 (System 43) computer [program size, 17 kbytes; data size, 50 kbytes (double precision, 64-bit)] and the ES4EC2 program was run in an IBM AT personal computer (512 kb, double precision, 64-bit). In order to evaluate the reliability of the proposed algorithm and of the ES4ECl and ES4EC2 programs, four systems showing some interesting difficulties were selected. System I. This had two components (A13+ and H’) and 10 species: log 0 [AlH_,] = -4.4, log P[AlH_,] = -9.3, log /3[AlH_,] = -15.0, log /3[AlH_,] = -23.0, log p[A12H_J = -7.7; log /3[A13H_,J = -13.9, log P[AL,HJ = -28.2, log p [Al,H_,,] = -54.7, log /3[A113H_32] = -100.0 and log p [H_,] = log K, = -14.0. The conditions tested were CA1 = 3.7 X lo-‘, 10S6 and 10-l mol dme3 at several pH values. Formation constants of system I were taken from various sources. Though the considered constants do not represent the best set to explain the hydrolysis of aluminium, the proposed system shows an interesting variety of mono- and polynuclear species, so that various problems can be encoun-

Fig. 1. Flow diagram of the ES4ECl program.

tered if the algorithm used is not robust. In particular, the presence of species for which /3 < 10-50 can lead to underflow problems. In the ES4EC program, formation constants > 10ESPL or ESPL, the formation log pi = ilog p;NruT, and an index IBi = 2 is recorded (otherwise IS, = 1); (2) the species concentrations are calculated from the equation log X, = 1Bi log PI + C Pij log

Cj

391

System II. This had five components [K’, Cu’+, Ni’+, cit3- (citrate), H’] and 18 species: log p[H(cit)] = 5.99, log @[Hz(cit)] = 10.42, log P[H3(cit)] = 13.41, log P[K(cit)] = 0.74, log P[K(cit)H] = 6.15, log fl[K(cit)H,] = 10.26, log p[Cu,(cit),] = 15.28, log fl[Cuz(cit)zH_,] = 11.90, log /3[Cu,(cit)ZH_2] = 6.66, log p[Cu(cit)H] = 10.01, log P[Cu,(cit)H_,] = 5.37, log fl[Ni(cit)] = 5.56, log P[Ni(cit)H] = 9.59, log p[Ni(cit)l] = 8.85, log P[Ni(cit),H] = 14.34, log fl[Ni2(cit)2H_2] = -3.90, log p[CuNi(cit)H_,] = 2.37, and log /I[H_,] = log K, = -13.97. The conditions tested were Cx = 0.75, ecu = 0.003, C,, = 0.003, C,it = 0.008 mol dme3 at several pH values. System II, taken from earlier work [48] is an interesting example of studies of complicated equilibria with the formation of polynuclear and mixed species. System III. The fourteen components were Na+, K’, Ca*+, ac- (acetate), mal*- (malonate), male*- (maleate), succ*- (succinate), mala*- (malate), tart*- (tartrate), pht*- (phthalate), oda*- (oxydiacetate), gly- (glycinate), his- (histidine), H’ and there were 75 species (see Table 1). The conditions tested were CNa = Cx = 110, Cca = 10, c, = 1, c* = 0.1, CL = 0.5 mmol dmb3 (L = other ligands) at several pH values. This rather large system is an example of a multicomponent solution in which several weak species are formed. The formation constants for this system were taken from earlier work [44-47,49-511.

Fig. 2. Distribution of species vs. pH for the system IIa. Percentage refer to citrate. For the analytical conditions, see Table 3, and for the indices of species, see footnote (b) of Table 4.

392

TABLE 1 Log p values for the complexes formed in System IIIa Species

[HLI

Metal

-

W,Ll

-

[H,L] [ML]

Na K

[MLH]

g K

[MLH,]

;z

aLogp[NaH,]

Ligands ac

ma1

male

succ

4.571

5.36 8.06 _ 0.57 0.68 1.64 5.15 5.31 5.90 -

6.02 7.79 0.83 0.84 1.76 5.99 6.00 6.67 _

5.26 4.74 9.31 8.02 _ 0.47 0.3 0.48 0.38 1.45 1.95 5.26 4.7 5.28 4.5 5.96 5.72 _ _

-

_ -0.27 -6.47 0.57 1

=-13.94;logfl[CaH,]

mala

=-12.92;logP[H_,]

tart

pht

4.01 6.87 _ 0.58 0.5 2.10 4.05 4.0 5 02 _’

5.09 3.97 7.91 6.79 _ 0.73 0.34 0.83 0.25 1.71 3.46 4.96 3.50 5.16 3.50 5.70 6.38 =logK,

oda

gly

his

9.55 11.95 9.10 1.05 9.08 9.85 -

9.06 15.07 16.82 0.95 9.35 15.4

z-13.84.

System IV. This concerned the hydrolysis of six metals, MI--M6 and H’, i.e., 7 components and 35 species (see Table 2). The conditions tested were CM = 1 mmol dme3 (for all metals) at several pH values. This system is an artificial example, for which the hydrolysis of six metal ions, simultaneously present in a solution and forming several species with a high degree of polymerization and hydrolysis was simulated. This system showed the same difficulties as System I, increased by the presence of more than one metal ion. The ES4ECl program was tested with all the above systems, and no divergence difficulties were encountered. The ES4EC2 program was tested with Systems II and III, and again no divergence problems were observed. The execution time for both programs (and machines) was very reasonable for all systems, as reported in Table 3. Moreover, it must be considered that all systems were tested by also calculating errors in the concentrations, and this increases the calculation time considerably. [For the systems here proposed, the ratio r,/~ (7 = execution time; 7, = execution time when calculating errors) varies from 9 (System IIIa) to 2 (System IIb), according to the complexity of the system.] For clarity, only total iterations and functions evaluations, are reported in Table 3, but it should be noted that the first three points require generally a higher number of iterations (2-5 times) with respect to the mean value, and this indicates that the adopted extrapolation method (to find initial free concentrations, see Eqns. 6 and 7) allows a significant time-saving. For the damping procedure (see equations for Rp’ and c@’ k,damped), the values p = 4 (ES4ECl) and p = 2 (ES4EC2) were chosen for all systems. The choice of the factor p is not critical (p = 24 is a good choice in all these cases) but can have some influence on the number of function evaluations, in particular when initial estimates of free concentra-

393 TABLE 2 Equilibrium constants for the hydrolysis of metal ions in System IV* Species

W

M2

M3

-3.00 4.83 -8.29 -20.0 -32.2 -

-3.40 4.80 -7.10 -16.6 -6.88 -

-7.71 -17.12 -28.06 -23.88 -20.88 -

-16.17 -

%og p[H_,] = log

-10.83 -15.54 -28.52 -24.31 -

43.61 -6.36

-

-

M.l

-

M6 4.97 -7.70 -9.30 -15.0 -23.0 -13.94 -

M6 -2.19 -2.95 -5.67 -12.5 -21.6 -6.30 -

K, = -13.85.

TABLE 3 Some characteristics and execution times for the proposed test systems (In all cases errors in free and species concentrations are calculated) Program

System*

IOPb

No. of points

Total no. of iterations

Total no. of function evaluations

TimeC (s)

ES4ECl

Ia Ib Ic

0 0 0

12 12 12

18 (1.5)d 13 (1.1) 26 (2.2)

142 (11.8)d 125 (10.4) 173 (14.4)

30 (6) 30 (6) 30 (7)

ES4ECl

IIa IIb

0 1

41 10

85 (2.1) 50 (5.0)

176 (4.3) 120 (12.0)

151 (26) 45 (12)

ES4EC2

IIa

0

41

85 (2.1)

174 (4.2)

249 (150)

ES4ECle

IIIa IIIb

0 1

20 10

32 (1.6) 70 (7.0)

105 (5.3) 162 (16.2)

341 (200) 226 (150)

ES4ECl

IVa IVb

0 1

14 21

48 (3.4) 82 (3.9)

156 (11.1) 252 (12.0)

79 (23) 140 (50)

aIa: CA1 = 3.7 x lo-$. Ib: CA* = 10-6. Ic: CA1 = 10-l mol dm-‘. pH l-12 (constant increment). IIa: pH 3-9 (constant increment). IIb: titrant C~H = 0.75 moldme3, V, = 100 cm); 10 additions of 0.4 cm3; pH 2.7-9.9. IIIa: pH l-10.5 (constant increment). IIIb: titrant C 0~ = 0.1 mol dm-); V, = 100cm’; 9 additionsof 1.3cm3; pH 2.9-11.3. IVa: pH 3.5-10 (constant increment). IVb : titrant CoB = 0.25 mol dmp3; V,,= 100 cm3; 20 additions of 0.5 cm”; pH 2-9.2. Bad initial estimates of free concentrations were used, in order to check the converence ability of the algorithm. g 0, distribution diagrams; 1, titration curves. ‘Including printer time (ES4ECl) or disk recording (ES4EC2); in parentheses, net time. dMean values for each point. %ame results for ES4EC2 (execution time increased as in the case of System II).

394 TABLE Some

4 calculation

details

for System

V (cm’)

Ib

9

10

0.4

0.127 11.31

2.652 4.50

2.578 2.73

7.832 2.40

2.659 7.22

4.47

3.02

0.8

0.129 11.09

2.754 4.21

2.627 2.67

7.320 2.56

2.884 6.71

5.69

6.41

8.99

4.18

2.96

4.87

0.131 10.85

2.913 3.96

2.703 2.65

6.825 2.78

3.123 6.22

5.37

7.72

3.93

2.94

4.20

4.45

6.57

1.6

0.133 10.61

3.145 3.75

2.802 2.69

6.378 3.06

3.364 5.77

2.0

0.134 10.36

3.464 3.59

2.927 2.77

3.72

2.98

3.77

3.78

5.66

5.968 3.39

3.609 5.36

2.4

0.136 10.08

3.904 3.45

3.56

3.06

3.58

3.36

4.99

3.103 2.91

5.541 3.81

3.894 4.94

2.8

0.137 9.68

3.42

3.21

3.61

3.10

4.44

4.551 3.33

3.398 3.20

5.028 4.50

4.294 4.43

3.2

3.31

3.49

3.88

2.96

3.91

0.139 8.93

5.502 3.35

3.983 3.96

4.298 6.01

5.041 3.70

3.33

4.26

4.72

3.06

3.26

3.6

0.141 7.15

7.918 4.50

4.611 6.89

3.670 10.73

6.822 3.07

4.48

7.20

7.90

4.45

2.87

4.0

0.143 4.02

11.644 7.38

6.407 12.89

3.420 19.85

9.946 2.82

7.36

13.20

14.85

8.28

3.58

3 for concentrations.

bFirst

1.2

%ee

text for formation

2

3

IIba

constants

4

5

and Table

6

(1) IOH1 ; (2) IWclt)l; (3) W,(clt)l; (4) W,(clt)l, (5) 13) [Ni(cit)]; (14) [Ni(cit)H]; (15) [Ni(cit),]; (16) 6 See ref. 9. (For the program EQUIL, nrt = rife.))

7

8

row:

---log ck, k = 1

[Ni(cit),H];

(17)

...4

(7) W(cit)H,l; [Ni,(cit),H_,]; (18)

W(cit)l; (6) IK(cit)Hl;

tions are not very good. In Table 4, some detail of System IIb calculations are reported and Fig. 2 shows a distribution diagram (System IIa). Comparisons with the well known COMICS [ 51 and EQUIL [ 91 computer programs were also made. These two computer programs were chosen because they use two very different calculation methods, have been widely used in several laboratories, and are the progenitors of two families of other programs written for both general and specific problems (see e.g., [28] ). Moreover, the comparisons can be extended to other computer programs through the testing results of Leggett [22] and Garcia and Madariaga [33]. For the comparison between ES4ECl and COMICS, six test systems proposed by Leggett [22] were also used. In Table 5, some details of the comparisons between ES4ECl and other computer programs are reported. All the programs gave the same results within 0.01%. In the case of System I, both EQUIL and COMICS failed. For System II, COMICS failed

395

11

12

13

14

15

16

17

18

3.13

5.11

4.85

3.48

9.39

6.56

19.40

13.21

2.95

4.58

4.39

3.24

8.42

5.81

18.03

11.89

2.85

4.16

3.97

3.06

7.50

5.14

16.71

10.65

2.88

3.93

3.62

2.95

6.71

4.58

15.53

9.61

3.03

3.92

3.33

2.91

6.01

4.13

14.47

8.74

3.33

4.08

3.08

2.95

5.33

3.74

13.40

7.93

3.86

4.47

2.87

3.13

4.60

3.41

12.16

7.05

5.03

5.29

2.72

3.73

3.73

3.28

10.38

5.83

8.40

7.31

2.72

5.51

3.10

4.43

6.82

3.86

15.00

11.39

4.27

10.18

4.40

8.85

3.66

2.63 Total

n&

nPe

n& (EQUILd)

12

37

22

5

7

6

5

7

6

3

5

5

3

5

5

3

5

I

3

5

13

4

6

12

7

9

10

7

11

8

52

97

94

(Eqn. l’);(l) K+;(2)Cu’*; (3) Ni’+; (4) cit3-; (5) H’. Second row: -log Xi, i = 1 . . . 18 (Eqn. 1’); (8) [Cu,(cit),l; (9) [Cu,(cit),H_,]; (10) [Cu,(cit),H_,]; (11) [Cu(cit)H]; (12) [Cu(cit)H_,l; [Ni(cit),H_,] (charges omitted). ‘nit, number of iterations; nf,, number of function evaluations.

at pH > 8. ES4ECl was always faster than the other two programs. [It should be observed that the EQUIL program (Newton-Raphson technique) is more efficient for dealing with a few mass-balance equations whilst the COMICS program is more efficient for dealing with large systems with uncomplicated equilibria (such as System III)]. On average, ES4ECl is five times faster than COMICS and three times faster than EQUIL (see Table 5). By considering the results of Leggett [22] and Garcia and Madariaga [33], ES4ECl is probably faster than the majority of other programs in this field. Future developments of this work may be considered. First of all, it is necessary to find an algorithm for the comparison of execution times among various machines. In fact, comparisons of various computer programs need an objective method for avoiding dependence on the computer used; then several computer programs could be compared, even in different laboratories. Although test routines are available, for the present problem it would

396 TABLE 6 Comparisons for execution times among different computer programs for the calculation of equilibrium concentrations System

IaC IIa IIb IIIa IIIb IVa IVb I’d II’ III’ IV’ V’ VI’ I’-vIe

Calculation times (8) ES4ECl*

EQUILb

COMICSb

4 11 6 22 43 10 19 9 7 9 8 10 7 -

failed -

failed 85 -

-

15

141 -

42

100

-

73

failed 59 13 25 82 57 10 307

Other computer programs

Leggett results [ 22 ]

Garcia and Madariaga results [ 33 ]

COMICS COMPLEX [19] EQUIL COMPLEX-80 [32] HALTAFALL [ 4 ] SOLGASWATER [ 271

lOOf 53 33 -

lOOf -

-

110 103 104

aNet time, without calculating errors. bThese computer programs were run in doubleprecision arithmetic. ‘Systems I-IV: see text and Table 3. dSystems I’-VI’: see Leggett [22]. eAverage results of Leggett [22] for systems I’-VI’; arbitrary times: r (EQUIL) = 100. *Average results of Leggett [22] and Garcia and Madariaga [33] (in both cases, six relatively simple systems were tested).

be better to build up specific algorithms that take into account the features of the problem under study. Further developments of the ES4EC program will be: (a) the use of formation constants at various ionic strengths and the simulation of titration curves at variable ionic strength; (b) the possibility of plotting results in a graphical form; (c) the possibility of dealing with more than one phase; (d) an analysis of various factors to be considered in the calculations relative to different systems e.g., hints must be given for choosing initial estimates and p values on the basis of the type of system; (e) the possibility of solving specific problems with a high degree of efficiency. Some of these developments are already under way in this laboratory. Results on these investiga-

397

tions, together with definitive versions of the program ES4EC, will soon be reported. We thank CNR (Rome) financial support.

and the Minister0 della Pubblica Istruzione for

REFERENCES 1 W. B. White, S. M. Johnson and G. B. Dantzig, J. Chem. Phys., 28 (1958) 751. 2 J. W. Swinnerton and W. W. Miller, J. Chem. Educ., 36 (1959) 485. 3 A. J. Bard and D. M. King, J. Chem. Educ., 42 (1965) 127. 4 N. Ingri, W. Kakolowicz, L. G. Sillen and B. Warnqvist, Talanta, 14 (1967) 1261; B. Elgquist, Talanta, 16 (1969) 1502; R. Ekelund, L. G. Sillen and 0. Wahlberg, Acta Chem. Stand., 24 (1970) 3073. 5 D. D. Perrin and I. G. Sayce, Talanta, 14 (1967) 833. 6 E. S. Johansen, Acta Chem. Stand., 21 (1967) 2273. 7 S. Gordon, Complex Chemical Equilibrium Calculations, NASA Spec. Publ. 239, 1970. 8 F. Van Zeggeren and H. S. Storey, The Computation of Chemical Equilibria, Cambridge University Press, Cambridge, 1970. 9 Ting PO I and G. H. Nancollas, Anal. Chem., 44 (1972) 1940. 10 I. K. Karpof and L. A. Kaz’min, Geochem. Int., 9 (1972) 252. 11 M. Bos and H. Q. J. Meershoek, Anal. Chim. Acta, 61 (1972) 185. 12 F. Mare1 and J. Morgan, Environ. Sci. Technol., 6 (1972) 58. 13 R. E. McDuff and F. M. M. Morel, Description and Use of the Chemical Equilibrium Program REDEQLB, W. M. Keck Laboratory of Environmental Engineering Science, California Institute of Technology, Pasadena, CA, 1973. 14 H. S. Dunsmore and D. Midgley, Anal. Chim. Acta, 72 (1974) 121. 15 A. H. Truesdell and B. F. Jones, J. Res. U. S. Geol. Surv., 2 (1974) 233. 16 D. A. Crerar, Geochim. Cosmochim. Acta, 39 (1975) 1375. 17 L. J. Walters, Jr., and T. J. Walery, Comput. Geosci., 1 (1975) 57. 18 L. N. Plummer, D. L. Paskhurst and D. R. Kosiur, MIXB: A Computer Program for Modeling Chemical Reactions in Natural Waters, U. S. Geol. Survey, Water-Resources Investigations, Reston, VA, 1975, pp. 75-161. 19 G. Ginzburg, Talanta, 23 (1976) 149. 20 R. Maggiore, S. Musumeci and S. Sammartano, Talanta, 23 (1976) 43. 21 D. J. M. Park, J. Chem. Phys., 65 (1976) 3085. 22 D. J. Leggett, Talanta, 24 (1977) 535. 23 P. M. May,P. W. Linder and D. R. Williams, J. Chem. Sot. Dalton Trans., (1977) 588. 24 G. Arena, C. Rigano and S. Sammartano, Ann. Chim. (Rome), 68 (1978) 693. 25 J. J. Fardy and R. N. Sylva, SIAS. A Computer Program for the Generalized Calculation of Speciation in Mixed Metal-Ligand Aqueous Systems, AAEC/E445, Lucas Heights, Australia, 1978. 26 J. W. Ball, E. A. Jenne and D. K. Nordstrom, in E. A. Jenne (Ed.), Chemical Modeling in Aqueous Systems, ACS Symp. Ser. No. 93, Am. Chem. Sot., Washington, D.C., 1979. 27 G. Eriksson, Anal. Chim. Acta, 112 (1979) 375. 28 D. K. Nordstrom et al., in E. A. Jenne (Ed.), Chemical Modeling in Aqueous Systems, ACS Symp. Ser. No. 93, Am. Chem. Sot., Washington, D.C., 1979; D. K. Nordstrom and J. W. Ball, in C. J. M. Kramer and J. C. Duinker (Eds.), Complexation of Trace Metals in Natural Waters, M. Nijhoff and W. Junk, The Hague, 1984. 29 L. Avery, J. Chem. Phys., 76 (1982) 3242.

398 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52

53 54 55 56

D. D. Clark and S. M. Schuster, Computers and Chem., 7 (1983) 25. V. S. Tripathi, Talanta, 30 (1983) 65. J. M. Madariaga and A. Garcia, Comput. and Chem., 8 (1984) 187. A. Garcia and J. M. Madariaga, Comput. Chem., 8 (1984) 193. J. Kostrowicki and A. Liwo, Comput. Chem., 8 (1984) 91. F. B. Hildebrand, Introduction to Numerical Analysis, McGraw-Hill, New York, 1956. L. C. W. Dixon, Nonlinear Optimization, The English Universities Press, London, 1972. D. D. Perrin, Nature, 206 (1965) 170. P. S. Hallman, D. D. Perrin and A. E. Watt, Biochem, J., 121 (1971) 549. D. Dyrssen, D. Jagner and F. Wengelin, Computer Calculation of Ionic Equilibria and Titration Procedures, Almqvist, Stockholm, 1968. T. Anfalt and D. Jagner, Anal. Chim. Acta, 47 (1969) 57 ; 57 (1971) 165. E. A. Jenne (Ed.), Chemical Modeling in Aqueous Systems, ACS Symp. Ser. No. 93, Am. Chem. Sot., Washington, D.C., 1979. C. Rigano, E. Rizzarelli and S. Sammartano, Thermochim. Acta, 33 (1979) 211. G. Arena, E. Rizzarelli, S. Sammartano and C. Rigano, Talanta, 26 (1979) 1. P. G. Daniele, A. De Robertis, S. Sammartano and C. Rigano, Ann. Chim. (Rome), 73 (1983) 619. P. G. Daniele, C. Rigano and S. Sammartano, Thermochim. Acta, 62 (1983) 101. P. G. Daniele, A. De Robertis, S. Sammartano and C. Rigano, Thermochim. Acta, 72 (1984) 305. A. De Robertis, C. Rigano and S. Sammartano, Thermochim. Acta, 74 (1984) 343. P. G. Daniele, G. Ostacoli, C. Rigano and S. Sammartano, Transition Met. Chem., 9 (1984) 385. A. De Robertis, C. De Stefano, R. Scarcella and C. Rigano, Thermochim. Acta, 80 (1984) 197. A. De Robertis, C. De Stefano, S. Sammartano, R. Scarcella and C. Rigano, J. Chem. Res., 42(S) (1985) 629 (M). P. G. Daniele, A. De Robertis, C. De Stefano, S. Sammartano and C. Rigano, J. Chem. Sot. Dalton Trans., (1985) 2356. P. Amico, G. Arena, P. G. Daniele, G. Ostacoli, E. Rizzarelli and S. Sammartano, in K. J. Irgolic and A. E. Martell (Eds.), Environmental Inorganic Chemistry, VCH Publ., Deerfield Beach, FL, U.S.A., 1985. I. S. Berezin and N. P. Zhidkov, Computing Methods, Pergamon, Oxford, 1965. I. Nagypal, I. Paka and L. Zekani, Talanta, 25 (1978) 549. W. E. Wentworth, J. Chem. Educ., 42 (1965) 96,162. M. Micheloni, A. Sabatini and A. Vacca, Inorg. Chim. Acta, 25 (1977) 41.