% Calculate Equilibrium Reactions % Moran Bercovici, v1.0 May 2008: First version % Moran Bercovici, v1.1 Feb 2009: Comments added function ITPCalculator clear all; close all; clc; format short g global F Rmu Temp Kw muH muOH global met2lit % INPUT DATA %-------------------------------------------------------------------------- %-------------------------------------------------------------------------- % Input sructure consists of [valence, mobility, pKa] triplets % the order or triplets doesn't matter (ok to go from z=-1 to z=2 or vice versa). % z mu [m^2/Vs] pKa INP={ [ -1 42.4E-9 4.756 ]; % ACETIC ACID (S) [ 1 29.5E-9 8.076 ]; % TRIS [ 2 59.2E-9 1.82 1 29.6E-9 6.04 -1 28.3E-9 9.33 ]; % Histidine [ 1 51.9E-9 13.7 ]; % Sodium }; % Total concentration for each species cTot=[0.05,0.1,0.2,0.01]; % Total concentration [M] %-------------------------------------------------------------------------- %-------------------------------------------------------------------------- % Define constants F=9.65E4; % Faraday's const.[C/mol] Rmu=8.31; % Universal gas const. [J/mol*K] Temp=300; % Temperature [K] Kw=1E-14; % Water equilibrium constant muH=362E-9/F; % Mobility of Hydronium % [m^2/s*V]/F --> [mol*s/Kg] muOH=205E-9/F; % Mobility of Hydroxide % [m^2/s*V]/F --> [mol*s/Kg] met2lit=1000; % Call Equilibrium Solution routine [cH,zMat,muMat,muEffVec,cizMat]=SolveEquilibrium(INP,cTot); % Print results to screen fprintf('pH = %5.5f (cH = %2.2e [M])\n',-log10(cH),cH); % Print pH fprintf(['cTotal [M]\t\t|\t',repmat('z\tc [M]\t\t\t\t',1,size(cizMat,2)),'\n']); % Print table header fprintf([repmat('-----------',1,15),'\n']); % loop on all species for ij=1:size(cizMat,1) % Create vector to be displayed DispVec(1:2:2*size(cizMat,2)-1)=zMat(ij,:); % Every odd entry will show the valence DispVec(2:2:2*size(cizMat,2))=cizMat(ij,:); % Every even entry will show the concentration corresponding to the valence % print: total concentration, valence, concentration, effective mobility fprintf(['%2.2e\t\t|\t',repmat('%d\t%2.2e\t\t\t',1,size(cizMat,2)),' Eff Mobility = %5.2e [m^2/Vs]\n'],cTot(ij),DispVec,muEffVec(ij)); end % for ij %-------------------------------------------------------------------------- function [cH,zMat,muMat,muEffVec,cizMat]=SolveEquilibrium(INP,cTot) global F Rmu Nspecies Temp Kw % INPUT: INP - Structure, Containing the valence, mobility and pKa of each species % cTot - Vecotr, containing the concentration [mol/lit] of each species % OUTPUT: cH - Scalar, concentration of hydronium ions ( pH = -log10(cH) ) % zMat - Matrix, every row corresponds to a difrernt species and contains valence values for that species % muMat - Matrix, every row corresponds to a difrernt species and contains mobility values for that species % muEff - Vector, effective mobility values for all species % cizCube - Matrix, concentration values. Every row correponds to a different species. % every columns corresponds to a differnt valence within that % specie. Nspecies=size(INP,1); % PREPARE L Matrix %------------------ % Calculate the number of columns required for the matrix. % this is determined by the maximum value of (p_i - n_i) for all species % p_i - n_i = difference between most positive and most negative valence. MaxCol=-Inf; for j=1:size(INP,1) MaxCol=max([MaxCol,max(INP{j}(1:3:end))-min(INP{j}(1:3:end))+1]); end % Initialize matrices to zero LMat=zeros(size(INP,1),MaxCol); zMat=LMat; muMat=LMat; KaMat=LMat; DMat=LMat; % Loop on species for j=1:Nspecies zList=INP{j}(1:3:end); % Get list of valences muList=INP{j}(2:3:end)./(F*abs(zList)); % Get list of normalized mobilities pKaList=INP{j}(3:3:end); % Get list of pKa KaList=10.^(-pKaList); % Create list of Ka DList=Rmu*Temp*muList; % Creale list of diffusivities using Nernst-Einstein relation [zList,Index]=sort(zList); % Sort the valence in increasing order, and get indices KaList=KaList(Index); % Use indices to sort the other list in a consistent order muList=muList(Index); DList=DList(Index); Ip1=find(zList==1); Im1=find(zList==-1); % Find indices where the valence value is +1 and -1 (at least one of these valences must always for any species used). zList = [zList(1:Im1),0,zList(Ip1:end)]; % Add the neutral state to the list of valence, between negative and positive values muList = [muList(1:Im1),0,muList(Ip1:end)]; % Add the value of mobility corresponding to the neutral state (zero) DList = [DList(1:Im1),mean(DList),DList(Ip1:end)]; % Assign the neutral state a diffusivity which is the average of all other diffusivities KaList = [KaList(1:Im1),1,KaList(Ip1:end)]; % 1 added to the list of Ka values for math reasons (see formulation for 'L' matrix) zMat(j,1:length(zList))=zList; % Put all lists into corresponding matrices. muMat(j,1:length(muList))=muList; % Each row in the matrix corresponds to a different species KaMat(j,1:length(KaList))=KaList; DMat(j,1:length(DList))=DList; zListArranged{j}=zList; % Get minimum and maximum valences for this species % and construct the matrix L nj=min(zList); pj=max(zList); for z=zList if z<0 LMat(j,z-nj+1)=prod(KaList(z-nj+1:-nj)); elseif z>0 LMat(j,z-nj+1)=1/prod(KaList(-nj+2:z-nj+1)); elseif z==0 LMat(j,z-nj+1)=1; end %if end % for z end %for ij % CONSTRUCT POLYNOMIALS %-------------------- % Construct the polynomial Q by multiplying all the polynomials in the % matrix L (all rows). Multiplying polynomials is equivalent to convolving % their vectors of coefficients Q1=1; for j=1:size(LMat,1) Q1=conv(Q1,LMat(j,:)); end %for j Q2=[-Kw 0 1]; Q=conv(Q1,Q2); % Construct the polynomials Pi % Loop on all species for i=1:Nspecies LMatMod=LMat; % Defined Modified L Matrix, initially identical to L Matrix LMatMod(i,:)=LMat(i,:).*zMat(i,:); % Modify just the row i % convolve all rows in the LMatMod to construct the polynomial Pi=1; for kl=1:size(LMatMod,1) Pi=conv(Pi,LMatMod(kl,:)); end %for j Pi=conv([0 1],Pi); PMat(i,:)=Pi; % Insert all polynomials Pi as rows in the matrix PMat end %for i cTotMat=repmat(cTot',1,size(PMat,2)); % All the enries in each row are identical and equal to the total concentration of that specie P=sum(cTotMat.*PMat,1); % The polynomial P is the sum of c_i * P_i Poly=PolSum(P,Q); % Construct the final polynomial as the sum of P and Q roo=roots(fliplr(Poly)); % Find all roots of the polynomial roo=roo(imag(roo)==0); % Eliminate complex solutions cH=roo(roo>0); % Choose the positive solution. Note that cH is in mol/lit PolDeg=size(LMat,2); % Polynomial degree cHPolVec=[1;cumprod(ones(PolDeg-1,1)*cH,1)]; % Create vector of cH Powers (cH^0, cH^1, ... , cH^(PolDeg-1)) cHMatPower=repmat(cHPolVec',[Nspecies,1]); % Replicate powers vector to create matrix M1Mat=repmat(cTot'./(LMat*cHPolVec),[1,PolDeg]); % Temporarty matrix (see formulation for c_i_z) cizMat=LMat.*cHMatPower.*M1Mat; % Calculate concentration of each ionic state (valence) of each species muEffVec=F*sum(cizMat.*muMat.*zMat,2)./cTot'; % Calculate effective mobility % Sum arbitrary length polynomials %----------------------------------------------------- function out=PolSum(P1,P2) % This essentially pads the shorter vecotr with zeros, and then sums the % two vectors Psize=max(length(P1),length(P2)); % get the length of the largest polynomial P1out=zeros(1,Psize); P2out=P1out; % Set temporary vectors of size Psize to zero. P1out(1:length(P1))=P1; % Inject P1 values into temporar vector P2out(1:length(P2))=P2; % Inject P2 values into temporar vector out=P1out+P2out; % sum both vectors