“Difficult” problems

1. “Difficult” problems Marcel Maeder Department of Chemistry University of Newcastle Australia Marcel.maeder@newcastle.edu.au

2. A “simple” problem What is the pH of a 1M solution of acetic acid?

3. What is the pH of a 1M solution of acetic acid? 2 equations and 3 unknowns: [AH], [A], [H]

4. Approximation 1: Acetic acid is a weak acid, only a very small fraction dissociates and acetic acid is the only source of protons

5. Approximation 2: acetic acid is the only source of protons

6. But, there are two sources of protons, acetic acid and the autodissociation of water

7. Activities In reality the law of mass action really applies to activities and not concentrations. More on that on Thursday

8. The “simple” problem of calculating the pH of a 1M acetic acid involves the solution of a cubic equation. This is possible, but the equation for that solution is very long and nobody ever uses it. Things are much more difficult for polyprotic acids. We need a general solution!

9. nomenclature X, Y, Z are components XxYyZz are species

10. Example: Cu + ethylene diamine, in water Species: M L H LH LH2 ML … Components: M (Cu) L (en) H

11. 3 components 11 species

12. The task is to compute all species concentrations, knowing the total concentrations of the components and the equilibrium constants There are 11 unknowns (the species concentrations) We need 11 equations! 3 equations for the total concentrations of the components:

13. 8 equations for the species concentrations

14. Generalised equations: Allow very compact MATLAB code: c_spec =beta.*prod(repmat(c',1,nspec).^Model,1); %species conc c_tot_calc=sum(Model.*repmat(c_spec,ncomp,1),2)';%comp ctot calc

15. How can we write a MATLAB program that can resolve these systems of many unknowns with many equations and that can deal with any model? The Newton-Raphson algorithm

16. The Taylor series expansion

17. f(x) f’(x) f(x+x) 0 x x+x 20 22 24 26 28 30 32 34

18. The function to minimise is: The difference between the known total concentrations and the calculated values:

20. The Jacobian J

25. Kinetics at non-constant temperature • Kinetics and equilibria, taking into account changes in ionic strength • Kinetics at non-constant pH

26. Why? • In traditional measurements temperature, ionic strength, pH are kept constant using buffers, excess inert salt, thermostatting • All these external means are only required to simplify the computations. If it is possible to take changes in these parameters into account, buffers etc are no more required

27. Kinetics at non-constant temperature Example: Differential equations cannot be solved explicitly → numerical integration: Euler, Runge-Kutta, MATLAB ODE solvers

28. ode_AeqBtoC.m function c_dot=ode_AeqBtoC(t,c,flag,k) % A <-> B % B --> C c_dot(1,1)=-k(1)*c(1)+k(2)*c(2); % A_dot c_dot(2,1)= k(1)*c(1)-k(2)*c(2)-k(3)*c(2); % B_dot c_dot(3,1)= k(3)*c(2); % C_dot Translation into MATLAB language, the function ode_AeqBtoC.m is used by the MATLAB ode-solvers This function computes the derivatives of the concentrations vs. time (c_dot) at time t for a given set of concentrations c. The rest is done by the ODE-solver

29. Fast first reaction and slow second reaction Difficult to observe both in one measurement AeqBtoC.m c0=[1;0;0]; % initial conc of A, Cat, B and C k=[.2;.1;.01]; % rate constants k1 and k2 times=[0:1:300]'; [times,C] = ode45('ode_AeqBtoC',times,c0,[],k); % call ode-solver figure(1); plot(times,C) % plotting C vs t xlabel('time');ylabel('conc.');

30. Arrhenius equation Rate constants are temperature dependent, only constant at constant temperature. Arrhenius equation describes rate constant k as a function of the temperature (Eyring equation could be used as well)

31. slope: -Ea/R ln(k) 1/T Arrhenius plot

32. ode_AeqBtoC_T.m function c_dot=ode_AeqBtoC_T(t,c,flag,k,temp,A,Ea,times) % A <-> B % B --> C R=8.314; % gas constant J K-1 mol-1 T=lolipop(times,temp,t,2,5); % interpolation to comp T at particular time t k=A.*exp(-Ea./(R*(T+273))); % rate constants at T c_dot(1,1)=-k(1)*c(1)+k(2)*c(2); % A_dot c_dot(2,1)= k(1)*c(1)-k(2)*c(2)-k(3)*c(2); % B_dot c_dot(3,1)= k(3)*c(2); % C_dot A function is required that computes the derivatives of the concentrations vs. time, at one particular time t At that time the rate constants are computed using the Eyring parameters, based on the temperature T at that time T has to be interpolated from the temperatures temp measured at the times times

33. AeqBtoC.m … contiuned A =[12.13 1.77 36.79]; % pre-exp factor for each rate constant Ea=[1e4 7e3 2e4]; % activation energies temp=20+0.2*times; % temperatures at measurement times [times,C] = ode45('ode_AeqBtoC_T',times,c0,[],k,temp,A,Ea,times); % ode-solver figure(2); plot(times,C) % plotting C vs t xlabel('time');ylabel('conc.');

34. at 20°C at 20-80°C Similar to temperature program in GC or non-isocratic mobile phase in HPLC

35. Modeling - Fitting • The core of the fitting algorithm is the modeling of the measurement and in most cases this is the computation of the matrix C • This is exactly what we have done so far for non-isothermal kinetics • The parameters used to model the concentration profiles can, at least potentially, also be fitted to appropriate data • In this case non-isothermal measurements allows the determination of the activation parameters of all rate constants • Much more convenient than classical method of repeating the measurement at different temperatures with subsequent Arrhenius analysis

36. Kinetics, taking into account changes in ionic strength

37. Kinetics, taking into account changes in ionic strength {A} = activity of species A

38. For Ionic compounds, activity coefficients  can be approximated in dilute solutions as (Debye-Hückel): With A parameter depending on dielectric constant of solvent, in water A~0.51 zi charge of species  Ionic strength of solution, computed as: Activity coefficients

39. As the concentrations change during the reaction, the ionic strength and all activity coefficients change as well.

40. Example → + 1M Fe(ClO4)3 1M Na2SO4 1M Na+ 0.5 M SO42- 0.5 M Fe3+ 1.5 M ClO4-

41. 1M Na+ 0.5 M SO42- 0.5 M Fe3+ 1.5 M ClO4- Note: Debye-Hückel equation not valid at such high concentrations !

42. The kinetics and the equilibrium position are dramatically affected by the very small activity coefficients !

43. c_dot=ode_AplusBeqC_I.m function c_dot=ode_AplusBeqC_I(t,c,flag,k,c_I,charges,mode,A) % A + B <--> C including activities [gamma,mu] = act_coeffs_I([c;c_I],charges,mode,A); % act coeff and ionic strength act=c.*gamma(1:length(c)); % activities instead of conc. % c_dot(1,1)=- k(1)*c(1)*c(2)+k(2)*c(3); % A_dot ignoring activities c_dot(1,1)=-k(1)*act(1)*act(2)+k(2)*act(3); % A_dot c_dot(2,1)= c_dot(1,1); % B_dot c_dot(3,1)= -c_dot(1,1); % C_dot function[gamma,mu] = act_coeffs (conc,charges,mode,A) % calculating activity coefficients mu = sum(1/2*(conc.*(charges.^2))); ifmode==1 gamma = ones(size(conc)); elseifmode==2 log_gamma = (-A*(charges.^2)*(mu^0.5))/(1+(mu^0.5)); % defining log g gamma = 10.^log_gamma; % defining activity coefficient, gamma end

44. AplusBeqC.m % AplusBeqC % A + B <--> C c0=[.01;.01;0]; % initial conc of Fe, SO4 and Fe(SO4) k=[10;.01]; % rate constants k1 and k2 times=[0:1:300]'; [times,C] = ode45('ode_AplusBeqC',times,c0,[],k); % call ode-solver figure(1); plot(times,C) % plotting C vs t xlabel('time');ylabel('conc.'); c_I=[2*c0(1); 3*c0(2)]; % Na ClO4 charges=[3; -2; 1; 1; -1]; % Fe SO4 Fe(SO4) Na ClO4 mode=2; A=0.51; [times,C] = ode45('ode_AplusBeqC_I',times,c0,[],k,c_I,charges,mode,A);, % ode-solver figure(2); plot(times,C) % plotting C vs t xlabel('time');ylabel('conc.'); 0.01M Na2SO4 0.01M Fe(ClO4)3 k1=10 k2=0.01 without with ionic strength Kinetics is much slower but also equilibrium position very different !

45. Equilibrium calculations are more complex:In kinetics, using an ODE solver, the only requirement is to compute the derivatives of the concentrations vs. time for a given set of concentrations. Thus, it is possible to calculate the activity coefficients at that moment and incorporate into the derivatives.In equilibria, the activity coefficients depend on the concentrations and the concentrations depend on the activities. An internal iterative process is required.

46. Titration of phosphoric acid:

47. Iterative approach: Guess -values Calculate concentrations Calculate new -values

48. Newton-Raphson algorithm The core is the computation of the difference d between actual and computed total concentrations c, and The derivatives of these with respect to the total concentrations

49. OLD c_spec =beta.*prod(repmat(c',1,nspec).^Model,1); % species conc c_tot_calc=sum(Model.*repmat(c_spec,ncomp,1),2)'; % comp ctot calc d =c_tot-c_tot_calc; % diff actual and calc total conc NEW, with activity coefficients it=1; while it<99 act_comp = c.*gamma(1:ncomp); act_spec = beta.*prod(repmat(act_comp',1,nspec).^Model,1); % species activities c_spec = act_spec./gamma(1:nspec); % species activities [gamma_new,mu] = act_coeffs_I([c_spec c_tot_I],charges,mode,A); % calc gamma for next it if all(abs(gamma-gamma_new)<1e-3) break end gamma=gamma_new; it=it+1; end c_tot_calc = sum(Model.*repmat(c_spec,ncomp,1),2)'; % comp tot d = c_tot-c_tot_calc; % diff actual and calc total conc

50. no yes act_comp = c.*gamma(1:ncomp); act_spec = beta.*prod(repmat(act_comp',1,nspec).^Model,1); c_spec = act_spec./gamma(1:nspec); [gamma_new,mu]= act_coeffs_I([c_spec c_tot_I],charges,mode,A); if all(abs(gamma-gamma_new)<1e-3)