source: src/main/java/agents/anac/y2015/agentBuyogV2/flanagan/physchem/DonnanConcn.java

Last change on this file was 127, checked in by Wouter Pasman, 6 years ago

#41 ROLL BACK of rev.126 . So this version is equal to rev. 125

File size: 22.8 KB
Line 
1/*
2* Classes DonnanConcn
3*
4* Class DonnanConcn impliments the interface
5* MinimisationFunction to the Class Minimisation
6* and contains a function that calculates the ionic
7* concentrations if the ionophore binds more than
8* one species of ion. The sum of the squares of the
9* equilibrium functions is minimised.
10*
11* WRITTEN BY: Michael Thomas Flanagan
12*
13* DATE: November 2004
14* UPDATE: 6 December 2004
15*
16* DOCUMENTATION:
17* See Michael Thomas Flanagan's JAVA library on-line web page:
18* http://www.ee.ucl.ac.uk/~mflanaga/java/Donnan.html
19* http://www.ee.ucl.ac.uk/~mflanaga/java/
20*
21* Copyright (c) November 2004
22*
23* PERMISSION TO COPY:
24* Permission to use, copy and modify this software and its documentation for
25* NON-COMMERCIAL purposes is granted, without fee, provided that an acknowledgement
26* to the author, Michael Thomas Flanagan at www.ee.ucl.ac.uk/~mflanaga, appears in all copies.
27*
28* Dr Michael Thomas Flanagan makes no representations about the suitability
29* or fitness of the software for any or for a particular purpose.
30* Michael Thomas Flanagan shall not be liable for any damages suffered
31* as a result of using, modifying or distributing this software or its derivatives.
32*
33***************************************************************************************/
34
35package agents.anac.y2015.agentBuyogV2.flanagan.physchem;
36
37import agents.anac.y2015.agentBuyogV2.flanagan.math.*;
38
39class DonnanConcn implements MinimisationFunction{
40
41 public int numOfIons = 0; // number of ionic species
42 public double[] concnA = null; // concentration of ions in partition A (moles per cubic metre)
43 public double[] concnB = null; // concentration of ions in partition B (moles per cubic metre)
44 public double[] molesT = null; // total moles of an ion
45 public double[] complex = null; // concentration of complex in partition B (moles per cubic metre)
46 public double[] excessConcnA = null; // excess concentration of ions in partition A in interfacial charge (moles per cubic metre)
47 public double[] excessConcnB = null; // excess concentration of ions in partition B in interfacial charge (moles per cubic metre)
48 public double[] excessComplex = null; // excess concentration of complex in partition B in interfacial charge (moles per cubic metre)
49 public double[] assocConsts = null; // association constants of the ions with the ionophore (cubic metres per mole)
50 public double[] partCoeffPot = null; // partition coefficients in presence of inter-partition, e.g. Donnan, potential
51 public int[] indexK = null; // ion index of ionic species with non zero assocConsts
52 public int nonZeroAssocK = 0; // number of ionic species with affinity for the ionophore
53 public double[] radii = null; // radii of ions
54 public double[] charges = null; // charge of ions
55 public double ionophoreConcn = 0.0D; // ionophore concentration - entered as molar, converted to moles per cubic metre
56 public double ionophoreRad = 0.0D; // ionophore radius (metres)
57 public double volumeA = 0.0D; // volume of partition A (cubic metres)
58 public double volumeB = 0.0D; // volume of partition B (cubic metres)
59 public double interfacialArea = 0.0D; // onterfacial area between Partitions A and B(square metres)
60 public double epsilonA = 0.0D; // relative electrical permittivity of partition A
61 public double epsilonB = 0.0D; // relative electrical permittivity of partition B
62 public double epsilonSternA = 0.0D; // relative electrical permittivity of partition A Stern layer
63 public double epsilonSternB = 0.0D; // relative electrical permittivity of partition B SternLayer
64 public double temp = 25.0-Fmath.T_ABS; // Temperature (degrees Kelvin) [Enter temperature in degrees Celsius]
65 public double diffPotentialA = 0.0D; // Double layer potential difference - compartment A [volts]
66 public double diffPotentialB = 0.0D; // Double layer potential difference - compartment B [volts]
67 public double sternPotential = 0.0D; // Stern potential difference [volts]
68 public double sternCap = 0.0D; // Stern layer capacitance [F]
69 public double sternDeltaA = 0.0D; // Stern layer thickness in partition A [m]
70 public double sternDeltaB = 0.0D; // Stern layer thickness in partition B [m]
71 public double chargeValue = 0; // Absolute value of the charge valency if all are the same, i.e. symmetric electrolytes of the same valency
72 public boolean chargeSame = true; // = false if chargeValue not the same for all ions
73 public double interfacialChargeDensity = 0.0D;// interface Charge Density [C per square metre]
74 public double interfacialCharge = 0.0D; // interface Charge [C]
75 public boolean includeIc = true; // = true - interface charge included in the calculation of the Donnan potential
76 // = false - interface charge ignored in the calculation of the Donnan Poptential
77 public double potential = 0; // current estimate of Donnan potential
78 private double penalty = 1.0E+50; // penalty function multiplier invoked if an estimated ion concentration is negative
79
80 // method that returns sum of squares of equilibrium equations
81 // x[0] transfers the current estimate of the concentrations in partition B.
82 public double function(double[] x){
83
84 double sumOfSquares = 0.0D;
85
86 // Calculate estimate of complex concentration
87 if(this.nonZeroAssocK>0 && this.ionophoreConcn>0.0D){
88 if(this.nonZeroAssocK==1){
89 complex[indexK[0]] = this.assocConsts[indexK[0]]*x[indexK[0]]*this.ionophoreConcn/(1.0D + this.assocConsts[indexK[0]]*x[indexK[0]]);
90 }
91 else{
92 double[] vec = new double[this.nonZeroAssocK];
93 double[][] mat = new double[this.nonZeroAssocK][this.nonZeroAssocK];
94
95 // set up simultaneous equations
96 for(int i=0; i<this.nonZeroAssocK; i++){
97 vec[i] = this.assocConsts[indexK[i]]*x[indexK[i]]*this.ionophoreConcn;
98 for(int j=0; j<this.nonZeroAssocK; j++){
99 mat[i][j] = this.assocConsts[indexK[i]]*x[indexK[i]];
100 if(i==j)mat[i][j] += 1.0D;
101 }
102 }
103
104 // solve simultaneous equations to obtain complex concentrations
105 Matrix matrix = new Matrix(mat);
106 vec = matrix.solveLinearSet(vec);
107 for(int i=0; i<this.nonZeroAssocK; i++){
108 this.complex[indexK[i]] = vec[i];
109 }
110 }
111 }
112
113 // Test to check whether interface charge is included in the calculation of the Donnan potential
114 if(this.includeIc){
115 // Calculate the excess charge in the interface region
116 double excess = Math.abs(this.interfaceCharge(x, this.potential));
117
118 // Calculate concentrations of ions involved in excess interfacial charge region
119 excessConcentrations(x, excess, this.potential);
120
121 // Calculate sum of squares of extended equilibrium functions
122 for(int i=0; i<this.numOfIons; i++){
123 double aa = x[i]*(this.volumeB + this.partCoeffPot[i]*this.volumeA) + this.excessConcnA[i]*this.volumeA + (this.excessConcnB[i] + this.complex[i] +this.excessComplex[i])*this.volumeB - this.molesT[i];
124 sumOfSquares += aa*aa;
125 if(x[i]<0.0D)sumOfSquares += x[i]*x[i]*this.penalty;
126 }
127 }
128 else{
129 // Calculate sum of squares of equilibrium functions
130 for(int i=0; i<this.numOfIons; i++){
131 double aa = x[i]*(this.volumeB + this.partCoeffPot[i]*this.volumeA) + this.complex[i]*this.volumeB - this.molesT[i];
132 sumOfSquares += aa*aa;
133 if(x[i]<0.0D)sumOfSquares += x[i]*x[i]*this.penalty;
134 }
135 }
136 return sumOfSquares;
137 }
138
139
140 // Calculates the concentrations of ions involved in excess interfacial charge region
141 public void excessConcentrations(double[] x, double excess, double potential){
142
143 if(potential==0.0D){
144 for(int i=0; i<this.numOfIons; i++){
145 this.excessConcnA[i] = 0.0D;
146 this.excessConcnB[i] = 0.0D;
147 this.excessComplex[i] = 0.0D;
148 }
149 }
150 else{
151 double sumA = 0.0D;
152 double sumB = 0.0D;
153 double sumC = 0.0D;
154 for(int i=0; i<this.numOfIons; i++){
155 if(potential>0.0D){
156 if(this.charges[i]>0.0D){
157 sumB += x[i]*Math.abs(charges[i]);
158 sumC += this.complex[i]*Math.abs(charges[i]);
159
160 }
161 else{
162 sumA += x[i]*partCoeffPot[i]*Math.abs(charges[i]);
163 }
164 }
165 else{
166 if(this.charges[i]<0.0D){
167 sumB += x[i]*Math.abs(charges[i]);
168 sumC += this.complex[i]*Math.abs(charges[i]);
169 }
170 else{
171 sumA += x[i]*partCoeffPot[i]*Math.abs(charges[i]);
172 }
173 }
174 }
175 double factorA = excess/(sumA*this.volumeA);
176 double factorB = excess/((sumB+sumC)*this.volumeB);
177 for(int i=0; i<this.numOfIons; i++){
178 if(potential>0.0D){
179 if(this.charges[i]>0.0D){
180 this.excessConcnB[i] = Math.abs(this.concnB[i]*factorB);
181 this.excessComplex[i] = Math.abs(this.complex[i]*factorB);
182
183 }
184 else{
185 this.excessConcnA[i] = Math.abs(this.concnA[i]*factorA);
186 }
187 }
188 else{
189 if(this.charges[i]<0.0D){
190 this.excessConcnB[i] = Math.abs(this.concnB[i]*factorB);
191 this.excessComplex[i] = Math.abs(this.complex[i]*factorB);
192 }
193 else{
194 this.excessConcnA[i] = Math.abs(this.concnA[i]*factorA);
195 }
196 }
197 }
198 }
199 }
200
201
202 // calculates the excess charge as (surface charge density) in the interfacial region on the Patrition B side
203 public double interfaceCharge(double[] ions, double potential){
204
205 if(potential==0){
206 this.interfacialCharge=0.0D;
207 this.interfacialChargeDensity=0.0D;
208 this.diffPotentialA = 0.0D;
209 this.diffPotentialB = 0.0D;
210 this.sternPotential = 0.0D;
211 }
212 else{
213 // bisection method
214 double sigmaM = 0.0D;
215 double funcM = 0.0D;
216 double sigmaL = 0.0D, funcL = 0.0D;
217 double sumAions = 0.0D;
218 double sumBions = 0.0D;
219 double aveCharge = 0.0D;
220 for(int i=0; i<this.numOfIons; i++){
221 sumBions += Math.abs(ions[i]*this.charges[i]);
222 sumAions += Math.abs(ions[i]*this.charges[i]*this.partCoeffPot[i]);
223 aveCharge += Math.abs(charges[i]);
224 }
225 aveCharge /= this.numOfIons;
226 sumBions /= (2.0D*aveCharge);
227 sumAions /= (2.0D*aveCharge);
228 double maxQ = 1.2D*Math.sqrt(8.0D*Fmath.N_AVAGADRO*sumBions*Fmath.K_BOLTZMANN*this.temp*Fmath.EPSILON_0*this.epsilonB)*Fmath.sinh(-aveCharge*Fmath.Q_ELECTRON*Math.abs(potential)/(2.0D*Fmath.K_BOLTZMANN*this.temp));
229 double sigmaH = maxQ, funcH = 0.0D;
230 double tolQ = Math.abs(potential)*1e-8;
231 int nIterQ = 10000;
232 boolean testQ = true;
233 int iExpandQ = 0, iBisectQ = 0;
234 double diffQ = 0.0D;
235
236 while(testQ){
237 funcL = icFunct(sigmaL, potential, ions);
238 funcH = icFunct(sigmaH, potential, ions);
239 if(funcH*funcL>0.0D){
240 iExpandQ++;
241 if(iExpandQ>10)throw new IllegalArgumentException("iExpandQ has reached its limit");
242 diffQ = sigmaH - sigmaL;
243 sigmaH += diffQ;
244 }
245 else{
246 testQ=false;
247 }
248 }
249 if(Math.abs(funcL)<=tolQ){
250 sigmaM = sigmaL;
251 }
252 else{
253 if(Math.abs(funcH)<=tolQ){
254 sigmaM = sigmaH;
255 }
256 else{
257 testQ=true;
258 while(testQ){
259 sigmaM = (sigmaL + sigmaH)/2.0D;
260 funcM = icFunct(sigmaM, potential, ions);
261 if(Math.abs(funcM)<=tolQ){
262 testQ = false;
263 }
264 else{
265 if(funcL*funcM>0.0D){
266 funcL = funcM;
267 sigmaL = sigmaM;
268 }
269 else{
270 funcH = funcM;
271 sigmaH = sigmaM;
272 }
273 }
274 iBisectQ++;
275 if(iBisectQ>nIterQ){
276 System.out.println("Class: DonnanConcn\nMethod: interfaceCharge");
277 System.out.println("Maximum iterations in bisection procedure exceeded\nCurrent value of interface charge returned");
278 testQ = false;
279 }
280 }
281 }
282 }
283 this.interfacialCharge = sigmaM;
284 this.interfacialChargeDensity = sigmaM/this.interfacialArea;
285 }
286
287 // return equivalent moles of total excess charge
288 return this.interfacialCharge/(-Fmath.Q_ELECTRON*Fmath.N_AVAGADRO);
289 }
290
291 // function for interfaceCharge bisection procedure
292 // returns the present estimate of the Donnan potential minus the sum of the potential
293 // differences calculated for the Stern layer and both ionic double layers
294 public double icFunct(double sigma, double potential, double[] ions){
295 double sigmaAbs = Math.abs(sigma);
296 double sgn = Fmath.sign(potential);
297
298 if(this.chargeSame){
299 double NtotalA = 0.0D;
300 double NtotalB = 0.0D;
301 for(int i=0; i<this.numOfIons; i++){
302 NtotalA += this.concnA[i];
303 NtotalB += this.concnB[i] + this.complex[i];
304 }
305 NtotalA /= 2;
306 NtotalB /= 2;
307 double preterm = (2.0D*Fmath.K_BOLTZMANN*this.temp)/(-this.chargeValue*Fmath.Q_ELECTRON);
308 this.diffPotentialA = sgn*preterm*Fmath.asinh(sigmaAbs/Math.sqrt(8.0D*Fmath.N_AVAGADRO*NtotalA*Fmath.K_BOLTZMANN*this.temp*Fmath.EPSILON_0*this.epsilonA));
309 this.diffPotentialB = sgn*preterm*Fmath.asinh(sigmaAbs/Math.sqrt(8.0D*Fmath.N_AVAGADRO*NtotalB*Fmath.K_BOLTZMANN*this.temp*Fmath.EPSILON_0*this.epsilonB));
310 }
311 else{
312 // bisection method for double layer potentials
313 // partition A
314 double phiAm = 0.0D;
315 double funcM = 0.0D;
316 double phiAl = 0.0D;
317 double funcL = 0.0D;
318 double maxPhiA = 1.1*potential;
319 double phiAh = maxPhiA;
320 double funcH = 0.0D;
321 double tolP = Math.abs(sigma)*1e-1;
322 int nIterP = 1000;
323 boolean testP = true;
324 int iExpandP = 0, iBisectP = 0;
325 double diffP = 0.0D;
326
327 while(testP){
328 funcL = phiAfunct(sigma, phiAl, ions);
329 funcH = phiAfunct(sigma, phiAh, ions);
330 if(funcH*funcL>0.0D){
331 iExpandP++;
332 if(iExpandP>10)throw new IllegalArgumentException("iExpandP (partition A) has reached its limit");
333 diffP = phiAh - phiAl;
334 phiAh += diffP;
335 }
336 else{
337 testP=false;
338 }
339 }
340
341 testP=true;
342 while(testP){
343 phiAm = (phiAl + phiAh)/2.0D;
344 funcM = phiAfunct(sigma, phiAm, ions);
345 if(Math.abs(funcM)<=tolP){
346 this.diffPotentialA = sgn*phiAm;
347 testP = false;
348 }
349 else{
350 if(funcL*funcM>0.0D){
351 funcL = funcM;
352 phiAl = phiAm;
353 }
354 else{
355 funcH = funcM;
356 phiAh = phiAm;
357 }
358 }
359 iBisectP++;
360 if(iBisectP>nIterP){
361 //System.out.println("Class: DonnanConcn\nMethod: icFunct");
362 //System.out.println("Maximum iterations in bisection A procedure exceeded\nCurrent value of interface charge returned");
363 System.out.println("phiA = " + phiAm + " sigma = " + sigma + " funcM = " + funcM + " tol = " + tolP);
364 this.diffPotentialA = sgn*phiAm;
365 testP = false;
366 }
367 }
368
369 // partition B
370 double phiBm = 0.0D;
371 double phiBl = 0.0D;
372 double maxPhiB = -1.1*potential;
373 double phiBh = maxPhiB;
374 tolP = Math.abs(potential)*1e-5;
375 if(tolP==0.0D)tolP=1e-6;
376 nIterP = 100000;
377 testP = true;
378 iExpandP = 0;
379 iBisectP = 0;
380 diffP = 0.0D;
381
382 while(testP){
383 funcL = phiAfunct(sigma, phiBl, ions);
384 funcH = phiAfunct(sigma, phiBh, ions);
385 if(funcH*funcL>0.0D){
386 iExpandP++;
387 if(iExpandP>10)throw new IllegalArgumentException("iExpandP (partition B) has reached its limit");
388 diffP = phiBh - phiBl;
389 phiBh += diffP;
390 }
391 else{
392 testP=false;
393 }
394 }
395
396 if(Math.abs(funcH)<=tolP){
397 phiBm = phiBh;
398 }
399 else{
400 testP=true;
401 while(testP){
402 phiBm = (phiBl + phiBh)/2.0D;
403 funcM = phiAfunct(sigma, phiBm, ions);
404 if(Math.abs(funcM)<=tolP){
405 testP = false;
406 }
407 else{
408 if(funcL*funcM>0.0D){
409 funcL = funcM;
410 phiBl = phiBm;
411 }
412 else{
413 funcH = funcM;
414 phiBh = phiBm;
415 }
416 }
417 iBisectP++;
418 if(iBisectP>nIterP){
419 System.out.println("Class: DonnanConcn\nMethod: icFunct");
420 System.out.println("Maximum iterations in bisection B procedure exceeded\nCurrent value of interface charge returned");
421 System.out.println("phiB = " + phiBm + " maxPhiB = " + maxPhiB + " funcM = " + funcM + " tol = " + tolP);
422 testP = false;
423 }
424 }
425 }
426 this.diffPotentialB = sgn*phiBm;
427 }
428
429 // Calculate Stern capacitance
430 sternCapacitance(ions, sigma, this.diffPotentialA, -this.diffPotentialB);
431 this.sternPotential = sgn*sigmaAbs/this.sternCap;
432
433 return potential - (this.diffPotentialA + this.diffPotentialB + this.sternPotential);
434 }
435
436 // Function returns estimated interface charge - interface charge calculated for passed potential
437 // partition A
438 public double phiAfunct(double sigma, double potential, double[] ions){
439 double sumAsigma = 0.0D;
440 double sgns = Fmath.sign(sigma);
441 double preterm1 = 2.0D*Fmath.EPSILON_0*this.epsilonA*Fmath.K_BOLTZMANN*this.temp*Fmath.N_AVAGADRO;
442 double preterm2 = potential*Fmath.Q_ELECTRON/(Fmath.K_BOLTZMANN*this.temp);
443
444 for(int i=0; i<this.numOfIons; i++){
445
446 sumAsigma += (preterm1*ions[i]*this.partCoeffPot[i])*(Math.exp(charges[i]*preterm2) - 1.0D);
447 }
448 if(sumAsigma<0.0){
449 sgns = - sgns;
450 sumAsigma = -sumAsigma;
451 }
452 double diffSigma = sigma - sgns*Math.sqrt(sumAsigma);
453 return diffSigma;
454 }
455
456 // Function returns estimated interface charge - interface charge calculated for passed potential
457 // partition B
458 public double phiBfunct(double sigma, double potential, double[] ions){
459 double sumBsigma = 0.0D;
460 double sgns = Fmath.sign(sigma);
461 double preterm1 = 2.0D*Fmath.EPSILON_0*this.epsilonB*Fmath.K_BOLTZMANN*this.temp*Fmath.N_AVAGADRO;
462 double preterm2 = potential*Fmath.Q_ELECTRON/(Fmath.K_BOLTZMANN*this.temp);
463
464 for(int i=0; i<this.numOfIons; i++){
465 sumBsigma += (preterm1*(ions[i]+this.complex[i]))*(Math.exp(charges[i]*preterm2) - 1.0D);
466 }
467 if(sumBsigma<0.0){
468 sgns = - sgns;
469 sumBsigma = -sumBsigma;
470 }
471 double diffSigma = sigma - sgns*Math.sqrt(sumBsigma);
472 return diffSigma;
473 }
474
475 // calculates the Stern capacitances and the Stern potential
476 // given a surface charge and given diffuse layer potentials
477 public void sternCapacitance(double[] ions, double sigma, double psiA, double psiB){
478
479 double denomA = 0.0D;
480 double denomB = 0.0D;
481 this.sternDeltaA = 0.0D;
482 this.sternDeltaB = 0.0D;
483 double preterm = -Fmath.Q_ELECTRON/(Fmath.K_BOLTZMANN*this.temp);
484 for(int i=0; i<this.numOfIons; i++){
485 this.sternDeltaA += radii[i]*ions[i]*this.partCoeffPot[i]*Math.exp(preterm*psiA*charges[i]);
486 this.sternDeltaB += (this.radii[i]*ions[i] + this.ionophoreRad*this.complex[i])*Math.exp(-preterm*psiB*charges[i]);
487 denomA += ions[i]*this.partCoeffPot[i]*Math.exp(preterm*psiA*charges[i]);
488 denomB += (ions[i] + this.complex[i])*Math.exp(-preterm*psiB*charges[i]);
489 }
490 this.sternDeltaA /= denomA;
491 this.sternDeltaB /= denomB;
492 this.sternCap = Fmath.EPSILON_0*this.epsilonSternA*this.epsilonSternB/(this.sternDeltaA*this.epsilonSternB + this.sternDeltaB*this.epsilonSternA);
493 }
494}
495
496
Note: See TracBrowser for help on using the repository browser.