ISLEC  Version 4.2
psc_model.c
Go to the documentation of this file.
1 /* src/psc_model.c
2  *
3  * Copyright (C) 2011-2018 Dongdong Li
4  *
5  * This program is free software; you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundataion; either version 3 of the License, or (at
8  * your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful, but
11  * WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABLITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13  * General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program; if not, write to the Free Software
17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18  */
19 
32 #define CONSTANT_MW 18.015 // molar mass of water molecular
33 #define CONSTANT_pi 3.1415926 // value of constant pai
34 #define CONSTANT_NA 6.02214129E+23 // value of constant Avogadro constant
35 #define CONSTANT_k 1.3806448E-16 // unit in esu
36 #define CONSTANT_e 4.8024E-10 // unit in esu
37 #define CONSTANT_pc 0.322 // critical pressure
38 #define CONSTANT_Tc 647.096 // critical temperature
39 
40 #define CONSTANT_b1 1.99274064
41 #define CONSTANT_b2 1.09965342
42 #define CONSTANT_b3 -0.510839303
43 #define CONSTANT_b4 -1.75493479
44 #define CONSTANT_b5 -45.5170352
45 #define CONSTANT_b6 -6.74694450E+5
46 
47 #define tao_T (1.0 - T/CONSTANT_Tc)
48 #define p_T (CONSTANT_pc)*(1 + CONSTANT_b1 * pow (tao_T, 1.0/3.0) + CONSTANT_b2 * pow (tao_T, 2.0/3.0) \
49  + CONSTANT_b3 * pow (tao_T, 5.0/3.0) + CONSTANT_b4 * pow (tao_T, 16.0/3.0) \
50  + CONSTANT_b5 * pow (tao_T, 43.0/3.0) + CONSTANT_b6 * pow (tao_T, 110.0/3.0)) // g.cm3
51 
52 #define CONSTANT_U1 3.4279E+2
53 #define CONSTANT_U2 -5.0866E-3
54 #define CONSTANT_U3 9.4690E-7
55 #define CONSTANT_U4 -2.0525
56 #define CONSTANT_U5 3.1159E+3
57 #define CONSTANT_U6 -1.8289E+2
58 #define CONSTANT_U7 -8.0325E+3
59 #define CONSTANT_U8 4.2142E+6
60 #define CONSTANT_U9 2.1417
61 
62 #define D1000_T (CONSTANT_U1 * exp (CONSTANT_U2 * T + CONSTANT_U3 * T * T))
63 #define C_T (CONSTANT_U4 + CONSTANT_U5/(CONSTANT_U6 + T))
64 #define B_T (CONSTANT_U7 + CONSTANT_U8/T + CONSTANT_U9 * T)
65 #define D_T (D1000_T + C_T * log ((B_T + P)/(B_T + 1000)))
66 
67 #define dp_T -(CONSTANT_pc/CONSTANT_Tc)*(CONSTANT_b1/3.0 * pow (tao_T, -2.0/3.0) \
68  + 2.0 * CONSTANT_b2/3.0 * pow (tao_T, -1.0/3.0) \
69  + 5.0 * CONSTANT_b3/3.0 * pow (tao_T, 2.0/3.0) \
70  + 16.0 * CONSTANT_b4/3.0 * pow (tao_T, 13/3.0) \
71  + 43.0 * CONSTANT_b5/3.0 * pow (tao_T, 40.0/3.0) \
72  + 110.0 * CONSTANT_b6/3.0 * pow (tao_T, 107/3.0))
73 
74 #define aw_T -(1/(p_T))*dp_T
75 #define dD1000_T CONSTANT_U1*exp(CONSTANT_U2*T + CONSTANT_U3*T*T)*(CONSTANT_U2 + 2*CONSTANT_U3*T)
76 #define dC_T -CONSTANT_U5/pow(CONSTANT_U6 + T, 2)
77 #define dB_T (-CONSTANT_U8/(T * T) + CONSTANT_U9)
78 #define dBP1000_T ((B_T + 1000)*dB_T - (B_T + P)*dB_T)/pow(B_T + 1000, 2)
79 #define dB1000P_T ((B_T + P)*dB_T - (B_T + 1000)*dB_T)/pow(B_T + P, 2)
80 #define dD_T (dD1000_T + dC_T * log ((B_T + P)/(B_T + 1000)) + C_T * (B_T + 1000)/(B_T + P)*dBP1000_T)
81 #define dlnD_T (1/D_T)*dD_T
82 #define d2D1000_T CONSTANT_U1*exp(CONSTANT_U2*T + CONSTANT_U3*T*T)*(pow(CONSTANT_U2 + 2*CONSTANT_U3*T, 2) + 2*CONSTANT_U3)
83 #define d2C_T 2*CONSTANT_U5/pow(CONSTANT_U6 + T, 3)
84 #define d2B_T 2*CONSTANT_U8/pow(T, 3)
85 #define d2BP1000_T (((B_T + 1000) * d2B_T - d2B_T * (B_T + P)) * (B_T + 1000) * (B_T + 1000) \
86  - 2 * (B_T + 1000) * dB_T * \
87  ((B_T + 1000) * dB_T - (B_T + P) * dB_T))/pow(B_T + 1000, 4)
88 #define d2D_T (d2D1000_T + d2C_T * log ((B_T + P)/(B_T + 1000)) \
89  + dC_T*(B_T+1000)/(B_T+P)*dBP1000_T + (dC_T * (B_T + 1000)/(B_T + P) \
90  + C_T * dB1000P_T) * dBP1000_T + C_T * (B_T+1000)/(B_T + P) * d2BP1000_T)
91 #define d2p_T (CONSTANT_pc/(CONSTANT_Tc * CONSTANT_Tc)) * ((-2.0*CONSTANT_b1/9.0)* pow (tao_T, -5.0/3.0) + \
92  (-2.0*CONSTANT_b2/9.0)* pow (tao_T, -4.0/3.0) + \
93  (10.0*CONSTANT_b3/9.0)* pow (tao_T, -1.0/3.0) + \
94  (208.0*CONSTANT_b4/9.0)* pow (tao_T, 10.0/3.0) + \
95  (1720.0*CONSTANT_b5/9.0)* pow (tao_T, 37.0/3.0) + \
96  (11770.0*CONSTANT_b6/9.0)* pow (tao_T, 104.0/3.0))
97 
98 #define ROU 13.0 // here used a constant p = 13.0 instead of (2150*sqrt(p_T/D_T/T))
99 
100 
101 static double aphi (double T, double P)
102 {
103  return (1.0/3.0) * sqrt(2.0 * CONSTANT_pi * CONSTANT_NA * p_T / 1000.0) * pow(CONSTANT_e*CONSTANT_e/(D_T*CONSTANT_k*T), 1.5);
104 }
105 
106 static double ax (double T, double P)
107 {
108  return sqrt(1000/CONSTANT_MW) * aphi (T, P);
109 }
110 
111 static double ah (double T, double P)
112 {
113  return -6.0 * R * T * aphi (T, P) * (1.0 + T * dlnD_T + T * aw_T/3.0);
114 }
115 
116 static double ahx (double T, double P)
117 {
118  return sqrt(1000/CONSTANT_MW) * ah(T, P);
119 }
120 
121 static double aj (double T, double P)
122 {
123  return 3 * aphi (T, P) * R * T * (1.0/T + 2 * dlnD_T
124  + 5 * T * dlnD_T * dlnD_T + 2 * T * dlnD_T * aw_T + 2.0/3.0 * aw_T
125  + T * aw_T * aw_T - (2 * T / D_T) * d2D_T - (2.0 * T * p_T / 3.0) *
126  (2.0 / pow(p_T, 3) * dp_T * dp_T - (1.0/pow(p_T, 2)) * d2p_T));
127 }
128 
129 static double ajx (double T, double P)
130 {
131  return sqrt(1000/CONSTANT_MW) * aj(T, P);
132 }
133 
134 /* g(x) */
135 static double g(double x)
136 {
137  if (x != 0)
138  {
139  return 2 * (1 - (1 + x) * exp(-x) ) / pow (x, 2);
140  }
141  else
142  {
143  return 0;
144  }
145 }
146 
147 /* g'(x) */
148 static double gp(double x)
149 {
150  if (x != 0)
151  {
152  return -2 * (1 - (1 + x + pow (x, 2) / 2.0) * exp (-x)) / pow (x, 2);
153  }
154  else
155  {
156  return 0;
157  }
158 }
159 
160 static double jj (double x)
161 {
162  return x * pow ( (4.0 + C1 * pow (x, -C2) * exp (-C3 * pow (x, C4))), -1.0);
163 }
164 
165 static double jp (double x)
166 {
167  return pow ( (4.0 + C1 * pow (x, -C2) * exp (-C3 * pow (x, C4))), -1.0) +
168  pow (4.0 + C1 * pow (x, -C2) * exp (-C3 * pow (x, C4)), -2.0) *
169  ( (C1 * x * exp (-C3 * pow (x, C4))) * (C2 * pow (x, -C2-1.0) +
170  C3 * C4 * pow (x, C4-1.0) * pow (x, -C2)));
171 }
172 
173 static double psc_rn (AQUEOUS_PHASE aq, int index, double *n, double T, double P)
174 {
175  if (aq.aqueous_species[index].charge != 0)
176  {
177  printf ("The species is not a neutral species!!!\n");
178  return 1.0;
179  }
180 
181  if (index >= aq.species_num)
182  {
183  printf ("The species dose not existed!!!\n");
184  return 1.0;
185  }
186 
187  double lnfn;
188  double x[aq.species_num];
189  double Z[aq.species_num];
190  double E[aq.species_num];
191  int i, j, k, l, m, s;
192 
193  double sumn;
194  double Ax, Ix, F;
195 
196  double aii[aq.species_num][aq.species_num];
197  double a1ii[aq.species_num][aq.species_num];
198  double Bii[aq.species_num][aq.species_num];
199  double B1ii[aq.species_num][aq.species_num];
200  double Wnn[aq.species_num][aq.species_num];
201  double Unn[aq.species_num][aq.species_num];
202  double Wnii[aq.species_num][aq.species_num][aq.species_num];
203  double Unii[aq.species_num][aq.species_num][aq.species_num];
204  double Vnii[aq.species_num][aq.species_num][aq.species_num];
205  double Wiii[aq.species_num][aq.species_num][aq.species_num];
206  double Qniii[aq.species_num][aq.species_num][aq.species_num][aq.species_num];
207  double Yniii[aq.species_num][aq.species_num][aq.species_num][aq.species_num];
208  double Ynnii[aq.species_num][aq.species_num][aq.species_num][aq.species_num];
209 
210  double Bii_a, Bii_b, Bii_c, Bii_d, Bii_e, Bii_f, \
211  B1ii_a, B1ii_b, B1ii_c, B1ii_d, B1ii_e, B1ii_f, \
212  Wnn_a, Wnn_b, Wnn_c, Wnn_d, Wnn_e, Wnn_f, \
213  Unn_a, Unn_b, Unn_c, Unn_d, Unn_e, Unn_f, \
214  Wnii_a, Wnii_b, Wnii_c, Wnii_d, Wnii_e, Wnii_f, \
215  Unii_a, Unii_b, Unii_c, Unii_d, Unii_e, Unii_f, \
216  Vnii_a, Vnii_b, Vnii_c, Vnii_d, Vnii_e, Vnii_f, \
217  Wiii_a, Wiii_b, Wiii_c, Wiii_d, Wiii_e, Wiii_f, \
218  Qniii_a, Qniii_b, Qniii_c, Qniii_d, Qniii_e, Qniii_f, \
219  Yniii_a, Yniii_b, Yniii_c, Yniii_d, Yniii_e, Yniii_f, \
220  Ynnii_a, Ynnii_b, Ynnii_c, Ynnii_d, Ynnii_e, Ynnii_f;
221 
222  Ax = ax (T, 1.0);
223 
224  for (k = 0; k < aq.species_num; k ++)
225  {
226  Z[k] = aq.aqueous_species[k].charge;
227  }
228 
229  sumn = 0;
230  for (k = 0; k < aq.species_num; k ++)
231  {
232  sumn += n[k];
233  }
234 
235  for (k = 0; k < aq.species_num; k ++)
236  {
237  x[k] = n[k] / sumn;
238  }
239 
240  F = 0;
241  for (k = 0; k < aq.species_num; k++)
242  {
243  F += 0.5 * x[k] * fabs(Z[k]);
244  }
245  F = 1.0/F;
246 
247  Ix = 0;
248  for (k = 0; k < aq.species_num; k ++)
249  {
250  Ix += 0.5 * (x[k] * Z[k] * Z[k]);
251  }
252 
253  for (k = 0; k < aq.species_num; k ++)
254  {
255  E[k] = 0;
256  }
257 
258  // INITALIZATION set all binary parameters zero
259  for (i = 0; i < aq.species_num; i ++)
260  {
261  for (j = 0; j < aq.species_num; j ++)
262  {
263  aii[i][j] = 0;
264  a1ii[i][j] = 0;
265  Bii[i][j] = 0;
266  B1ii[i][j] = 0;
267  Wnn[i][j] = 0;
268  Unn[i][j] = 0;
269  }
270  }
271 
272  // INITALIZATION set all ternary parameters zero
273  for (i = 0; i < aq.species_num; i ++)
274  {
275  for (j = 0; j < aq.species_num; j ++)
276  {
277  for (l = 0; l < aq.species_num; l ++)
278  {
279  Wnii[i][j][l] = 0;
280  Unii[i][j][l] = 0;
281  Vnii[i][j][l] = 0;
282  Wiii[i][j][l] = 0;
283  }
284  }
285  }
286 
287  // INITALIZATION set all quaternary parameters zero
288  for (i = 0; i < aq.species_num; i ++)
289  {
290  for (j = 0; j < aq.species_num; j ++)
291  {
292  for (l = 0; l < aq.species_num; l ++)
293  {
294  for (m = 0; m < aq.species_num; m ++)
295  {
296  Qniii[i][j][l][m] = 0;
297  Yniii[i][j][l][m] = 0;
298  Ynnii[i][j][l][m] = 0;
299  }
300  }
301  }
302  }
303 
304  for (k = 0; k < aq.ii_bin_param_num; k ++)
305  {
306  for (m = 0; m < aq.species_num; m ++)
307  {
308  if (strcmp (aq.ii_bin_params[k].species_s1, aq.aqueous_species[m].species_symbol) == 0)
309  {
310  i = m;
311  }
312  }
313 
314  for (m = 0; m < aq.species_num; m ++)
315  {
316  if (strcmp (aq.ii_bin_params[k].species_s2, aq.aqueous_species[m].species_symbol) == 0)
317  {
318  j = m;
319  }
320  }
321 
322  double aii_value = aq.ii_bin_params[k].alpha;
323  double a1ii_value = aq.ii_bin_params[k].alpha1;
324 
325  Bii_a = aq.ii_bin_params[k].B0ss[0];
326  Bii_b = aq.ii_bin_params[k].B0ss[1];
327  Bii_c = aq.ii_bin_params[k].B0ss[2];
328  Bii_d = aq.ii_bin_params[k].B0ss[3];
329  Bii_e = aq.ii_bin_params[k].B0ss[4];
330  Bii_f = aq.ii_bin_params[k].B0ss[5];
331  B1ii_a = aq.ii_bin_params[k].B1ss[0];
332  B1ii_b = aq.ii_bin_params[k].B1ss[1];
333  B1ii_c = aq.ii_bin_params[k].B1ss[2];
334  B1ii_d = aq.ii_bin_params[k].B1ss[3];
335  B1ii_e = aq.ii_bin_params[k].B1ss[4];
336  B1ii_f = aq.ii_bin_params[k].B1ss[5];
337 
338  double Bii_value = Bii_a + Bii_b * T + Bii_c * T * log(T) + Bii_d * T * T + Bii_e * T * T * T + Bii_f/T;
339  double B1ii_value = B1ii_a + B1ii_b * T + B1ii_c * T * log(T) + B1ii_d * T * T + B1ii_e * T * T * T + B1ii_f/T;
340 
341  aii[i][j] = aii_value; aii[j][i] = aii[i][j];
342  a1ii[i][j] = a1ii_value; a1ii[j][i] = a1ii[i][j];
343  Bii[i][j] = Bii_value; Bii[j][i] = Bii[i][j];
344  B1ii[i][j] = B1ii_value; B1ii[j][i] = B1ii[i][j];
345  }
346 
347  for (k = 0; k < aq.nn_bin_param_num; k ++)
348  {
349  for (m = 0; m < aq.species_num; m ++)
350  {
351  if (strcmp (aq.nn_bin_params[k].species_s1, aq.aqueous_species[m].species_symbol) == 0)
352  {
353  i = m;
354  }
355  }
356 
357  for (m = 0; m < aq.species_num; m ++)
358  {
359  if (strcmp (aq.nn_bin_params[k].species_s2, aq.aqueous_species[m].species_symbol) == 0)
360  {
361  j = m;
362  }
363  }
364 
365  Wnn_a = aq.nn_bin_params[k].Wss[0];
366  Wnn_b = aq.nn_bin_params[k].Wss[1];
367  Wnn_c = aq.nn_bin_params[k].Wss[2];
368  Wnn_d = aq.nn_bin_params[k].Wss[3];
369  Wnn_e = aq.nn_bin_params[k].Wss[4];
370  Wnn_f = aq.nn_bin_params[k].Wss[5];
371  Unn_a = aq.nn_bin_params[k].Uss[0];
372  Unn_b = aq.nn_bin_params[k].Uss[1];
373  Unn_c = aq.nn_bin_params[k].Uss[2];
374  Unn_d = aq.nn_bin_params[k].Uss[3];
375  Unn_e = aq.nn_bin_params[k].Uss[4];
376  Unn_f = aq.nn_bin_params[k].Uss[5];
377 
378  double Wnn_value = Wnn_a + Wnn_b * T + Wnn_c * T * log(T) + Wnn_d * T * T + Wnn_e * T * T * T + Wnn_f/T;
379  double Unn_value = Unn_a + Unn_b * T + Unn_c * T * log(T) + Unn_d * T * T + Unn_e * T * T * T + Unn_f/T;
380 
381  Wnn[i][j] = Wnn_value; Wnn[j][i] = Wnn[i][j];
382  Unn[i][j] = Unn_value; Unn[j][i] = Unn[i][j];
383  }
384 
385  for (k = 0; k < aq.nii_ter_param_num; k++)
386  {
387  for (m = 0; m < aq.species_num; m ++)
388  {
389  if (strcmp (aq.nii_ter_params[k].species_s1, aq.aqueous_species[m].species_symbol) == 0)
390  {
391  i = m;
392  }
393  }
394 
395  for (m = 0; m < aq.species_num; m ++)
396  {
397  if (strcmp (aq.nii_ter_params[k].species_s2, aq.aqueous_species[m].species_symbol) == 0)
398  {
399  j = m;
400  }
401  }
402 
403  for (m = 0; m < aq.species_num; m ++)
404  {
405  if (strcmp (aq.nii_ter_params[k].species_s3, aq.aqueous_species[m].species_symbol) == 0)
406  {
407  l = m;
408  }
409  }
410 
411  Wnii_a = aq.nii_ter_params[k].Wsss[0];
412  Wnii_b = aq.nii_ter_params[k].Wsss[1];
413  Wnii_c = aq.nii_ter_params[k].Wsss[2];
414  Wnii_d = aq.nii_ter_params[k].Wsss[3];
415  Wnii_e = aq.nii_ter_params[k].Wsss[4];
416  Wnii_f = aq.nii_ter_params[k].Wsss[5];
417  Unii_a = aq.nii_ter_params[k].Usss[0];
418  Unii_b = aq.nii_ter_params[k].Usss[1];
419  Unii_c = aq.nii_ter_params[k].Usss[2];
420  Unii_d = aq.nii_ter_params[k].Usss[3];
421  Unii_e = aq.nii_ter_params[k].Usss[4];
422  Unii_f = aq.nii_ter_params[k].Usss[5];
423  Vnii_a = aq.nii_ter_params[k].Vsss[0];
424  Vnii_b = aq.nii_ter_params[k].Vsss[1];
425  Vnii_c = aq.nii_ter_params[k].Vsss[2];
426  Vnii_d = aq.nii_ter_params[k].Vsss[3];
427  Vnii_e = aq.nii_ter_params[k].Vsss[4];
428  Vnii_f = aq.nii_ter_params[k].Vsss[5];
429 
430  double Wnii_value = Wnii_a + Wnii_b * T + Wnii_c * T * log(T) + Wnii_d * T * T + Wnii_e * T * T * T + Wnii_f/T;
431  double Unii_value = Unii_a + Unii_b * T + Unii_c * T * log(T) + Unii_d * T * T + Unii_e * T * T * T + Unii_f/T;
432  double Vnii_value = Vnii_a + Vnii_b * T + Vnii_c * T * log(T) + Vnii_d * T * T + Vnii_e * T * T * T + Vnii_f/T;
433 
434  Wnii[i][j][l] = Wnii_value;
435  Wnii[i][l][j] = Wnii[i][j][l];
436  Wnii[j][i][l] = Wnii[i][j][l];
437  Wnii[j][l][i] = Wnii[i][j][l];
438  Wnii[l][i][j] = Wnii[i][j][l];
439  Wnii[l][j][i] = Wnii[i][j][l];
440 
441  Unii[i][j][l] = Unii_value;
442  Unii[i][l][j] = Unii[i][j][l];
443  Unii[j][i][l] = Unii[i][j][l];
444  Unii[j][l][i] = Unii[i][j][l];
445  Unii[l][i][j] = Unii[i][j][l];
446  Unii[l][j][i] = Unii[i][j][l];
447 
448  Vnii[i][j][l] = Vnii_value;
449  Vnii[i][l][j] = Vnii[i][j][l];
450  Vnii[j][i][l] = Vnii[i][j][l];
451  Vnii[j][l][i] = Vnii[i][j][l];
452  Vnii[l][i][j] = Vnii[i][j][l];
453  Vnii[l][j][i] = Vnii[i][j][l];
454  }
455 
456  for (k = 0; k < aq.iii_ter_param_num; k++)
457  {
458  for (m = 0; m < aq.species_num; m ++)
459  {
460  if (strcmp (aq.iii_ter_params[k].species_s1, aq.aqueous_species[m].species_symbol) == 0)
461  {
462  i = m;
463  }
464  }
465 
466  for (m = 0; m < aq.species_num; m ++)
467  {
468  if (strcmp (aq.iii_ter_params[k].species_s2, aq.aqueous_species[m].species_symbol) == 0)
469  {
470  j = m;
471  }
472  }
473 
474  for (m = 0; m < aq.species_num; m ++)
475  {
476  if (strcmp (aq.iii_ter_params[k].species_s3, aq.aqueous_species[m].species_symbol) == 0)
477  {
478  l = m;
479  }
480  }
481 
482  Wiii_a = aq.iii_ter_params[k].Wsss[0];
483  Wiii_b = aq.iii_ter_params[k].Wsss[1];
484  Wiii_c = aq.iii_ter_params[k].Wsss[2];
485  Wiii_d = aq.iii_ter_params[k].Wsss[3];
486  Wiii_e = aq.iii_ter_params[k].Wsss[4];
487  Wiii_f = aq.iii_ter_params[k].Wsss[5];
488 
489  double Wiii_value = Wiii_a + Wiii_b * T + Wiii_c * T * log(T) + Wiii_d * T * T + Wiii_e * T * T * T + Wiii_f/T;
490 
491  Wiii[i][j][l] = Wiii_value;
492  Wiii[i][l][j] = Wiii[i][j][l];
493  Wiii[j][i][l] = Wiii[i][j][l];
494  Wiii[j][l][i] = Wiii[i][j][l];
495  Wiii[l][i][j] = Wiii[i][j][l];
496  Wiii[l][j][i] = Wiii[i][j][l];
497  }
498 
499  for (k = 0; k < aq.niii_qua_param_num; k++)
500  {
501  for (m = 0; m < aq.species_num; m ++)
502  {
503  if (strcmp (aq.niii_qua_params[k].species_s1, aq.aqueous_species[m].species_symbol) == 0)
504  {
505  i = m;
506  }
507  }
508 
509  for (m = 0; m < aq.species_num; m ++)
510  {
511  if (strcmp (aq.niii_qua_params[k].species_s2, aq.aqueous_species[m].species_symbol) == 0)
512  {
513  j = m;
514  }
515  }
516 
517  for (m = 0; m < aq.species_num; m ++)
518  {
519  if (strcmp (aq.niii_qua_params[k].species_s3, aq.aqueous_species[m].species_symbol) == 0)
520  {
521  l = m;
522  }
523  }
524 
525  for (m = 0; m < aq.species_num; m ++)
526  {
527  if (strcmp (aq.niii_qua_params[k].species_s4, aq.aqueous_species[m].species_symbol) == 0)
528  {
529  s = m;
530  }
531  }
532  Qniii_a = aq.niii_qua_params[k].Qssss[0];
533  Qniii_b = aq.niii_qua_params[k].Qssss[1];
534  Qniii_c = aq.niii_qua_params[k].Qssss[2];
535  Qniii_d = aq.niii_qua_params[k].Qssss[3];
536  Qniii_e = aq.niii_qua_params[k].Qssss[4];
537  Qniii_f = aq.niii_qua_params[k].Qssss[5];
538  Yniii_a = aq.niii_qua_params[k].Yssss[0];
539  Yniii_b = aq.niii_qua_params[k].Yssss[1];
540  Yniii_c = aq.niii_qua_params[k].Yssss[2];
541  Yniii_d = aq.niii_qua_params[k].Yssss[3];
542  Yniii_e = aq.niii_qua_params[k].Yssss[4];
543  Yniii_f = aq.niii_qua_params[k].Yssss[5];
544 
545  double Qniii_value = Qniii_a + Qniii_b * T + Qniii_c * T * log(T) + Qniii_d * T * T + Qniii_e * T * T * T + Qniii_f/T;
546  double Yniii_value = Yniii_a + Yniii_b * T + Yniii_c * T * log(T) + Yniii_d * T * T + Yniii_e * T * T * T + Yniii_f/T;
547 
548  Qniii[i][j][l][s] = Qniii_value;
549  Qniii[i][j][s][l] = Qniii[i][j][l][s];
550  Qniii[i][l][j][s] = Qniii[i][j][l][s];
551  Qniii[i][l][s][j] = Qniii[i][j][l][s];
552  Qniii[i][s][j][l] = Qniii[i][j][l][s];
553  Qniii[i][s][l][j] = Qniii[i][j][l][s];
554 
555  Qniii[j][i][l][s] = Qniii[i][j][l][s];
556  Qniii[j][i][s][l] = Qniii[i][j][l][s];
557  Qniii[j][l][i][s] = Qniii[i][j][l][s];
558  Qniii[j][l][s][i] = Qniii[i][j][l][s];
559  Qniii[j][s][i][l] = Qniii[i][j][l][s];
560  Qniii[j][s][l][i] = Qniii[i][j][l][s];
561 
562  Qniii[l][i][j][s] = Qniii[i][j][l][s];
563  Qniii[l][i][s][j] = Qniii[i][j][l][s];
564  Qniii[l][j][i][s] = Qniii[i][j][l][s];
565  Qniii[l][j][s][i] = Qniii[i][j][l][s];
566  Qniii[l][s][i][j] = Qniii[i][j][l][s];
567  Qniii[l][s][j][i] = Qniii[i][j][l][s];
568 
569  Qniii[s][i][j][l] = Qniii[i][j][l][s];
570  Qniii[s][i][l][j] = Qniii[i][j][l][s];
571  Qniii[s][j][l][i] = Qniii[i][j][l][s];
572  Qniii[s][j][i][l] = Qniii[i][j][l][s];
573  Qniii[s][l][i][j] = Qniii[i][j][l][s];
574  Qniii[s][l][j][i] = Qniii[i][j][l][s];
575 
576  Yniii[i][j][l][s] = Yniii_value;
577  Yniii[i][j][s][l] = Yniii[i][j][l][s];
578  Yniii[i][l][j][s] = Yniii[i][j][l][s];
579  Yniii[i][l][s][j] = Yniii[i][j][l][s];
580  Yniii[i][s][j][l] = Yniii[i][j][l][s];
581  Yniii[i][s][l][j] = Yniii[i][j][l][s];
582 
583  Yniii[j][i][l][s] = Yniii[i][j][l][s];
584  Yniii[j][i][s][l] = Yniii[i][j][l][s];
585  Yniii[j][l][i][s] = Yniii[i][j][l][s];
586  Yniii[j][l][s][i] = Yniii[i][j][l][s];
587  Yniii[j][s][i][l] = Yniii[i][j][l][s];
588  Yniii[j][s][l][i] = Yniii[i][j][l][s];
589 
590  Yniii[l][i][j][s] = Yniii[i][j][l][s];
591  Yniii[l][i][s][j] = Yniii[i][j][l][s];
592  Yniii[l][j][i][s] = Yniii[i][j][l][s];
593  Yniii[l][j][s][i] = Yniii[i][j][l][s];
594  Yniii[l][s][i][j] = Yniii[i][j][l][s];
595  Yniii[l][s][j][i] = Yniii[i][j][l][s];
596 
597  Yniii[s][i][j][l] = Yniii[i][j][l][s];
598  Yniii[s][i][l][j] = Yniii[i][j][l][s];
599  Yniii[s][j][l][i] = Yniii[i][j][l][s];
600  Yniii[s][j][i][l] = Yniii[i][j][l][s];
601  Yniii[s][l][i][j] = Yniii[i][j][l][s];
602  Yniii[s][l][j][i] = Yniii[i][j][l][s];
603  }
604 
605  for (k = 0; k < aq.nnii_qua_param_num; k++)
606  {
607  for (m = 0; m < aq.species_num; m ++)
608  {
609  if (strcmp (aq.nnii_qua_params[k].species_s1, aq.aqueous_species[m].species_symbol) == 0)
610  {
611  i = m;
612  }
613  }
614 
615  for (m = 0; m < aq.species_num; m ++)
616  {
617  if (strcmp (aq.nnii_qua_params[k].species_s2, aq.aqueous_species[m].species_symbol) == 0)
618  {
619  j = m;
620  }
621  }
622 
623  for (m = 0; m < aq.species_num; m ++)
624  {
625  if (strcmp (aq.nnii_qua_params[k].species_s3, aq.aqueous_species[m].species_symbol) == 0)
626  {
627  l = m;
628  }
629  }
630 
631  for (m = 0; m < aq.species_num; m ++)
632  {
633  if (strcmp (aq.nnii_qua_params[k].species_s4, aq.aqueous_species[m].species_symbol) == 0)
634  {
635  s = m;
636  }
637  }
638 
639  Ynnii_a = aq.nnii_qua_params[k].Yssss[0];
640  Ynnii_b = aq.nnii_qua_params[k].Yssss[1];
641  Ynnii_c = aq.nnii_qua_params[k].Yssss[2];
642  Ynnii_d = aq.nnii_qua_params[k].Yssss[3];
643  Ynnii_e = aq.nnii_qua_params[k].Yssss[4];
644  Ynnii_f = aq.nnii_qua_params[k].Yssss[5];
645 
646  double Ynnii_value = Ynnii_a + Ynnii_b * T + Ynnii_c * T * log(T) + Ynnii_d * T * T + Ynnii_e * T * T * T + Ynnii_f/T;
647 
648  Ynnii[i][j][l][s] = Ynnii_value;
649  Ynnii[i][j][s][l] = Ynnii[i][j][l][s];
650  Ynnii[i][l][j][s] = Ynnii[i][j][l][s];
651  Ynnii[i][l][s][j] = Ynnii[i][j][l][s];
652  Ynnii[i][s][j][l] = Ynnii[i][j][l][s];
653  Ynnii[i][s][l][j] = Ynnii[i][j][l][s];
654 
655  Ynnii[j][i][l][s] = Ynnii[i][j][l][s];
656  Ynnii[j][i][s][l] = Ynnii[i][j][l][s];
657  Ynnii[j][l][i][s] = Ynnii[i][j][l][s];
658  Ynnii[j][l][s][i] = Ynnii[i][j][l][s];
659  Ynnii[j][s][i][l] = Ynnii[i][j][l][s];
660  Ynnii[j][s][l][i] = Ynnii[i][j][l][s];
661 
662  Ynnii[l][i][j][s] = Ynnii[i][j][l][s];
663  Ynnii[l][i][s][j] = Ynnii[i][j][l][s];
664  Ynnii[l][j][i][s] = Ynnii[i][j][l][s];
665  Ynnii[l][j][s][i] = Ynnii[i][j][l][s];
666  Ynnii[l][s][i][j] = Ynnii[i][j][l][s];
667  Ynnii[l][s][j][i] = Ynnii[i][j][l][s];
668 
669  Ynnii[s][i][j][l] = Ynnii[i][j][l][s];
670  Ynnii[s][i][l][j] = Ynnii[i][j][l][s];
671  Ynnii[s][j][l][i] = Ynnii[i][j][l][s];
672  Ynnii[s][j][i][l] = Ynnii[i][j][l][s];
673  Ynnii[s][l][i][j] = Ynnii[i][j][l][s];
674  Ynnii[s][l][j][i] = Ynnii[i][j][l][s];
675  }
676 
677  for (i = 0; i < aq.species_num; i ++)
678  {
679  E[i] = x[i] * fabs (Z[i]);
680  double sum = 0;
681 
682  for (k = 0; k < aq.species_num; k++)
683  {
684  if (Z[i] * Z[k] > 0) {sum = sum + x[k] * fabs (Z[k]);}
685  }
686 
687  if (sum != 0) {E[i] = E[i] / sum;}
688  }
689 
690  // DH term
691  lnfn = 2.0 * Ax * Ix * sqrt (Ix) / (1.0 + ROU * sqrt (Ix));
692 
693 
694  // cation-anion binary interaction
695  for (i = 0; i < aq.species_num; i ++)
696  {
697  for (j = 0; j < aq.species_num; j ++)
698  {
699  if (Z[i]*Z[j] != 0)
700  {
701  lnfn = lnfn - 0.5 * (x[i] * x[j] * Bii[i][j] * exp (-aii[i][j] * sqrt (Ix)) + x[i] * x[j] * B1ii[i][j] * exp (-a1ii[i][j] * sqrt (Ix)));
702  }
703  }
704  }
705 
706  // cation-cation or anion-anion ternary HOE term
707  for (i = 0; i < aq.species_num; i ++)
708  {
709  for (j = 0; j < aq.species_num; j ++)
710  {
711  if (Z[i] * Z[j] > 0 && i != j && Ix != 0)
712  {
713  double xij = 6.0 * Z[i] * Z[j] * Ax * sqrt (Ix);
714  double xii = 6.0 * Z[i] * Z[i] * Ax * sqrt (Ix);
715  double xjj = 6.0 * Z[j] * Z[j] * Ax * sqrt (Ix);
716 
717  double vij = Z[i] * Z[j] / (4.0 * Ix) * (jj(xij) - 0.5 * jj(xii) - 0.5 * jj(xjj));
718  double vpij = - vij/Ix + (Z[i]*Z[j]/(8.0 * Ix * Ix)) * (xij * jp(xij) - 0.5 * xii * jp(xii) - 0.5 * xjj * jp(xjj));
719 
720  lnfn = lnfn + 0.5 * (-2.0 * x[i] * x[j] * (vij + Ix * vpij));
721  }
722  }
723  }
724 
725  // neutral-caion-anion ternary interaction
726  for (i = 0; i < aq.species_num; i ++)
727  {
728  for (j = 0; j < aq.species_num; j ++)
729  {
730  double sum;
731 
732  if (Z[i]*Z[j] < 0)
733  {
734  double Wnii_sum = 0;
735  double Unii_sum = 0;
736  double Vnii_sum = 0;
737 
738  for (k = 0; k < aq.species_num; k ++)
739  {
740  if (aq.aqueous_species[k].charge == 0)
741  {
742  Wnii_sum += x[k] * Wnii[i][j][k];
743  Unii_sum += 2 * x[k] * Unii[i][j][k];
744  Vnii_sum += 3 * x[k]* x[k] * Vnii[i][j][k];
745  }
746  }
747 
748  lnfn = lnfn + 0.5 * ((1/F) * E[i] * E[j] * (fabs(Z[i]) + fabs(Z[j])) / (fabs(Z[i]) * fabs(Z[j])) * (Wnii[index][i][j] - Wnii_sum) +
749  x[i] * x[j] * pow (fabs(Z[i]) + fabs(Z[j]), 2) / (fabs(Z[i]) * fabs(Z[j])) * (Unii[index][i][j] - Unii_sum) +
750  4 * x[i] * x[j] * (2 * x[index] * Vnii[index][i][j] - Vnii_sum));
751  }
752  }
753  }
754 
755  // ion-ion-ion ternary interacton
756  for (i = 0; i < aq.species_num; i ++)
757  {
758  double sum = 0;
759 
760  for (j = 0; j < aq.species_num; j ++)
761  {
762  for (l = 0; l < aq.species_num; l ++)
763  {
764  if ((Z[j] * Z[l] > 0) && (Z[j] * Z[i] < 0))
765  {
766  lnfn = lnfn + (- 2 * E[i] * x[j] * x[l] * Wiii[i][j][l]) * 0.5;
767  }
768  }
769  }
770  }
771 
772  // neutral-ion-ion-ion quaternary interacton
773  for (i = 0; i < aq.species_num; i ++)
774  {
775  double sum = 0;
776 
777  for (j = 0; j < aq.species_num; j ++)
778  {
779  for (l = 0; l < aq.species_num; l ++)
780  {
781  if ((Z[j] * Z[l] > 0) && (Z[j] * Z[i] < 0))
782  {
783  double Qniii_sum = 0;
784  double Yniii_sum = 0;
785 
786  for (k = 0; k < aq.species_num; k ++)
787  {
788  if (aq.aqueous_species[k].charge == 0)
789  {
790  Qniii_sum += 2 * x[k] * Qniii[i][j][l][k];
791  Yniii_sum += 3 * x[k] * x[k] * Yniii[i][j][l][k];
792  }
793  }
794 
795  lnfn = lnfn + 0.5 * (4 * E[i] * x[j] * x[l] * (Qniii[index][i][j][l] - Qniii_sum)
796  + 4 * E[i] * x[j] * x[l] * (2 * x[index] * Yniii[index][i][j][l] - Yniii_sum));
797  }
798  }
799  }
800  }
801 
802  // neutral-neutral binary interaction
803  for (i = 0; i < aq.species_num; i ++)
804  {
805  if (Z[i] == 0 && i != index)
806  {
807  lnfn = lnfn + x[i] * ((1 - x[index]) * Wnn[index][i] + (2 * (x[index] - x[i]) * (1 - x[index]) + x[i]) * Unn[index][i]);
808  }
809  }
810 
811  for (i = 0; i < aq.species_num; i ++)
812  {
813  for (j = 0; j < i; j ++)
814  {
815  if (Z[i] == 0 && Z[j] == 0 && i != index && j != index && i != j)
816  {
817  lnfn = lnfn - x[i] * x[j] * (Wnn[i][j] + 2 * (x[i] - x[j]) * Unn[i][j]);
818  }
819  }
820  }
821 
822  for (i = 0; i < aq.species_num; i ++)
823  {
824  for (j = 0; j < aq.species_num; j ++)
825  {
826  double sum1 = 0;
827  double sum2 = 0;
828 
829  for (k = 0; k < aq.species_num; k ++)
830  {
831  if (Z[i] * Z[j] < 0 && Z[k] == 0 && k != index)
832  {
833  sum1 += x[k] * Ynnii[index][i][j][k];
834  }
835  }
836 
837  for (k = 0; k < aq.species_num; k ++)
838  {
839  for (l = 0; l < aq.species_num; l ++)
840  {
841  if (Z[i] * Z[j] < 0 && Z[k] == 0 && Z[l] == 0 && k != l)
842  {
843  sum2 += 2.0 * 0.5 * x[k] * x[l] * Ynnii[i][j][k][l];
844  }
845  }
846  }
847 
848  if (Z[i] * Z[j] < 0)
849  {
850  lnfn = lnfn + 0.5 * (1/F) * E[i] * E[j] * (fabs(Z[i]) + fabs(Z[j]))/(fabs(Z[i]) * fabs(Z[j])) * (sum1 - sum2);
851  }
852  }
853  }
854 
855  for (k = 0; k < aq.species_num; k ++)
856  {
857  if (Z[k] == 0 && strcmp (aq.aqueous_species[k].species_symbol, WATER_SYMBOL) == 0)
858  {
859  break;
860  }
861  }
862 
863  lnfn = lnfn - (Wnn[index][k] - Unn[index][k]);
864 
865  if (index != k) {return exp(lnfn) * x[k];}
866 
867  return exp(lnfn);
868 }
869 
870 static double psc_an (AQUEOUS_PHASE aq, int index, double *n, double T, double P)
871 {
872  if (aq.aqueous_species[index].charge != 0)
873  {
874  printf ("The species is not a neutral species!!!\n");
875  return 1.0;
876  }
877 
878  if (index >= aq.species_num)
879  {
880  printf ("The species dose not existed!!!\n");
881  return 1.0;
882  }
883 
884 
885  double lnfn;
886  double x[aq.species_num];
887  double Z[aq.species_num];
888  double E[aq.species_num];
889  int i, j, k, l, m, s;
890 
891  double sumn;
892  double Ax, Ix, F;
893 
894  double aii[aq.species_num][aq.species_num];
895  double a1ii[aq.species_num][aq.species_num];
896  double Bii[aq.species_num][aq.species_num];
897  double B1ii[aq.species_num][aq.species_num];
898  double Wnn[aq.species_num][aq.species_num];
899  double Unn[aq.species_num][aq.species_num];
900  double Wnii[aq.species_num][aq.species_num][aq.species_num];
901  double Unii[aq.species_num][aq.species_num][aq.species_num];
902  double Vnii[aq.species_num][aq.species_num][aq.species_num];
903  double Wiii[aq.species_num][aq.species_num][aq.species_num];
904  double Qniii[aq.species_num][aq.species_num][aq.species_num][aq.species_num];
905  double Yniii[aq.species_num][aq.species_num][aq.species_num][aq.species_num];
906  double Ynnii[aq.species_num][aq.species_num][aq.species_num][aq.species_num];
907 
908  double Bii_a, Bii_b, Bii_c, Bii_d, Bii_e, Bii_f, \
909  B1ii_a, B1ii_b, B1ii_c, B1ii_d, B1ii_e, B1ii_f, \
910  Wnn_a, Wnn_b, Wnn_c, Wnn_d, Wnn_e, Wnn_f, \
911  Unn_a, Unn_b, Unn_c, Unn_d, Unn_e, Unn_f, \
912  Wnii_a, Wnii_b, Wnii_c, Wnii_d, Wnii_e, Wnii_f, \
913  Unii_a, Unii_b, Unii_c, Unii_d, Unii_e, Unii_f, \
914  Vnii_a, Vnii_b, Vnii_c, Vnii_d, Vnii_e, Vnii_f, \
915  Wiii_a, Wiii_b, Wiii_c, Wiii_d, Wiii_e, Wiii_f, \
916  Qniii_a, Qniii_b, Qniii_c, Qniii_d, Qniii_e, Qniii_f, \
917  Yniii_a, Yniii_b, Yniii_c, Yniii_d, Yniii_e, Yniii_f, \
918  Ynnii_a, Ynnii_b, Ynnii_c, Ynnii_d, Ynnii_e, Ynnii_f;
919 
920  Ax = ax (T, 1.0);
921 
922  for (k = 0; k < aq.species_num; k ++)
923  {
924  Z[k] = aq.aqueous_species[k].charge;
925  }
926 
927  sumn = 0;
928  for (k = 0; k < aq.species_num; k ++)
929  {
930  sumn += n[k];
931  }
932 
933  for (k = 0; k < aq.species_num; k ++)
934  {
935  x[k] = n[k] / sumn;
936  }
937 
938  F = 0;
939  for (k = 0; k < aq.species_num; k++)
940  {
941  F += 0.5 * x[k] * fabs(Z[k]);
942  }
943  F = 1.0/F;
944 
945  Ix = 0;
946  for (k = 0; k < aq.species_num; k ++)
947  {
948  Ix += 0.5 * (x[k] * Z[k] * Z[k]);
949  }
950 
951  for (k = 0; k < aq.species_num; k ++)
952  {
953  E[k] = 0;
954  }
955 
956  // INITALIZATION set all binary parameters zero
957  for (i = 0; i < aq.species_num; i ++)
958  {
959  for (j = 0; j < aq.species_num; j ++)
960  {
961  aii[i][j] = 0;
962  a1ii[i][j] = 0;
963  Bii[i][j] = 0;
964  B1ii[i][j] = 0;
965  Wnn[i][j] = 0;
966  Unn[i][j] = 0;
967  }
968  }
969 
970  // INITALIZATION set all ternary parameters zero
971  for (i = 0; i < aq.species_num; i ++)
972  {
973  for (j = 0; j < aq.species_num; j ++)
974  {
975  for (l = 0; l < aq.species_num; l ++)
976  {
977  Wnii[i][j][l] = 0;
978  Unii[i][j][l] = 0;
979  Vnii[i][j][l] = 0;
980  Wiii[i][j][l] = 0;
981  }
982  }
983  }
984 
985  // INITALIZATION set all quaternary parameters zero
986  for (i = 0; i < aq.species_num; i ++)
987  {
988  for (j = 0; j < aq.species_num; j ++)
989  {
990  for (l = 0; l < aq.species_num; l ++)
991  {
992  for (m = 0; m < aq.species_num; m ++)
993  {
994  Qniii[i][j][l][m] = 0;
995  Yniii[i][j][l][m] = 0;
996  Ynnii[i][j][l][m] = 0;
997  }
998  }
999  }
1000  }
1001 
1002  for (k = 0; k < aq.ii_bin_param_num; k ++)
1003  {
1004  for (m = 0; m < aq.species_num; m ++)
1005  {
1006  if (strcmp (aq.ii_bin_params[k].species_s1, aq.aqueous_species[m].species_symbol) == 0)
1007  {
1008  i = m;
1009  }
1010  }
1011 
1012  for (m = 0; m < aq.species_num; m ++)
1013  {
1014  if (strcmp (aq.ii_bin_params[k].species_s2, aq.aqueous_species[m].species_symbol) == 0)
1015  {
1016  j = m;
1017  }
1018  }
1019 
1020  double aii_value = aq.ii_bin_params[k].alpha;
1021  double a1ii_value = aq.ii_bin_params[k].alpha1;
1022 
1023  Bii_a = aq.ii_bin_params[k].B0ss[0];
1024  Bii_b = aq.ii_bin_params[k].B0ss[1];
1025  Bii_c = aq.ii_bin_params[k].B0ss[2];
1026  Bii_d = aq.ii_bin_params[k].B0ss[3];
1027  Bii_e = aq.ii_bin_params[k].B0ss[4];
1028  Bii_f = aq.ii_bin_params[k].B0ss[5];
1029  B1ii_a = aq.ii_bin_params[k].B1ss[0];
1030  B1ii_b = aq.ii_bin_params[k].B1ss[1];
1031  B1ii_c = aq.ii_bin_params[k].B1ss[2];
1032  B1ii_d = aq.ii_bin_params[k].B1ss[3];
1033  B1ii_e = aq.ii_bin_params[k].B1ss[4];
1034  B1ii_f = aq.ii_bin_params[k].B1ss[5];
1035 
1036  double Bii_value = Bii_a + Bii_b * T + Bii_c * T * log(T) + Bii_d * T * T + Bii_e * T * T * T + Bii_f/T;
1037  double B1ii_value = B1ii_a + B1ii_b * T + B1ii_c * T * log(T) + B1ii_d * T * T + B1ii_e * T * T * T + B1ii_f/T;
1038 
1039  aii[i][j] = aii_value; aii[j][i] = aii[i][j];
1040  a1ii[i][j] = a1ii_value; a1ii[j][i] = a1ii[i][j];
1041  Bii[i][j] = Bii_value; Bii[j][i] = Bii[i][j];
1042  B1ii[i][j] = B1ii_value; B1ii[j][i] = B1ii[i][j];
1043  }
1044 
1045  for (k = 0; k < aq.nn_bin_param_num; k ++)
1046  {
1047  for (m = 0; m < aq.species_num; m ++)
1048  {
1049  if (strcmp (aq.nn_bin_params[k].species_s1, aq.aqueous_species[m].species_symbol) == 0)
1050  {
1051  i = m;
1052  }
1053  }
1054 
1055  for (m = 0; m < aq.species_num; m ++)
1056  {
1057  if (strcmp (aq.nn_bin_params[k].species_s2, aq.aqueous_species[m].species_symbol) == 0)
1058  {
1059  j = m;
1060  }
1061  }
1062 
1063  Wnn_a = aq.nn_bin_params[k].Wss[0];
1064  Wnn_b = aq.nn_bin_params[k].Wss[1];
1065  Wnn_c = aq.nn_bin_params[k].Wss[2];
1066  Wnn_d = aq.nn_bin_params[k].Wss[3];
1067  Wnn_e = aq.nn_bin_params[k].Wss[4];
1068  Wnn_f = aq.nn_bin_params[k].Wss[5];
1069  Unn_a = aq.nn_bin_params[k].Uss[0];
1070  Unn_b = aq.nn_bin_params[k].Uss[1];
1071  Unn_c = aq.nn_bin_params[k].Uss[2];
1072  Unn_d = aq.nn_bin_params[k].Uss[3];
1073  Unn_e = aq.nn_bin_params[k].Uss[4];
1074  Unn_f = aq.nn_bin_params[k].Uss[5];
1075 
1076  double Wnn_value = Wnn_a + Wnn_b * T + Wnn_c * T * log(T) + Wnn_d * T * T + Wnn_e * T * T * T + Wnn_f/T;
1077  double Unn_value = Unn_a + Unn_b * T + Unn_c * T * log(T) + Unn_d * T * T + Unn_e * T * T * T + Unn_f/T;
1078 
1079  Wnn[i][j] = Wnn_value; Wnn[j][i] = Wnn[i][j];
1080  Unn[i][j] = Unn_value; Unn[j][i] = Unn[i][j];
1081  }
1082 
1083  for (k = 0; k < aq.nii_ter_param_num; k++)
1084  {
1085  for (m = 0; m < aq.species_num; m ++)
1086  {
1087  if (strcmp (aq.nii_ter_params[k].species_s1, aq.aqueous_species[m].species_symbol) == 0)
1088  {
1089  i = m;
1090  }
1091  }
1092 
1093  for (m = 0; m < aq.species_num; m ++)
1094  {
1095  if (strcmp (aq.nii_ter_params[k].species_s2, aq.aqueous_species[m].species_symbol) == 0)
1096  {
1097  j = m;
1098  }
1099  }
1100 
1101  for (m = 0; m < aq.species_num; m ++)
1102  {
1103  if (strcmp (aq.nii_ter_params[k].species_s3, aq.aqueous_species[m].species_symbol) == 0)
1104  {
1105  l = m;
1106  }
1107  }
1108 
1109  Wnii_a = aq.nii_ter_params[k].Wsss[0];
1110  Wnii_b = aq.nii_ter_params[k].Wsss[1];
1111  Wnii_c = aq.nii_ter_params[k].Wsss[2];
1112  Wnii_d = aq.nii_ter_params[k].Wsss[3];
1113  Wnii_e = aq.nii_ter_params[k].Wsss[4];
1114  Wnii_f = aq.nii_ter_params[k].Wsss[5];
1115  Unii_a = aq.nii_ter_params[k].Usss[0];
1116  Unii_b = aq.nii_ter_params[k].Usss[1];
1117  Unii_c = aq.nii_ter_params[k].Usss[2];
1118  Unii_d = aq.nii_ter_params[k].Usss[3];
1119  Unii_e = aq.nii_ter_params[k].Usss[4];
1120  Unii_f = aq.nii_ter_params[k].Usss[5];
1121  Vnii_a = aq.nii_ter_params[k].Vsss[0];
1122  Vnii_b = aq.nii_ter_params[k].Vsss[1];
1123  Vnii_c = aq.nii_ter_params[k].Vsss[2];
1124  Vnii_d = aq.nii_ter_params[k].Vsss[3];
1125  Vnii_e = aq.nii_ter_params[k].Vsss[4];
1126  Vnii_f = aq.nii_ter_params[k].Vsss[5];
1127 
1128  double Wnii_value = Wnii_a + Wnii_b * T + Wnii_c * T * log(T) + Wnii_d * T * T + Wnii_e * T * T * T + Wnii_f/T;
1129  double Unii_value = Unii_a + Unii_b * T + Unii_c * T * log(T) + Unii_d * T * T + Unii_e * T * T * T + Unii_f/T;
1130  double Vnii_value = Vnii_a + Vnii_b * T + Vnii_c * T * log(T) + Vnii_d * T * T + Vnii_e * T * T * T + Vnii_f/T;
1131 
1132  Wnii[i][j][l] = Wnii_value;
1133  Wnii[i][l][j] = Wnii[i][j][l];
1134  Wnii[j][i][l] = Wnii[i][j][l];
1135  Wnii[j][l][i] = Wnii[i][j][l];
1136  Wnii[l][i][j] = Wnii[i][j][l];
1137  Wnii[l][j][i] = Wnii[i][j][l];
1138 
1139  Unii[i][j][l] = Unii_value;
1140  Unii[i][l][j] = Unii[i][j][l];
1141  Unii[j][i][l] = Unii[i][j][l];
1142  Unii[j][l][i] = Unii[i][j][l];
1143  Unii[l][i][j] = Unii[i][j][l];
1144  Unii[l][j][i] = Unii[i][j][l];
1145 
1146  Vnii[i][j][l] = Vnii_value;
1147  Vnii[i][l][j] = Vnii[i][j][l];
1148  Vnii[j][i][l] = Vnii[i][j][l];
1149  Vnii[j][l][i] = Vnii[i][j][l];
1150  Vnii[l][i][j] = Vnii[i][j][l];
1151  Vnii[l][j][i] = Vnii[i][j][l];
1152  }
1153 
1154  for (k = 0; k < aq.iii_ter_param_num; k++)
1155  {
1156  for (m = 0; m < aq.species_num; m ++)
1157  {
1158  if (strcmp (aq.iii_ter_params[k].species_s1, aq.aqueous_species[m].species_symbol) == 0)
1159  {
1160  i = m;
1161  }
1162  }
1163 
1164  for (m = 0; m < aq.species_num; m ++)
1165  {
1166  if (strcmp (aq.iii_ter_params[k].species_s2, aq.aqueous_species[m].species_symbol) == 0)
1167  {
1168  j = m;
1169  }
1170  }
1171 
1172  for (m = 0; m < aq.species_num; m ++)
1173  {
1174  if (strcmp (aq.iii_ter_params[k].species_s3, aq.aqueous_species[m].species_symbol) == 0)
1175  {
1176  l = m;
1177  }
1178  }
1179 
1180  Wiii_a = aq.iii_ter_params[k].Wsss[0];
1181  Wiii_b = aq.iii_ter_params[k].Wsss[1];
1182  Wiii_c = aq.iii_ter_params[k].Wsss[2];
1183  Wiii_d = aq.iii_ter_params[k].Wsss[3];
1184  Wiii_e = aq.iii_ter_params[k].Wsss[4];
1185  Wiii_f = aq.iii_ter_params[k].Wsss[5];
1186 
1187  double Wiii_value = Wiii_a + Wiii_b * T + Wiii_c * T * log(T) + Wiii_d * T * T + Wiii_e * T * T * T + Wiii_f/T;
1188 
1189  Wiii[i][j][l] = Wiii_value;
1190  Wiii[i][l][j] = Wiii[i][j][l];
1191  Wiii[j][i][l] = Wiii[i][j][l];
1192  Wiii[j][l][i] = Wiii[i][j][l];
1193  Wiii[l][i][j] = Wiii[i][j][l];
1194  Wiii[l][j][i] = Wiii[i][j][l];
1195  }
1196 
1197  for (k = 0; k < aq.niii_qua_param_num; k++)
1198  {
1199  for (m = 0; m < aq.species_num; m ++)
1200  {
1201  if (strcmp (aq.niii_qua_params[k].species_s1, aq.aqueous_species[m].species_symbol) == 0)
1202  {
1203  i = m;
1204  }
1205  }
1206 
1207  for (m = 0; m < aq.species_num; m ++)
1208  {
1209  if (strcmp (aq.niii_qua_params[k].species_s2, aq.aqueous_species[m].species_symbol) == 0)
1210  {
1211  j = m;
1212  }
1213  }
1214 
1215  for (m = 0; m < aq.species_num; m ++)
1216  {
1217  if (strcmp (aq.niii_qua_params[k].species_s3, aq.aqueous_species[m].species_symbol) == 0)
1218  {
1219  l = m;
1220  }
1221  }
1222 
1223  for (m = 0; m < aq.species_num; m ++)
1224  {
1225  if (strcmp (aq.niii_qua_params[k].species_s4, aq.aqueous_species[m].species_symbol) == 0)
1226  {
1227  s = m;
1228  }
1229  }
1230  Qniii_a = aq.niii_qua_params[k].Qssss[0];
1231  Qniii_b = aq.niii_qua_params[k].Qssss[1];
1232  Qniii_c = aq.niii_qua_params[k].Qssss[2];
1233  Qniii_d = aq.niii_qua_params[k].Qssss[3];
1234  Qniii_e = aq.niii_qua_params[k].Qssss[4];
1235  Qniii_f = aq.niii_qua_params[k].Qssss[5];
1236  Yniii_a = aq.niii_qua_params[k].Yssss[0];
1237  Yniii_b = aq.niii_qua_params[k].Yssss[1];
1238  Yniii_c = aq.niii_qua_params[k].Yssss[2];
1239  Yniii_d = aq.niii_qua_params[k].Yssss[3];
1240  Yniii_e = aq.niii_qua_params[k].Yssss[4];
1241  Yniii_f = aq.niii_qua_params[k].Yssss[5];
1242 
1243  double Qniii_value = Qniii_a + Qniii_b * T + Qniii_c * T * log(T) + Qniii_d * T * T + Qniii_e * T * T * T + Qniii_f/T;
1244  double Yniii_value = Yniii_a + Yniii_b * T + Yniii_c * T * log(T) + Yniii_d * T * T + Yniii_e * T * T * T + Yniii_f/T;
1245 
1246  Qniii[i][j][l][s] = Qniii_value;
1247  Qniii[i][j][s][l] = Qniii[i][j][l][s];
1248  Qniii[i][l][j][s] = Qniii[i][j][l][s];
1249  Qniii[i][l][s][j] = Qniii[i][j][l][s];
1250  Qniii[i][s][j][l] = Qniii[i][j][l][s];
1251  Qniii[i][s][l][j] = Qniii[i][j][l][s];
1252 
1253  Qniii[j][i][l][s] = Qniii[i][j][l][s];
1254  Qniii[j][i][s][l] = Qniii[i][j][l][s];
1255  Qniii[j][l][i][s] = Qniii[i][j][l][s];
1256  Qniii[j][l][s][i] = Qniii[i][j][l][s];
1257  Qniii[j][s][i][l] = Qniii[i][j][l][s];
1258  Qniii[j][s][l][i] = Qniii[i][j][l][s];
1259 
1260  Qniii[l][i][j][s] = Qniii[i][j][l][s];
1261  Qniii[l][i][s][j] = Qniii[i][j][l][s];
1262  Qniii[l][j][i][s] = Qniii[i][j][l][s];
1263  Qniii[l][j][s][i] = Qniii[i][j][l][s];
1264  Qniii[l][s][i][j] = Qniii[i][j][l][s];
1265  Qniii[l][s][j][i] = Qniii[i][j][l][s];
1266 
1267  Qniii[s][i][j][l] = Qniii[i][j][l][s];
1268  Qniii[s][i][l][j] = Qniii[i][j][l][s];
1269  Qniii[s][j][l][i] = Qniii[i][j][l][s];
1270  Qniii[s][j][i][l] = Qniii[i][j][l][s];
1271  Qniii[s][l][i][j] = Qniii[i][j][l][s];
1272  Qniii[s][l][j][i] = Qniii[i][j][l][s];
1273 
1274  Yniii[i][j][l][s] = Yniii_value;
1275  Yniii[i][j][s][l] = Yniii[i][j][l][s];
1276  Yniii[i][l][j][s] = Yniii[i][j][l][s];
1277  Yniii[i][l][s][j] = Yniii[i][j][l][s];
1278  Yniii[i][s][j][l] = Yniii[i][j][l][s];
1279  Yniii[i][s][l][j] = Yniii[i][j][l][s];
1280 
1281  Yniii[j][i][l][s] = Yniii[i][j][l][s];
1282  Yniii[j][i][s][l] = Yniii[i][j][l][s];
1283  Yniii[j][l][i][s] = Yniii[i][j][l][s];
1284  Yniii[j][l][s][i] = Yniii[i][j][l][s];
1285  Yniii[j][s][i][l] = Yniii[i][j][l][s];
1286  Yniii[j][s][l][i] = Yniii[i][j][l][s];
1287 
1288  Yniii[l][i][j][s] = Yniii[i][j][l][s];
1289  Yniii[l][i][s][j] = Yniii[i][j][l][s];
1290  Yniii[l][j][i][s] = Yniii[i][j][l][s];
1291  Yniii[l][j][s][i] = Yniii[i][j][l][s];
1292  Yniii[l][s][i][j] = Yniii[i][j][l][s];
1293  Yniii[l][s][j][i] = Yniii[i][j][l][s];
1294 
1295  Yniii[s][i][j][l] = Yniii[i][j][l][s];
1296  Yniii[s][i][l][j] = Yniii[i][j][l][s];
1297  Yniii[s][j][l][i] = Yniii[i][j][l][s];
1298  Yniii[s][j][i][l] = Yniii[i][j][l][s];
1299  Yniii[s][l][i][j] = Yniii[i][j][l][s];
1300  Yniii[s][l][j][i] = Yniii[i][j][l][s];
1301  }
1302 
1303  for (k = 0; k < aq.nnii_qua_param_num; k++)
1304  {
1305  for (m = 0; m < aq.species_num; m ++)
1306  {
1307  if (strcmp (aq.nnii_qua_params[k].species_s1, aq.aqueous_species[m].species_symbol) == 0)
1308  {
1309  i = m;
1310  }
1311  }
1312 
1313  for (m = 0; m < aq.species_num; m ++)
1314  {
1315  if (strcmp (aq.nnii_qua_params[k].species_s2, aq.aqueous_species[m].species_symbol) == 0)
1316  {
1317  j = m;
1318  }
1319  }
1320 
1321  for (m = 0; m < aq.species_num; m ++)
1322  {
1323  if (strcmp (aq.nnii_qua_params[k].species_s3, aq.aqueous_species[m].species_symbol) == 0)
1324  {
1325  l = m;
1326  }
1327  }
1328 
1329  for (m = 0; m < aq.species_num; m ++)
1330  {
1331  if (strcmp (aq.nnii_qua_params[k].species_s4, aq.aqueous_species[m].species_symbol) == 0)
1332  {
1333  s = m;
1334  }
1335  }
1336 
1337  Ynnii_a = aq.nnii_qua_params[k].Yssss[0];
1338  Ynnii_b = aq.nnii_qua_params[k].Yssss[1];
1339  Ynnii_c = aq.nnii_qua_params[k].Yssss[2];
1340  Ynnii_d = aq.nnii_qua_params[k].Yssss[3];
1341  Ynnii_e = aq.nnii_qua_params[k].Yssss[4];
1342  Ynnii_f = aq.nnii_qua_params[k].Yssss[5];
1343 
1344  double Ynnii_value = Ynnii_a + Ynnii_b * T + Ynnii_c * T * log(T) + Ynnii_d * T * T + Ynnii_e * T * T * T + Ynnii_f/T;
1345 
1346  Ynnii[i][j][l][s] = Ynnii_value;
1347  Ynnii[i][j][s][l] = Ynnii[i][j][l][s];
1348  Ynnii[i][l][j][s] = Ynnii[i][j][l][s];
1349  Ynnii[i][l][s][j] = Ynnii[i][j][l][s];
1350  Ynnii[i][s][j][l] = Ynnii[i][j][l][s];
1351  Ynnii[i][s][l][j] = Ynnii[i][j][l][s];
1352 
1353  Ynnii[j][i][l][s] = Ynnii[i][j][l][s];
1354  Ynnii[j][i][s][l] = Ynnii[i][j][l][s];
1355  Ynnii[j][l][i][s] = Ynnii[i][j][l][s];
1356  Ynnii[j][l][s][i] = Ynnii[i][j][l][s];
1357  Ynnii[j][s][i][l] = Ynnii[i][j][l][s];
1358  Ynnii[j][s][l][i] = Ynnii[i][j][l][s];
1359 
1360  Ynnii[l][i][j][s] = Ynnii[i][j][l][s];
1361  Ynnii[l][i][s][j] = Ynnii[i][j][l][s];
1362  Ynnii[l][j][i][s] = Ynnii[i][j][l][s];
1363  Ynnii[l][j][s][i] = Ynnii[i][j][l][s];
1364  Ynnii[l][s][i][j] = Ynnii[i][j][l][s];
1365  Ynnii[l][s][j][i] = Ynnii[i][j][l][s];
1366 
1367  Ynnii[s][i][j][l] = Ynnii[i][j][l][s];
1368  Ynnii[s][i][l][j] = Ynnii[i][j][l][s];
1369  Ynnii[s][j][l][i] = Ynnii[i][j][l][s];
1370  Ynnii[s][j][i][l] = Ynnii[i][j][l][s];
1371  Ynnii[s][l][i][j] = Ynnii[i][j][l][s];
1372  Ynnii[s][l][j][i] = Ynnii[i][j][l][s];
1373  }
1374 
1375  for (i = 0; i < aq.species_num; i ++)
1376  {
1377  E[i] = x[i] * fabs (Z[i]);
1378  double sum = 0;
1379 
1380  for (k = 0; k < aq.species_num; k++)
1381  {
1382  if (Z[i] * Z[k] > 0) {sum = sum + x[k] * fabs (Z[k]);}
1383  }
1384 
1385  if (sum != 0) {E[i] = E[i] / sum;}
1386  }
1387 
1388  // DH term
1389  lnfn = 2.0 * Ax * Ix * sqrt (Ix) / (1.0 + ROU * sqrt (Ix));
1390 
1391 
1392  // cation-anion binary interaction
1393  for (i = 0; i < aq.species_num; i ++)
1394  {
1395  for (j = 0; j < aq.species_num; j ++)
1396  {
1397  if (Z[i]*Z[j] != 0)
1398  {
1399  lnfn = lnfn - 0.5 * (x[i] * x[j] * Bii[i][j] * exp (-aii[i][j] * sqrt (Ix)) + x[i] * x[j] * B1ii[i][j] * exp (-a1ii[i][j] * sqrt (Ix)));
1400  }
1401  }
1402  }
1403 
1404  // cation-cation or anion-anion ternary HOE term
1405  for (i = 0; i < aq.species_num; i ++)
1406  {
1407  for (j = 0; j < aq.species_num; j ++)
1408  {
1409  if (Z[i] * Z[j] > 0 && i != j && Ix != 0)
1410  {
1411  double xij = 6.0 * Z[i] * Z[j] * Ax * sqrt (Ix);
1412  double xii = 6.0 * Z[i] * Z[i] * Ax * sqrt (Ix);
1413  double xjj = 6.0 * Z[j] * Z[j] * Ax * sqrt (Ix);
1414 
1415  double vij = Z[i] * Z[j] / (4.0 * Ix) * (jj(xij) - 0.5 * jj(xii) - 0.5 * jj(xjj));
1416  double vpij = - vij/Ix + (Z[i]*Z[j]/(8.0 * Ix * Ix)) * (xij * jp(xij) - 0.5 * xii * jp(xii) - 0.5 * xjj * jp(xjj));
1417 
1418  lnfn = lnfn + 0.5 * (-2.0 * x[i] * x[j] * (vij + Ix * vpij));
1419  }
1420  }
1421  }
1422 
1423  // neutral-caion-anion ternary interaction
1424  for (i = 0; i < aq.species_num; i ++)
1425  {
1426  for (j = 0; j < aq.species_num; j ++)
1427  {
1428  double sum;
1429 
1430  if (Z[i]*Z[j] < 0)
1431  {
1432  double Wnii_sum = 0;
1433  double Unii_sum = 0;
1434  double Vnii_sum = 0;
1435 
1436  for (k = 0; k < aq.species_num; k ++)
1437  {
1438  if (aq.aqueous_species[k].charge == 0)
1439  {
1440  Wnii_sum += x[k] * Wnii[i][j][k];
1441  Unii_sum += 2 * x[k] * Unii[i][j][k];
1442  Vnii_sum += 3 * x[k]* x[k] * Vnii[i][j][k];
1443  }
1444  }
1445 
1446  lnfn = lnfn + 0.5 * ((1/F) * E[i] * E[j] * (fabs(Z[i]) + fabs(Z[j])) / (fabs(Z[i]) * fabs(Z[j])) * (Wnii[index][i][j] - Wnii_sum) +
1447  x[i] * x[j] * pow (fabs(Z[i]) + fabs(Z[j]), 2) / (fabs(Z[i]) * fabs(Z[j])) * (Unii[index][i][j] - Unii_sum) +
1448  4 * x[i] * x[j] * (2 * x[index] * Vnii[index][i][j] - Vnii_sum));
1449  }
1450  }
1451  }
1452 
1453  // ion-ion-ion ternary interacton
1454  for (i = 0; i < aq.species_num; i ++)
1455  {
1456  double sum = 0;
1457 
1458  for (j = 0; j < aq.species_num; j ++)
1459  {
1460  for (l = 0; l < aq.species_num; l ++)
1461  {
1462  if ((Z[j] * Z[l] > 0) && (Z[j] * Z[i] < 0))
1463  {
1464  lnfn = lnfn + (- 2 * E[i] * x[j] * x[l] * Wiii[i][j][l]) * 0.5;
1465  }
1466  }
1467  }
1468  }
1469 
1470  // neutral-ion-ion-ion quaternary interacton
1471  for (i = 0; i < aq.species_num; i ++)
1472  {
1473  double sum = 0;
1474 
1475  for (j = 0; j < aq.species_num; j ++)
1476  {
1477  for (l = 0; l < aq.species_num; l ++)
1478  {
1479  if ((Z[j] * Z[l] > 0) && (Z[j] * Z[i] < 0))
1480  {
1481  double Qniii_sum = 0;
1482  double Yniii_sum = 0;
1483 
1484  for (k = 0; k < aq.species_num; k ++)
1485  {
1486  if (aq.aqueous_species[k].charge == 0)
1487  {
1488  Qniii_sum += 2 * x[k] * Qniii[i][j][l][k];
1489  Yniii_sum += 3 * x[k] * x[k] * Yniii[i][j][l][k];
1490  }
1491  }
1492 
1493  lnfn = lnfn + 0.5 * (4 * E[i] * x[j] * x[l] * (Qniii[index][i][j][l] - Qniii_sum)
1494  + 4 * E[i] * x[j] * x[l] * (2 * x[index] * Yniii[index][i][j][l] - Yniii_sum));
1495  }
1496  }
1497  }
1498  }
1499 
1500  // neutral-neutral binary interaction
1501  for (i = 0; i < aq.species_num; i ++)
1502  {
1503  if (Z[i] == 0 && i != index)
1504  {
1505  lnfn = lnfn + x[i] * ((1 - x[index]) * Wnn[index][i] + (2 * (x[index] - x[i]) * (1 - x[index]) + x[i]) * Unn[index][i]);
1506  }
1507  }
1508 
1509  for (i = 0; i < aq.species_num; i ++)
1510  {
1511  for (j = 0; j < i; j ++)
1512  {
1513  if (Z[i] == 0 && Z[j] == 0 && i != index && j != index && i != j)
1514  {
1515  lnfn = lnfn - x[i] * x[j] * (Wnn[i][j] + 2 * (x[i] - x[j]) * Unn[i][j]);
1516  }
1517  }
1518  }
1519 
1520  for (i = 0; i < aq.species_num; i ++)
1521  {
1522  for (j = 0; j < aq.species_num; j ++)
1523  {
1524  double sum1 = 0;
1525  double sum2 = 0;
1526 
1527  for (k = 0; k < aq.species_num; k ++)
1528  {
1529  if (Z[i] * Z[j] < 0 && Z[k] == 0 && k != index)
1530  {
1531  sum1 += x[k] * Ynnii[index][i][j][k];
1532  }
1533  }
1534 
1535  for (k = 0; k < aq.species_num; k ++)
1536  {
1537  for (l = 0; l < aq.species_num; l ++)
1538  {
1539  if (Z[i] * Z[j] < 0 && Z[k] == 0 && Z[l] == 0 && k != l)
1540  {
1541  sum2 += 2.0 * 0.5 * x[k] * x[l] * Ynnii[i][j][k][l];
1542  }
1543  }
1544  }
1545 
1546  if (Z[i] * Z[j] < 0)
1547  {
1548  lnfn = lnfn + 0.5 * (1/F) * E[i] * E[j] * (fabs(Z[i]) + fabs(Z[j]))/(fabs(Z[i]) * fabs(Z[j])) * (sum1 - sum2);
1549  }
1550  }
1551  }
1552 
1553  for (k = 0; k < aq.species_num; k ++)
1554  {
1555  if (Z[k] == 0 && strcmp (aq.aqueous_species[k].species_symbol, WATER_SYMBOL) == 0)
1556  {
1557  break;
1558  }
1559  }
1560 
1561  lnfn = lnfn - (Wnn[index][k] - Unn[index][k]);
1562 
1563  if (index != k) {return exp(lnfn) * x[k] * (n[index] / n[k] * 1000/CONSTANT_MW);}
1564 
1565  return exp(lnfn) * x[index];
1566 }
1567 
1568 static double psc_ri (AQUEOUS_PHASE aq, int index, double *n, double T, double P)
1569 {
1570  if (aq.aqueous_species[index].charge == 0)
1571  {
1572  printf ("The species is not a ion species!!!\n");
1573  return 1.0;
1574  }
1575 
1576  if (index >= aq.species_num)
1577  {
1578  printf ("The species dose not existed!!!\n");
1579  return 1.0;
1580  }
1581 
1582  double lnfi;
1583  double x[aq.species_num];
1584  double Z[aq.species_num];
1585  double E[aq.species_num];
1586  double iEj[aq.species_num][aq.species_num];
1587  int i, j, k, l, m, s;
1588 
1589  double sumn;
1590  double Ax, Ix, F;
1591 
1592  double aii[aq.species_num][aq.species_num];
1593  double a1ii[aq.species_num][aq.species_num];
1594  double Bii[aq.species_num][aq.species_num];
1595  double B1ii[aq.species_num][aq.species_num];
1596  double Wnn[aq.species_num][aq.species_num];
1597  double Unn[aq.species_num][aq.species_num];
1598  double Wnii[aq.species_num][aq.species_num][aq.species_num];
1599  double Unii[aq.species_num][aq.species_num][aq.species_num];
1600  double Vnii[aq.species_num][aq.species_num][aq.species_num];
1601  double Wiii[aq.species_num][aq.species_num][aq.species_num];
1602  double Qniii[aq.species_num][aq.species_num][aq.species_num][aq.species_num];
1603  double Yniii[aq.species_num][aq.species_num][aq.species_num][aq.species_num];
1604  double Ynnii[aq.species_num][aq.species_num][aq.species_num][aq.species_num];
1605 
1606  double Bii_a, Bii_b, Bii_c, Bii_d, Bii_e, Bii_f, \
1607  B1ii_a, B1ii_b, B1ii_c, B1ii_d, B1ii_e, B1ii_f, \
1608  Wnn_a, Wnn_b, Wnn_c, Wnn_d, Wnn_e, Wnn_f, \
1609  Unn_a, Unn_b, Unn_c, Unn_d, Unn_e, Unn_f, \
1610  Wnii_a, Wnii_b, Wnii_c, Wnii_d, Wnii_e, Wnii_f, \
1611  Unii_a, Unii_b, Unii_c, Unii_d, Unii_e, Unii_f, \
1612  Vnii_a, Vnii_b, Vnii_c, Vnii_d, Vnii_e, Vnii_f, \
1613  Wiii_a, Wiii_b, Wiii_c, Wiii_d, Wiii_e, Wiii_f, \
1614  Qniii_a, Qniii_b, Qniii_c, Qniii_d, Qniii_e, Qniii_f, \
1615  Yniii_a, Yniii_b, Yniii_c, Yniii_d, Yniii_e, Yniii_f, \
1616  Ynnii_a, Ynnii_b, Ynnii_c, Ynnii_d, Ynnii_e, Ynnii_f;
1617 
1618  Ax = ax (T, 1.0);
1619 
1620  for (k = 0; k < aq.species_num; k ++)
1621  {
1622  Z[k] = aq.aqueous_species[k].charge;
1623  }
1624 
1625  sumn = 0;
1626  for (k = 0; k < aq.species_num; k ++)
1627  {
1628  sumn += n[k];
1629  }
1630 
1631  for (k = 0; k < aq.species_num; k ++)
1632  {
1633  x[k] = n[k] / sumn;
1634  }
1635 
1636  F = 0;
1637  for (k = 0; k < aq.species_num; k++)
1638  {
1639  F += 0.5 * x[k] * fabs(Z[k]);
1640  }
1641  F = 1.0/F;
1642 
1643  Ix = 0;
1644  for (k = 0; k < aq.species_num; k ++)
1645  {
1646  Ix += 0.5 * (x[k] * Z[k] * Z[k]);
1647  }
1648 
1649  for (k = 0; k < aq.species_num; k ++)
1650  {
1651  E[k] = 0;
1652  }
1653 
1654  // INITALIZATION set all binary parameters zero
1655  for (i = 0; i < aq.species_num; i ++)
1656  {
1657  for (j = 0; j < aq.species_num; j ++)
1658  {
1659  aii[i][j] = 0;
1660  a1ii[i][j] = 0;
1661  Bii[i][j] = 0;
1662  B1ii[i][j] = 0;
1663  Wnn[i][j] = 0;
1664  Unn[i][j] = 0;
1665  }
1666  }
1667 
1668  // INITALIZATION set all ternary parameters zero
1669  for (i = 0; i < aq.species_num; i ++)
1670  {
1671  for (j = 0; j < aq.species_num; j ++)
1672  {
1673  for (l = 0; l < aq.species_num; l ++)
1674  {
1675  Wnii[i][j][l] = 0;
1676  Unii[i][j][l] = 0;
1677  Vnii[i][j][l] = 0;
1678  Wiii[i][j][l] = 0;
1679  }
1680  }
1681  }
1682 
1683  // INITALIZATION set all quaternary parameters zero
1684  for (i = 0; i < aq.species_num; i ++)
1685  {
1686  for (j = 0; j < aq.species_num; j ++)
1687  {
1688  for (l = 0; l < aq.species_num; l ++)
1689  {
1690  for (m = 0; m < aq.species_num; m ++)
1691  {
1692  Qniii[i][j][l][m] = 0;
1693  Yniii[i][j][l][m] = 0;
1694  Ynnii[i][j][l][m] = 0;
1695  }
1696  }
1697  }
1698  }
1699 
1700  for (k = 0; k < aq.ii_bin_param_num; k ++)
1701  {
1702  for (m = 0; m < aq.species_num; m ++)
1703  {
1704  if (strcmp (aq.ii_bin_params[k].species_s1, aq.aqueous_species[m].species_symbol) == 0)
1705  {
1706  i = m;
1707  }
1708  }
1709 
1710  for (m = 0; m < aq.species_num; m ++)
1711  {
1712  if (strcmp (aq.ii_bin_params[k].species_s2, aq.aqueous_species[m].species_symbol) == 0)
1713  {
1714  j = m;
1715  }
1716  }
1717 
1718  double aii_value = aq.ii_bin_params[k].alpha;
1719  double a1ii_value = aq.ii_bin_params[k].alpha1;
1720 
1721  Bii_a = aq.ii_bin_params[k].B0ss[0];
1722  Bii_b = aq.ii_bin_params[k].B0ss[1];
1723  Bii_c = aq.ii_bin_params[k].B0ss[2];
1724  Bii_d = aq.ii_bin_params[k].B0ss[3];
1725  Bii_e = aq.ii_bin_params[k].B0ss[4];
1726  Bii_f = aq.ii_bin_params[k].B0ss[5];
1727  B1ii_a = aq.ii_bin_params[k].B1ss[0];
1728  B1ii_b = aq.ii_bin_params[k].B1ss[1];
1729  B1ii_c = aq.ii_bin_params[k].B1ss[2];
1730  B1ii_d = aq.ii_bin_params[k].B1ss[3];
1731  B1ii_e = aq.ii_bin_params[k].B1ss[4];
1732  B1ii_f = aq.ii_bin_params[k].B1ss[5];
1733 
1734  double Bii_value = Bii_a + Bii_b * T + Bii_c * T * log(T) + Bii_d * T * T + Bii_e * T * T * T + Bii_f/T;
1735  double B1ii_value = B1ii_a + B1ii_b * T + B1ii_c * T * log(T) + B1ii_d * T * T + B1ii_e * T * T * T + B1ii_f/T;
1736 
1737  aii[i][j] = aii_value; aii[j][i] = aii[i][j];
1738  a1ii[i][j] = a1ii_value; a1ii[j][i] = a1ii[i][j];
1739  Bii[i][j] = Bii_value; Bii[j][i] = Bii[i][j];
1740  B1ii[i][j] = B1ii_value; B1ii[j][i] = B1ii[i][j];
1741  }
1742 
1743  for (k = 0; k < aq.nn_bin_param_num; k ++)
1744  {
1745  for (m = 0; m < aq.species_num; m ++)
1746  {
1747  if (strcmp (aq.nn_bin_params[k].species_s1, aq.aqueous_species[m].species_symbol) == 0)
1748  {
1749  i = m;
1750  }
1751  }
1752 
1753  for (m = 0; m < aq.species_num; m ++)
1754  {
1755  if (strcmp (aq.nn_bin_params[k].species_s2, aq.aqueous_species[m].species_symbol) == 0)
1756  {
1757  j = m;
1758  }
1759  }
1760 
1761  Wnn_a = aq.nn_bin_params[k].Wss[0];
1762  Wnn_b = aq.nn_bin_params[k].Wss[1];
1763  Wnn_c = aq.nn_bin_params[k].Wss[2];
1764  Wnn_d = aq.nn_bin_params[k].Wss[3];
1765  Wnn_e = aq.nn_bin_params[k].Wss[4];
1766  Wnn_f = aq.nn_bin_params[k].Wss[5];
1767  Unn_a = aq.nn_bin_params[k].Uss[0];
1768  Unn_b = aq.nn_bin_params[k].Uss[1];
1769  Unn_c = aq.nn_bin_params[k].Uss[2];
1770  Unn_d = aq.nn_bin_params[k].Uss[3];
1771  Unn_e = aq.nn_bin_params[k].Uss[4];
1772  Unn_f = aq.nn_bin_params[k].Uss[5];
1773 
1774  double Wnn_value = Wnn_a + Wnn_b * T + Wnn_c * T * log(T) + Wnn_d * T * T + Wnn_e * T * T * T + Wnn_f/T;
1775  double Unn_value = Unn_a + Unn_b * T + Unn_c * T * log(T) + Unn_d * T * T + Unn_e * T * T * T + Unn_f/T;
1776 
1777  Wnn[i][j] = Wnn_value; Wnn[j][i] = Wnn[i][j];
1778  Unn[i][j] = Unn_value; Unn[j][i] = Unn[i][j];
1779  }
1780 
1781  for (k = 0; k < aq.nii_ter_param_num; k++)
1782  {
1783  for (m = 0; m < aq.species_num; m ++)
1784  {
1785  if (strcmp (aq.nii_ter_params[k].species_s1, aq.aqueous_species[m].species_symbol) == 0)
1786  {
1787  i = m;
1788  }
1789  }
1790 
1791  for (m = 0; m < aq.species_num; m ++)
1792  {
1793  if (strcmp (aq.nii_ter_params[k].species_s2, aq.aqueous_species[m].species_symbol) == 0)
1794  {
1795  j = m;
1796  }
1797  }
1798 
1799  for (m = 0; m < aq.species_num; m ++)
1800  {
1801  if (strcmp (aq.nii_ter_params[k].species_s3, aq.aqueous_species[m].species_symbol) == 0)
1802  {
1803  l = m;
1804  }
1805  }
1806 
1807  Wnii_a = aq.nii_ter_params[k].Wsss[0];
1808  Wnii_b = aq.nii_ter_params[k].Wsss[1];
1809  Wnii_c = aq.nii_ter_params[k].Wsss[2];
1810  Wnii_d = aq.nii_ter_params[k].Wsss[3];
1811  Wnii_e = aq.nii_ter_params[k].Wsss[4];
1812  Wnii_f = aq.nii_ter_params[k].Wsss[5];
1813  Unii_a = aq.nii_ter_params[k].Usss[0];
1814  Unii_b = aq.nii_ter_params[k].Usss[1];
1815  Unii_c = aq.nii_ter_params[k].Usss[2];
1816  Unii_d = aq.nii_ter_params[k].Usss[3];
1817  Unii_e = aq.nii_ter_params[k].Usss[4];
1818  Unii_f = aq.nii_ter_params[k].Usss[5];
1819  Vnii_a = aq.nii_ter_params[k].Vsss[0];
1820  Vnii_b = aq.nii_ter_params[k].Vsss[1];
1821  Vnii_c = aq.nii_ter_params[k].Vsss[2];
1822  Vnii_d = aq.nii_ter_params[k].Vsss[3];
1823  Vnii_e = aq.nii_ter_params[k].Vsss[4];
1824  Vnii_f = aq.nii_ter_params[k].Vsss[5];
1825 
1826  double Wnii_value = Wnii_a + Wnii_b * T + Wnii_c * T * log(T) + Wnii_d * T * T + Wnii_e * T * T * T + Wnii_f/T;
1827  double Unii_value = Unii_a + Unii_b * T + Unii_c * T * log(T) + Unii_d * T * T + Unii_e * T * T * T + Unii_f/T;
1828  double Vnii_value = Vnii_a + Vnii_b * T + Vnii_c * T * log(T) + Vnii_d * T * T + Vnii_e * T * T * T + Vnii_f/T;
1829 
1830  Wnii[i][j][l] = Wnii_value;
1831  Wnii[i][l][j] = Wnii[i][j][l];
1832  Wnii[j][i][l] = Wnii[i][j][l];
1833  Wnii[j][l][i] = Wnii[i][j][l];
1834  Wnii[l][i][j] = Wnii[i][j][l];
1835  Wnii[l][j][i] = Wnii[i][j][l];
1836 
1837  Unii[i][j][l] = Unii_value;
1838  Unii[i][l][j] = Unii[i][j][l];
1839  Unii[j][i][l] = Unii[i][j][l];
1840  Unii[j][l][i] = Unii[i][j][l];
1841  Unii[l][i][j] = Unii[i][j][l];
1842  Unii[l][j][i] = Unii[i][j][l];
1843 
1844  Vnii[i][j][l] = Vnii_value;
1845  Vnii[i][l][j] = Vnii[i][j][l];
1846  Vnii[j][i][l] = Vnii[i][j][l];
1847  Vnii[j][l][i] = Vnii[i][j][l];
1848  Vnii[l][i][j] = Vnii[i][j][l];
1849  Vnii[l][j][i] = Vnii[i][j][l];
1850  }
1851 
1852  for (k = 0; k < aq.iii_ter_param_num; k++)
1853  {
1854  for (m = 0; m < aq.species_num; m ++)
1855  {
1856  if (strcmp (aq.iii_ter_params[k].species_s1, aq.aqueous_species[m].species_symbol) == 0)
1857  {
1858  i = m;
1859  }
1860  }
1861 
1862  for (m = 0; m < aq.species_num; m ++)
1863  {
1864  if (strcmp (aq.iii_ter_params[k].species_s2, aq.aqueous_species[m].species_symbol) == 0)
1865  {
1866  j = m;
1867  }
1868  }
1869 
1870  for (m = 0; m < aq.species_num; m ++)
1871  {
1872  if (strcmp (aq.iii_ter_params[k].species_s3, aq.aqueous_species[m].species_symbol) == 0)
1873  {
1874  l = m;
1875  }
1876  }
1877 
1878  Wiii_a = aq.iii_ter_params[k].Wsss[0];
1879  Wiii_b = aq.iii_ter_params[k].Wsss[1];
1880  Wiii_c = aq.iii_ter_params[k].Wsss[2];
1881  Wiii_d = aq.iii_ter_params[k].Wsss[3];
1882  Wiii_e = aq.iii_ter_params[k].Wsss[4];
1883  Wiii_f = aq.iii_ter_params[k].Wsss[5];
1884 
1885  double Wiii_value = Wiii_a + Wiii_b * T + Wiii_c * T * log(T) + Wiii_d * T * T + Wiii_e * T * T * T + Wiii_f/T;
1886 
1887  Wiii[i][j][l] = Wiii_value;
1888  Wiii[i][l][j] = Wiii[i][j][l];
1889  Wiii[j][i][l] = Wiii[i][j][l];
1890  Wiii[j][l][i] = Wiii[i][j][l];
1891  Wiii[l][i][j] = Wiii[i][j][l];
1892  Wiii[l][j][i] = Wiii[i][j][l];
1893  }
1894 
1895  for (k = 0; k < aq.niii_qua_param_num; k++)
1896  {
1897  for (m = 0; m < aq.species_num; m ++)
1898  {
1899  if (strcmp (aq.niii_qua_params[k].species_s1, aq.aqueous_species[m].species_symbol) == 0)
1900  {
1901  i = m;
1902  }
1903  }
1904 
1905  for (m = 0; m < aq.species_num; m ++)
1906  {
1907  if (strcmp (aq.niii_qua_params[k].species_s2, aq.aqueous_species[m].species_symbol) == 0)
1908  {
1909  j = m;
1910  }
1911  }
1912 
1913  for (m = 0; m < aq.species_num; m ++)
1914  {
1915  if (strcmp (aq.niii_qua_params[k].species_s3, aq.aqueous_species[m].species_symbol) == 0)
1916  {
1917  l = m;
1918  }
1919  }
1920 
1921  for (m = 0; m < aq.species_num; m ++)
1922  {
1923  if (strcmp (aq.niii_qua_params[k].species_s4, aq.aqueous_species[m].species_symbol) == 0)
1924  {
1925  s = m;
1926  }
1927  }
1928  Qniii_a = aq.niii_qua_params[k].Qssss[0];
1929  Qniii_b = aq.niii_qua_params[k].Qssss[1];
1930  Qniii_c = aq.niii_qua_params[k].Qssss[2];
1931  Qniii_d = aq.niii_qua_params[k].Qssss[3];
1932  Qniii_e = aq.niii_qua_params[k].Qssss[4];
1933  Qniii_f = aq.niii_qua_params[k].Qssss[5];
1934  Yniii_a = aq.niii_qua_params[k].Yssss[0];
1935  Yniii_b = aq.niii_qua_params[k].Yssss[1];
1936  Yniii_c = aq.niii_qua_params[k].Yssss[2];
1937  Yniii_d = aq.niii_qua_params[k].Yssss[3];
1938  Yniii_e = aq.niii_qua_params[k].Yssss[4];
1939  Yniii_f = aq.niii_qua_params[k].Yssss[5];
1940 
1941  double Qniii_value = Qniii_a + Qniii_b * T + Qniii_c * T * log(T) + Qniii_d * T * T + Qniii_e * T * T * T + Qniii_f/T;
1942  double Yniii_value = Yniii_a + Yniii_b * T + Yniii_c * T * log(T) + Yniii_d * T * T + Yniii_e * T * T * T + Yniii_f/T;
1943 
1944  Qniii[i][j][l][s] = Qniii_value;
1945  Qniii[i][j][s][l] = Qniii[i][j][l][s];
1946  Qniii[i][l][j][s] = Qniii[i][j][l][s];
1947  Qniii[i][l][s][j] = Qniii[i][j][l][s];
1948  Qniii[i][s][j][l] = Qniii[i][j][l][s];
1949  Qniii[i][s][l][j] = Qniii[i][j][l][s];
1950 
1951  Qniii[j][i][l][s] = Qniii[i][j][l][s];
1952  Qniii[j][i][s][l] = Qniii[i][j][l][s];
1953  Qniii[j][l][i][s] = Qniii[i][j][l][s];
1954  Qniii[j][l][s][i] = Qniii[i][j][l][s];
1955  Qniii[j][s][i][l] = Qniii[i][j][l][s];
1956  Qniii[j][s][l][i] = Qniii[i][j][l][s];
1957 
1958  Qniii[l][i][j][s] = Qniii[i][j][l][s];
1959  Qniii[l][i][s][j] = Qniii[i][j][l][s];
1960  Qniii[l][j][i][s] = Qniii[i][j][l][s];
1961  Qniii[l][j][s][i] = Qniii[i][j][l][s];
1962  Qniii[l][s][i][j] = Qniii[i][j][l][s];
1963  Qniii[l][s][j][i] = Qniii[i][j][l][s];
1964 
1965  Qniii[s][i][j][l] = Qniii[i][j][l][s];
1966  Qniii[s][i][l][j] = Qniii[i][j][l][s];
1967  Qniii[s][j][l][i] = Qniii[i][j][l][s];
1968  Qniii[s][j][i][l] = Qniii[i][j][l][s];
1969  Qniii[s][l][i][j] = Qniii[i][j][l][s];
1970  Qniii[s][l][j][i] = Qniii[i][j][l][s];
1971 
1972  Yniii[i][j][l][s] = Yniii_value;
1973  Yniii[i][j][s][l] = Yniii[i][j][l][s];
1974  Yniii[i][l][j][s] = Yniii[i][j][l][s];
1975  Yniii[i][l][s][j] = Yniii[i][j][l][s];
1976  Yniii[i][s][j][l] = Yniii[i][j][l][s];
1977  Yniii[i][s][l][j] = Yniii[i][j][l][s];
1978 
1979  Yniii[j][i][l][s] = Yniii[i][j][l][s];
1980  Yniii[j][i][s][l] = Yniii[i][j][l][s];
1981  Yniii[j][l][i][s] = Yniii[i][j][l][s];
1982  Yniii[j][l][s][i] = Yniii[i][j][l][s];
1983  Yniii[j][s][i][l] = Yniii[i][j][l][s];
1984  Yniii[j][s][l][i] = Yniii[i][j][l][s];
1985 
1986  Yniii[l][i][j][s] = Yniii[i][j][l][s];
1987  Yniii[l][i][s][j] = Yniii[i][j][l][s];
1988  Yniii[l][j][i][s] = Yniii[i][j][l][s];
1989  Yniii[l][j][s][i] = Yniii[i][j][l][s];
1990  Yniii[l][s][i][j] = Yniii[i][j][l][s];
1991  Yniii[l][s][j][i] = Yniii[i][j][l][s];
1992 
1993  Yniii[s][i][j][l] = Yniii[i][j][l][s];
1994  Yniii[s][i][l][j] = Yniii[i][j][l][s];
1995  Yniii[s][j][l][i] = Yniii[i][j][l][s];
1996  Yniii[s][j][i][l] = Yniii[i][j][l][s];
1997  Yniii[s][l][i][j] = Yniii[i][j][l][s];
1998  Yniii[s][l][j][i] = Yniii[i][j][l][s];
1999  }
2000 
2001  for (k = 0; k < aq.nnii_qua_param_num; k++)
2002  {
2003  for (m = 0; m < aq.species_num; m ++)
2004  {
2005  if (strcmp (aq.nnii_qua_params[k].species_s1, aq.aqueous_species[m].species_symbol) == 0)
2006  {
2007  i = m;
2008  }
2009  }
2010 
2011  for (m = 0; m < aq.species_num; m ++)
2012  {
2013  if (strcmp (aq.nnii_qua_params[k].species_s2, aq.aqueous_species[m].species_symbol) == 0)
2014  {
2015  j = m;
2016  }
2017  }
2018 
2019  for (m = 0; m < aq.species_num; m ++)
2020  {
2021  if (strcmp (aq.nnii_qua_params[k].species_s3, aq.aqueous_species[m].species_symbol) == 0)
2022  {
2023  l = m;
2024  }
2025  }
2026 
2027  for (m = 0; m < aq.species_num; m ++)
2028  {
2029  if (strcmp (aq.nnii_qua_params[k].species_s4, aq.aqueous_species[m].species_symbol) == 0)
2030  {
2031  s = m;
2032  }
2033  }
2034 
2035  Ynnii_a = aq.nnii_qua_params[k].Yssss[0];
2036  Ynnii_b = aq.nnii_qua_params[k].Yssss[1];
2037  Ynnii_c = aq.nnii_qua_params[k].Yssss[2];
2038  Ynnii_d = aq.nnii_qua_params[k].Yssss[3];
2039  Ynnii_e = aq.nnii_qua_params[k].Yssss[4];
2040  Ynnii_f = aq.nnii_qua_params[k].Yssss[5];
2041 
2042  double Ynnii_value = Ynnii_a + Ynnii_b * T + Ynnii_c * T * log(T) + Ynnii_d * T * T + Ynnii_e * T * T * T + Ynnii_f/T;
2043 
2044  Ynnii[i][j][l][s] = Ynnii_value;
2045  Ynnii[i][j][s][l] = Ynnii[i][j][l][s];
2046  Ynnii[i][l][j][s] = Ynnii[i][j][l][s];
2047  Ynnii[i][l][s][j] = Ynnii[i][j][l][s];
2048  Ynnii[i][s][j][l] = Ynnii[i][j][l][s];
2049  Ynnii[i][s][l][j] = Ynnii[i][j][l][s];
2050 
2051  Ynnii[j][i][l][s] = Ynnii[i][j][l][s];
2052  Ynnii[j][i][s][l] = Ynnii[i][j][l][s];
2053  Ynnii[j][l][i][s] = Ynnii[i][j][l][s];
2054  Ynnii[j][l][s][i] = Ynnii[i][j][l][s];
2055  Ynnii[j][s][i][l] = Ynnii[i][j][l][s];
2056  Ynnii[j][s][l][i] = Ynnii[i][j][l][s];
2057 
2058  Ynnii[l][i][j][s] = Ynnii[i][j][l][s];
2059  Ynnii[l][i][s][j] = Ynnii[i][j][l][s];
2060  Ynnii[l][j][i][s] = Ynnii[i][j][l][s];
2061  Ynnii[l][j][s][i] = Ynnii[i][j][l][s];
2062  Ynnii[l][s][i][j] = Ynnii[i][j][l][s];
2063  Ynnii[l][s][j][i] = Ynnii[i][j][l][s];
2064 
2065  Ynnii[s][i][j][l] = Ynnii[i][j][l][s];
2066  Ynnii[s][i][l][j] = Ynnii[i][j][l][s];
2067  Ynnii[s][j][l][i] = Ynnii[i][j][l][s];
2068  Ynnii[s][j][i][l] = Ynnii[i][j][l][s];
2069  Ynnii[s][l][i][j] = Ynnii[i][j][l][s];
2070  Ynnii[s][l][j][i] = Ynnii[i][j][l][s];
2071  }
2072 
2073  for (i = 0; i < aq.species_num; i ++)
2074  {
2075  E[i] = x[i] * fabs (Z[i]);
2076  double sum = 0;
2077 
2078  for (k = 1; k < aq.species_num; k++)
2079  {
2080  if (Z[i] * Z[k] > 0) {sum = sum + x[k] * fabs (Z[k]);}
2081  }
2082 
2083  E[i] = E[i] / sum;
2084  }
2085 
2086  double sumc = 0, suma = 0;
2087  for (i = 0; i < aq.species_num; i ++)
2088  {
2089  if (Z[i] > 0) {sumc += x[i] * fabs (Z[i]);}
2090  if (Z[i] < 0) {suma += x[i] * fabs (Z[i]);}
2091  }
2092 
2093  for (i = 0; i < aq.species_num; i ++)
2094  {
2095  for (j = 0; j < aq.species_num; j ++)
2096  {
2097  if (Z[i] * Z[j] > 0)
2098  {
2099  if (i == j && Z[i] > 0) {iEj[i][j] = fabs (Z[i]) / (sumc) * (1.0 - E[i]);}
2100  if (i != j && Z[i] > 0) {iEj[i][j] = -fabs (Z[i]) * E[j] / (sumc);}
2101  if (i == j && Z[i] < 0) {iEj[i][j] = fabs (Z[i]) / (suma) * (1.0 - E[i]);}
2102  if (i != j && Z[i] < 0) {iEj[i][j] = -fabs (Z[i]) * E[j] / (suma);}
2103  }
2104  }
2105  }
2106 
2107  // DH term
2108  lnfi = -Z[index] * Z[index] * Ax * ((2.0 / ROU) * log (1.0 + ROU * sqrt (Ix)) +
2109  sqrt (Ix) * (1.0 - 2.0 * Ix / (Z[index] * Z[index])) / (1.0 + ROU * sqrt (Ix)));
2110 
2111  // cation-anion binary interaction
2112  for (i = 0; i < aq.species_num; i ++)
2113  {
2114  if (Z[i]*Z[index] != 0)
2115  {
2116  lnfi = lnfi + (x[i] * Bii[i][index] * g (aii[i][index] * sqrt (Ix)) + x[i] * B1ii[i][index] * g (a1ii[i][index] * sqrt (Ix)));
2117  }
2118  }
2119 
2120  for (i = 0; i < aq.species_num; i ++)
2121  {
2122  for (j = 0; j < aq.species_num; j ++)
2123  {
2124  if (Z[i] * Z[j] < 0)
2125  {
2126  lnfi = lnfi - 0.5 * (x[i] * x[j] * Bii[i][j] * (Z[index] * Z[index] * g (aii[i][j] * sqrt (Ix)) / (2.0 * Ix) +
2127  (1.0 - Z[index] * Z[index] / (2.0 * Ix)) * exp (-aii[i][j] * sqrt (Ix))) +
2128  x[i] * x[j] * B1ii[i][j] * (Z[index] * Z[index] * g (a1ii[i][j] * sqrt (Ix)) / (2.0 * Ix) +
2129  (1.0 - Z[index] * Z[index] / (2.0 * Ix)) * exp (-a1ii[i][j] * sqrt (Ix))));
2130  }
2131  }
2132  }
2133 
2134  // cation-cation or anion-anion binary interaction HOE term
2135  for (i = 0; i < aq.species_num; i ++)
2136  {
2137  if (Z[i] * Z[index] > 0 && i != index && Ix != 0)
2138  {
2139  double xij = 6.0 * Z[i] * Z[index] * Ax * sqrt (Ix);
2140  double xii = 6.0 * Z[i] * Z[i] * Ax * sqrt (Ix);
2141  double xjj = 6.0 * Z[index] * Z[index] * Ax * sqrt (Ix);
2142 
2143  double vij = Z[i] * Z[index] / (4.0 * Ix) * (jj(xij) - 0.5 * jj(xii) - 0.5 * jj(xjj));
2144  double vpij = - vij/Ix + (Z[i]*Z[index]/(8.0 * Ix * Ix)) * (xij * jp(xij) - 0.5 * xii * jp(xii) - 0.5 * xjj * jp(xjj));
2145 
2146  lnfi = lnfi + (2.0 * x[i] * (vij - x[index] * (vij + vpij * (Ix - Z[index] * Z[index] / 2.0))));
2147  }
2148  }
2149 
2150  for (i = 0; i < aq.species_num; i ++)
2151  {
2152  for (j = 0; j < aq.species_num; j ++)
2153  {
2154  if (Z[i] * Z[j] > 0 && i != j && i != index && j != index && Ix != 0)
2155  {
2156  double xij = 6.0 * Z[i] * Z[j] * Ax * sqrt (Ix);
2157  double xii = 6.0 * Z[i] * Z[i] * Ax * sqrt (Ix);
2158  double xjj = 6.0 * Z[j] * Z[j] * Ax * sqrt (Ix);
2159 
2160  double vij = Z[i] * Z[j] / (4.0 * Ix) * (jj(xij) - 0.5 * jj(xii) - 0.5 * jj(xjj));
2161  double vpij = - vij/Ix + (Z[i]*Z[j]/(8.0 * Ix * Ix)) * (xij * jp(xij) - 0.5 * xii * jp(xii) - 0.5 * xjj * jp(xjj));
2162 
2163  lnfi = lnfi - 0.5 * (2.0 * x[i] * x[j] * (vij + vpij * (Ix - Z[index] * Z[index] / 2.0)));
2164  }
2165  }
2166  }
2167 
2168  // neutral-cation-anion ternary interaction
2169  for (k = 0; k < aq.species_num; k ++)
2170  {
2171  double sum1 = 0;
2172 
2173  for (i = 0; i < aq.species_num; i ++)
2174  {
2175  if (Z[i] * Z[index] < 0 && Z[k] == 0)
2176  {
2177  double sum = 0;
2178 
2179  for (j = 0; j < aq.species_num; j ++)
2180  {
2181  if (Z[j] != 0)
2182  {
2183  sum += E[j] * (fabs (Z[i]) + fabs (Z[j])) / (fabs (Z[i]) * fabs (Z[j])) * Wnii[k][i][j];
2184  }
2185  }
2186 
2187  sum1 = sum1 + E[i] * ((fabs (Z[index]) + fabs (Z[i])) / fabs (Z[i]) * Wnii[k][i][index] - (fabs (Z[index]) / 2.0 + 1.0 / F) * sum);
2188  }
2189  }
2190 
2191  if (Z[k] == 0)
2192  {
2193  lnfi = lnfi + x[k] * sum1;
2194  }
2195  }
2196 
2197  for (k = 0; k < aq.species_num; k ++)
2198  {
2199  double sum1 = 0;
2200 
2201  for (i = 0; i < aq.species_num; i ++)
2202  {
2203  if (Z[i] * Z[index] < 0 && Z[k] == 0)
2204  {
2205  double sum = 0;
2206 
2207  for (j = 0; j < aq.species_num; j ++)
2208  {
2209  if (Z[j] != 0)
2210  {
2211  sum += 2.0 * x[j] * pow(fabs (Z[i]) + fabs (Z[j]), 2) / (fabs (Z[i]) * fabs (Z[j])) * Unii[k][i][j];
2212  }
2213  }
2214 
2215  sum1 = sum1 + x[i] * (pow (fabs (Z[index]) + fabs (Z[i]), 2) / (fabs (Z[i]) * fabs (Z[index])) * Unii[k][i][index] - sum);
2216  }
2217  }
2218 
2219  if (Z[k] == 0)
2220  {
2221  lnfi = lnfi + x[k] * sum1;
2222  }
2223  }
2224 
2225  for (k = 0; k < aq.species_num; k ++)
2226  {
2227  double sum1 = 0;
2228 
2229  for (i = 0; i < aq.species_num; i ++)
2230  {
2231  if (Z[i] * Z[index] < 0 && Z[k] == 0)
2232  {
2233  double sum = 0;
2234 
2235  for (j = 0; j < aq.species_num; j ++)
2236  {
2237  if (Z[j] != 0)
2238  {
2239  sum += 3.0 * x[j] * Vnii[k][i][j];
2240  }
2241  }
2242 
2243  sum1 = sum1 + x[i] * (Vnii[k][i][index] - sum);
2244  }
2245  }
2246 
2247  if (Z[k] == 0)
2248  {
2249  lnfi = lnfi + 4.0 * x[k] * x[k] * sum1;
2250  }
2251  }
2252 
2253  for (i = 0; i < aq.species_num; i ++)
2254  {
2255  for (k = 0; k < aq.species_num; k ++)
2256  {
2257  if (Z[k] == 0 && strcmp (aq.aqueous_species[k].species_symbol, WATER_SYMBOL) == 0)
2258  {
2259  break;
2260  }
2261  }
2262 
2263  if (Z[i] * Z[index] < 0)
2264  {
2265  double sum = 0;
2266 
2267  for (j = 0; j < aq.species_num; j ++)
2268  {
2269  if (j != index)
2270  {
2271  if (Z[j] != 0)
2272  {
2273  sum += E[j] * (fabs (Z[i]) + fabs (Z[j])) / (fabs (Z[i]) * fabs (Z[j])) * Wnii[k][i][j];
2274  }
2275  }
2276  }
2277  lnfi = lnfi - E[i] * ((1.0 - E[index] / 2.0) * (fabs (Z[index]) + fabs (Z[i])) / (fabs (Z[i])) * Wnii[k][i][index] - (fabs(Z[index]) / 2.0) * sum);
2278  }
2279  }
2280 
2281  // cation-anion-anion or cation-cation-anion ternary interacton
2282  for (i = 1; i < aq.species_num; i ++)
2283  {
2284  double sum = 0;
2285 
2286  for (j = 0; j < aq.species_num; j ++)
2287  {
2288  if (Z[j] * Z[index] > 0) {sum += x[j] * Wiii[index][i][j];}
2289  }
2290 
2291  for (j = 0; j < aq.species_num; j ++)
2292  {
2293  for (l = 0; l < aq.species_num; l ++)
2294  {
2295  if ((Z[j] * Z[l] > 0) && (Z[j] * Z[i] < 0))
2296  {
2297  sum = sum - 0.5 * x[j] * x[l] * Wiii[i][j][l];
2298  }
2299  }
2300  }
2301 
2302  if (Z[i] * Z[index] < 0)
2303  {
2304  lnfi = lnfi + 2.0 * E[i] * sum;
2305  }
2306  }
2307 
2308  for (i = 0; i < aq.species_num; i ++)
2309  {
2310  double sum = 0;
2311 
2312  for (j = 0; j < aq.species_num; j ++)
2313  {
2314  for (l = 0; l < aq.species_num; l ++)
2315  {
2316  if ((Z[j] * Z[l] > 0) && (Z[j] * Z[i] < 0))
2317  {
2318  sum = sum + 0.5 * x[j] * x[l] * Wiii[i][j][l];
2319  }
2320  }
2321  }
2322 
2323  if (Z[i] * Z[index] > 0)
2324  {
2325  lnfi = lnfi - 2.0 * (E[i] - iEj[index][i]) * sum;
2326  }
2327  }
2328 
2329  // neutral-ion-ion-ion quaternary interaction
2330  for (k = 0; k < aq.species_num; k ++)
2331  {
2332  double sum1 = 0;
2333 
2334  for (i = 0; i < aq.species_num; i ++)
2335  {
2336  double sum = 0;
2337 
2338  for (j = 0; j < aq.species_num; j ++)
2339  {
2340  if (Z[j] * Z[index] > 0 && Z[i] * Z[index] < 0 && j != index)
2341  {
2342  sum += x[j] * Qniii[k][index][i][j];
2343  }
2344  }
2345 
2346  for (j = 0; j < aq.species_num; j ++)
2347  {
2348  for (l = 0; l < aq.species_num; l ++)
2349  {
2350  if ((Z[j] * Z[l] > 0) && (Z[j] * Z[i] < 0))
2351  {
2352  sum = sum - 0.5 * 2.0 * x[j] * x[l] * Qniii[k][i][j][l];
2353  }
2354  }
2355  }
2356 
2357  if (Z[i] * Z[index] < 0)
2358  {
2359  sum1 = sum1 + E[i] * sum;
2360  }
2361  }
2362 
2363  if (Z[k] == 0)
2364  {
2365  lnfi = lnfi + 4.0 * x[k] * sum1;
2366  }
2367  }
2368 
2369  for (k = 0; k < aq.species_num; k ++)
2370  {
2371  double sum1 = 0;
2372 
2373  for (i = 0; i < aq.species_num; i ++)
2374  {
2375  double sum = 0;
2376 
2377  for (j = 0; j < aq.species_num; j ++)
2378  {
2379  for (l = 0; l < aq.species_num; l ++)
2380  {
2381  if ((Z[j] * Z[l] > 0) && (Z[j] * Z[i] < 0))
2382  {
2383  sum = sum + 0.5 * x[j] * x[l] * Qniii[k][i][j][l];
2384  }
2385  }
2386  }
2387 
2388  if (Z[i] * Z[index] > 0)
2389  {
2390  sum1 = sum1 + (2.0 * E[i] - iEj[index][i]) * sum;
2391  }
2392  }
2393 
2394  if (Z[k] == 0)
2395  {
2396  lnfi = lnfi - 4.0 * x[k] * sum1;
2397  }
2398  }
2399 
2401  for (k = 0; k < aq.species_num; k ++)
2402  {
2403  double sum1 = 0;
2404 
2405  for (i = 0; i < aq.species_num; i ++)
2406  {
2407  double sum = 0;
2408 
2409  for (j = 0; j < aq.species_num; j ++)
2410  {
2411  if (Z[j] * Z[index] > 0 && j != index)
2412  {
2413  sum += x[j] * Yniii[k][index][i][j];
2414  }
2415  }
2416 
2417  for (j = 1; j < aq.species_num; j ++)
2418  {
2419  for (l = 1; l < aq.species_num; l ++)
2420  {
2421  if ((Z[j] * Z[l] > 0) && (Z[j] * Z[i] < 0))
2422  {
2423  sum = sum - 0.5 * 3.0 * x[j] * x[l] * Yniii[k][i][j][l];
2424  }
2425  }
2426  }
2427 
2428  if (Z[i] * Z[index] < 0)
2429  {
2430  sum1 = sum1 + E[i] * sum;
2431  }
2432  }
2433 
2434  if (Z[k] == 0)
2435  {
2436  lnfi = lnfi + 4.0 * x[k] * x[k] * sum1;
2437  }
2438  }
2439 
2440  for (k = 0; k < aq.species_num; k ++)
2441  {
2442  double sum1 = 0;
2443 
2444  for (i = 0; i < aq.species_num; i ++)
2445  {
2446  double sum = 0;
2447 
2448  for (j = 0; j < aq.species_num; j ++)
2449  {
2450  for (l = 0; l < aq.species_num; l ++)
2451  {
2452  if ((Z[j] * Z[l] > 0) && (Z[j] * Z[i] < 0))
2453  {
2454  sum = sum + 0.5 * x[j] * x[l] * Yniii[k][i][j][l];
2455  }
2456  }
2457  }
2458 
2459  if (Z[i] * Z[index] > 0)
2460  {
2461  sum1 = sum1 + (3.0 * E[i] - iEj[index][i]) * sum;
2462  }
2463  }
2464 
2465  if (Z[k] == 0)
2466  {
2467  lnfi = lnfi - 4.0 * x[k] * x[k] * sum1;
2468  }
2469  }
2471 
2472  // neutral-neutral-cation-anion quaterary interaction
2473  for (i = 0; i < aq.species_num; i ++)
2474  {
2475  for (j = 0; j < aq.species_num; j ++)
2476  {
2477  double sum1 = 0;
2478  for (k = 0; k < aq.species_num; k ++)
2479  {
2480  double sum = 0;
2481  for (l = 0; l < aq.species_num; l ++)
2482  {
2483  if (Z[k] * Z[l] < 0 && Z[i] == 0 && Z[j] == 0)
2484  {
2485  sum += E[k] * (fabs (Z[l]) + fabs (Z[k])) / (fabs (Z[l]) * fabs (Z[k])) * Ynnii[i][j][k][l];
2486  }
2487  }
2488 
2489  if (Z[k] * Z[index] < 0)
2490  {
2491  sum1 += E[k] * ((fabs (Z[index]) + fabs (Z[k])) / fabs (Z[k]) * Ynnii[index][i][j][k] - (fabs (Z[index])/2 + 2/F) * sum);
2492  }
2493  }
2494 
2495  if (Z[i] == 0 && Z[j] == 0 && i != j)
2496  {
2497  lnfi = lnfi + 0.5 * x[i] * x[j] * sum1;
2498  }
2499  }
2500  }
2501 
2502  // neutral-neutral binary interaction
2503  for (i = 0; i < aq.species_num; i ++)
2504  {
2505  for (j = 0; j < i; j ++)
2506  {
2507  if (Z[i] == 0 && Z[j] == 0 && i != j)
2508  {
2509  lnfi = lnfi - 2.0 * x[i] * x[j] * (Wnn[i][j] + 2.0 * (x[i] - x[j]) * Unn[i][j]);
2510  }
2511  }
2512  }
2513 
2514  for (k = 0; k < aq.species_num; k ++)
2515  {
2516  if (Z[k] == 0 && strcmp (aq.aqueous_species[k].species_symbol, WATER_SYMBOL) == 0)
2517  {
2518  break;
2519  }
2520  }
2521 
2522  return exp(lnfi) * x[k];
2523 }
2524 
2525 static double psc_ai (AQUEOUS_PHASE aq, int index, double *n, double T, double P)
2526 {
2527  if (aq.aqueous_species[index].charge == 0)
2528  {
2529  printf ("The species is not a ion species!!!\n");
2530  return 1.0;
2531  }
2532 
2533  if (index >= aq.species_num)
2534  {
2535  printf ("The species dose not existed!!!\n");
2536  return 1.0;
2537  }
2538 
2539  double lnfi;
2540  double x[aq.species_num];
2541  double Z[aq.species_num];
2542  double E[aq.species_num];
2543  double iEj[aq.species_num][aq.species_num];
2544  int i, j, k, l, m, s;
2545 
2546  double sumn;
2547  double Ax, Ix, F;
2548 
2549  double aii[aq.species_num][aq.species_num];
2550  double a1ii[aq.species_num][aq.species_num];
2551  double Bii[aq.species_num][aq.species_num];
2552  double B1ii[aq.species_num][aq.species_num];
2553  double Wnn[aq.species_num][aq.species_num];
2554  double Unn[aq.species_num][aq.species_num];
2555  double Wnii[aq.species_num][aq.species_num][aq.species_num];
2556  double Unii[aq.species_num][aq.species_num][aq.species_num];
2557  double Vnii[aq.species_num][aq.species_num][aq.species_num];
2558  double Wiii[aq.species_num][aq.species_num][aq.species_num];
2559  double Qniii[aq.species_num][aq.species_num][aq.species_num][aq.species_num];
2560  double Yniii[aq.species_num][aq.species_num][aq.species_num][aq.species_num];
2561  double Ynnii[aq.species_num][aq.species_num][aq.species_num][aq.species_num];
2562 
2563  double Bii_a, Bii_b, Bii_c, Bii_d, Bii_e, Bii_f, \
2564  B1ii_a, B1ii_b, B1ii_c, B1ii_d, B1ii_e, B1ii_f, \
2565  Wnn_a, Wnn_b, Wnn_c, Wnn_d, Wnn_e, Wnn_f, \
2566  Unn_a, Unn_b, Unn_c, Unn_d, Unn_e, Unn_f, \
2567  Wnii_a, Wnii_b, Wnii_c, Wnii_d, Wnii_e, Wnii_f, \
2568  Unii_a, Unii_b, Unii_c, Unii_d, Unii_e, Unii_f, \
2569  Vnii_a, Vnii_b, Vnii_c, Vnii_d, Vnii_e, Vnii_f, \
2570  Wiii_a, Wiii_b, Wiii_c, Wiii_d, Wiii_e, Wiii_f, \
2571  Qniii_a, Qniii_b, Qniii_c, Qniii_d, Qniii_e, Qniii_f, \
2572  Yniii_a, Yniii_b, Yniii_c, Yniii_d, Yniii_e, Yniii_f, \
2573  Ynnii_a, Ynnii_b, Ynnii_c, Ynnii_d, Ynnii_e, Ynnii_f;
2574 
2575  Ax = ax (T, 1.0);
2576 
2577  for (k = 0; k < aq.species_num; k ++)
2578  {
2579  Z[k] = aq.aqueous_species[k].charge;
2580  }
2581 
2582  sumn = 0;
2583  for (k = 0; k < aq.species_num; k ++)
2584  {
2585  sumn += n[k];
2586  }
2587 
2588  for (k = 0; k < aq.species_num; k ++)
2589  {
2590  x[k] = n[k] / sumn;
2591  }
2592 
2593  F = 0;
2594  for (k = 0; k < aq.species_num; k++)
2595  {
2596  F += 0.5 * x[k] * fabs(Z[k]);
2597  }
2598  F = 1.0/F;
2599 
2600  Ix = 0;
2601  for (k = 0; k < aq.species_num; k ++)
2602  {
2603  Ix += 0.5 * (x[k] * Z[k] * Z[k]);
2604  }
2605 
2606  for (k = 0; k < aq.species_num; k ++)
2607  {
2608  E[k] = 0;
2609  }
2610 
2611  // INITALIZATION set all binary parameters zero
2612  for (i = 0; i < aq.species_num; i ++)
2613  {
2614  for (j = 0; j < aq.species_num; j ++)
2615  {
2616  aii[i][j] = 0;
2617  a1ii[i][j] = 0;
2618  Bii[i][j] = 0;
2619  B1ii[i][j] = 0;
2620  Wnn[i][j] = 0;
2621  Unn[i][j] = 0;
2622  }
2623  }
2624 
2625  // INITALIZATION set all ternary parameters zero
2626  for (i = 0; i < aq.species_num; i ++)
2627  {
2628  for (j = 0; j < aq.species_num; j ++)
2629  {
2630  for (l = 0; l < aq.species_num; l ++)
2631  {
2632  Wnii[i][j][l] = 0;
2633  Unii[i][j][l] = 0;
2634  Vnii[i][j][l] = 0;
2635  Wiii[i][j][l] = 0;
2636  }
2637  }
2638  }
2639 
2640  // INITALIZATION set all quaternary parameters zero
2641  for (i = 0; i < aq.species_num; i ++)
2642  {
2643  for (j = 0; j < aq.species_num; j ++)
2644  {
2645  for (l = 0; l < aq.species_num; l ++)
2646  {
2647  for (m = 0; m < aq.species_num; m ++)
2648  {
2649  Qniii[i][j][l][m] = 0;
2650  Yniii[i][j][l][m] = 0;
2651  Ynnii[i][j][l][m] = 0;
2652  }
2653  }
2654  }
2655  }
2656 
2657  for (k = 0; k < aq.ii_bin_param_num; k ++)
2658  {
2659  for (m = 0; m < aq.species_num; m ++)
2660  {
2661  if (strcmp (aq.ii_bin_params[k].species_s1, aq.aqueous_species[m].species_symbol) == 0)
2662  {
2663  i = m;
2664  }
2665  }
2666 
2667  for (m = 0; m < aq.species_num; m ++)
2668  {
2669  if (strcmp (aq.ii_bin_params[k].species_s2, aq.aqueous_species[m].species_symbol) == 0)
2670  {
2671  j = m;
2672  }
2673  }
2674 
2675  double aii_value = aq.ii_bin_params[k].alpha;
2676  double a1ii_value = aq.ii_bin_params[k].alpha1;
2677 
2678  Bii_a = aq.ii_bin_params[k].B0ss[0];
2679  Bii_b = aq.ii_bin_params[k].B0ss[1];
2680  Bii_c = aq.ii_bin_params[k].B0ss[2];
2681  Bii_d = aq.ii_bin_params[k].B0ss[3];
2682  Bii_e = aq.ii_bin_params[k].B0ss[4];
2683  Bii_f = aq.ii_bin_params[k].B0ss[5];
2684  B1ii_a = aq.ii_bin_params[k].B1ss[0];
2685  B1ii_b = aq.ii_bin_params[k].B1ss[1];
2686  B1ii_c = aq.ii_bin_params[k].B1ss[2];
2687  B1ii_d = aq.ii_bin_params[k].B1ss[3];
2688  B1ii_e = aq.ii_bin_params[k].B1ss[4];
2689  B1ii_f = aq.ii_bin_params[k].B1ss[5];
2690 
2691  double Bii_value = Bii_a + Bii_b * T + Bii_c * T * log(T) + Bii_d * T * T + Bii_e * T * T * T + Bii_f/T;
2692  double B1ii_value = B1ii_a + B1ii_b * T + B1ii_c * T * log(T) + B1ii_d * T * T + B1ii_e * T * T * T + B1ii_f/T;
2693 
2694  aii[i][j] = aii_value; aii[j][i] = aii[i][j];
2695  a1ii[i][j] = a1ii_value; a1ii[j][i] = a1ii[i][j];
2696  Bii[i][j] = Bii_value; Bii[j][i] = Bii[i][j];
2697  B1ii[i][j] = B1ii_value; B1ii[j][i] = B1ii[i][j];
2698  }
2699 
2700  for (k = 0; k < aq.nn_bin_param_num; k ++)
2701  {
2702  for (m = 0; m < aq.species_num; m ++)
2703  {
2704  if (strcmp (aq.nn_bin_params[k].species_s1, aq.aqueous_species[m].species_symbol) == 0)
2705  {
2706  i = m;
2707  }
2708  }
2709 
2710  for (m = 0; m < aq.species_num; m ++)
2711  {
2712  if (strcmp (aq.nn_bin_params[k].species_s2, aq.aqueous_species[m].species_symbol) == 0)
2713  {
2714  j = m;
2715  }
2716  }
2717 
2718  Wnn_a = aq.nn_bin_params[k].Wss[0];
2719  Wnn_b = aq.nn_bin_params[k].Wss[1];
2720  Wnn_c = aq.nn_bin_params[k].Wss[2];
2721  Wnn_d = aq.nn_bin_params[k].Wss[3];
2722  Wnn_e = aq.nn_bin_params[k].Wss[4];
2723  Wnn_f = aq.nn_bin_params[k].Wss[5];
2724  Unn_a = aq.nn_bin_params[k].Uss[0];
2725  Unn_b = aq.nn_bin_params[k].Uss[1];
2726  Unn_c = aq.nn_bin_params[k].Uss[2];
2727  Unn_d = aq.nn_bin_params[k].Uss[3];
2728  Unn_e = aq.nn_bin_params[k].Uss[4];
2729  Unn_f = aq.nn_bin_params[k].Uss[5];
2730 
2731  double Wnn_value = Wnn_a + Wnn_b * T + Wnn_c * T * log(T) + Wnn_d * T * T + Wnn_e * T * T * T + Wnn_f/T;
2732  double Unn_value = Unn_a + Unn_b * T + Unn_c * T * log(T) + Unn_d * T * T + Unn_e * T * T * T + Unn_f/T;
2733 
2734  Wnn[i][j] = Wnn_value; Wnn[j][i] = Wnn[i][j];
2735  Unn[i][j] = Unn_value; Unn[j][i] = Unn[i][j];
2736  }
2737 
2738  for (k = 0; k < aq.nii_ter_param_num; k++)
2739  {
2740  for (m = 0; m < aq.species_num; m ++)
2741  {
2742  if (strcmp (aq.nii_ter_params[k].species_s1, aq.aqueous_species[m].species_symbol) == 0)
2743  {
2744  i = m;
2745  }
2746  }
2747 
2748  for (m = 0; m < aq.species_num; m ++)
2749  {
2750  if (strcmp (aq.nii_ter_params[k].species_s2, aq.aqueous_species[m].species_symbol) == 0)
2751  {
2752  j = m;
2753  }
2754  }
2755 
2756  for (m = 0; m < aq.species_num; m ++)
2757  {
2758  if (strcmp (aq.nii_ter_params[k].species_s3, aq.aqueous_species[m].species_symbol) == 0)
2759  {
2760  l = m;
2761  }
2762  }
2763 
2764  Wnii_a = aq.nii_ter_params[k].Wsss[0];
2765  Wnii_b = aq.nii_ter_params[k].Wsss[1];
2766  Wnii_c = aq.nii_ter_params[k].Wsss[2];
2767  Wnii_d = aq.nii_ter_params[k].Wsss[3];
2768  Wnii_e = aq.nii_ter_params[k].Wsss[4];
2769  Wnii_f = aq.nii_ter_params[k].Wsss[5];
2770  Unii_a = aq.nii_ter_params[k].Usss[0];
2771  Unii_b = aq.nii_ter_params[k].Usss[1];
2772  Unii_c = aq.nii_ter_params[k].Usss[2];
2773  Unii_d = aq.nii_ter_params[k].Usss[3];
2774  Unii_e = aq.nii_ter_params[k].Usss[4];
2775  Unii_f = aq.nii_ter_params[k].Usss[5];
2776  Vnii_a = aq.nii_ter_params[k].Vsss[0];
2777  Vnii_b = aq.nii_ter_params[k].Vsss[1];
2778  Vnii_c = aq.nii_ter_params[k].Vsss[2];
2779  Vnii_d = aq.nii_ter_params[k].Vsss[3];
2780  Vnii_e = aq.nii_ter_params[k].Vsss[4];
2781  Vnii_f = aq.nii_ter_params[k].Vsss[5];
2782 
2783  double Wnii_value = Wnii_a + Wnii_b * T + Wnii_c * T * log(T) + Wnii_d * T * T + Wnii_e * T * T * T + Wnii_f/T;
2784  double Unii_value = Unii_a + Unii_b * T + Unii_c * T * log(T) + Unii_d * T * T + Unii_e * T * T * T + Unii_f/T;
2785  double Vnii_value = Vnii_a + Vnii_b * T + Vnii_c * T * log(T) + Vnii_d * T * T + Vnii_e * T * T * T + Vnii_f/T;
2786 
2787  Wnii[i][j][l] = Wnii_value;
2788  Wnii[i][l][j] = Wnii[i][j][l];
2789  Wnii[j][i][l] = Wnii[i][j][l];
2790  Wnii[j][l][i] = Wnii[i][j][l];
2791  Wnii[l][i][j] = Wnii[i][j][l];
2792  Wnii[l][j][i] = Wnii[i][j][l];
2793 
2794  Unii[i][j][l] = Unii_value;
2795  Unii[i][l][j] = Unii[i][j][l];
2796  Unii[j][i][l] = Unii[i][j][l];
2797  Unii[j][l][i] = Unii[i][j][l];
2798  Unii[l][i][j] = Unii[i][j][l];
2799  Unii[l][j][i] = Unii[i][j][l];
2800 
2801  Vnii[i][j][l] = Vnii_value;
2802  Vnii[i][l][j] = Vnii[i][j][l];
2803  Vnii[j][i][l] = Vnii[i][j][l];
2804  Vnii[j][l][i] = Vnii[i][j][l];
2805  Vnii[l][i][j] = Vnii[i][j][l];
2806  Vnii[l][j][i] = Vnii[i][j][l];
2807  }
2808 
2809  for (k = 0; k < aq.iii_ter_param_num; k++)
2810  {
2811  for (m = 0; m < aq.species_num; m ++)
2812  {
2813  if (strcmp (aq.iii_ter_params[k].species_s1, aq.aqueous_species[m].species_symbol) == 0)
2814  {
2815  i = m;
2816  }
2817  }
2818 
2819  for (m = 0; m < aq.species_num; m ++)
2820  {
2821  if (strcmp (aq.iii_ter_params[k].species_s2, aq.aqueous_species[m].species_symbol) == 0)
2822  {
2823  j = m;
2824  }
2825  }
2826 
2827  for (m = 0; m < aq.species_num; m ++)
2828  {
2829  if (strcmp (aq.iii_ter_params[k].species_s3, aq.aqueous_species[m].species_symbol) == 0)
2830  {
2831  l = m;
2832  }
2833  }
2834 
2835  Wiii_a = aq.iii_ter_params[k].Wsss[0];
2836  Wiii_b = aq.iii_ter_params[k].Wsss[1];
2837  Wiii_c = aq.iii_ter_params[k].Wsss[2];
2838  Wiii_d = aq.iii_ter_params[k].Wsss[3];
2839  Wiii_e = aq.iii_ter_params[k].Wsss[4];
2840  Wiii_f = aq.iii_ter_params[k].Wsss[5];
2841 
2842  double Wiii_value = Wiii_a + Wiii_b * T + Wiii_c * T * log(T) + Wiii_d * T * T + Wiii_e * T * T * T + Wiii_f/T;
2843 
2844  Wiii[i][j][l] = Wiii_value;
2845  Wiii[i][l][j] = Wiii[i][j][l];
2846  Wiii[j][i][l] = Wiii[i][j][l];
2847  Wiii[j][l][i] = Wiii[i][j][l];
2848  Wiii[l][i][j] = Wiii[i][j][l];
2849  Wiii[l][j][i] = Wiii[i][j][l];
2850  }
2851 
2852  for (k = 0; k < aq.niii_qua_param_num; k++)
2853  {
2854  for (m = 0; m < aq.species_num; m ++)
2855  {
2856  if (strcmp (aq.niii_qua_params[k].species_s1, aq.aqueous_species[m].species_symbol) == 0)
2857  {
2858  i = m;
2859  }
2860  }
2861 
2862  for (m = 0; m < aq.species_num; m ++)
2863  {
2864  if (strcmp (aq.niii_qua_params[k].species_s2, aq.aqueous_species[m].species_symbol) == 0)
2865  {
2866  j = m;
2867  }
2868  }
2869 
2870  for (m = 0; m < aq.species_num; m ++)
2871  {
2872  if (strcmp (aq.niii_qua_params[k].species_s3, aq.aqueous_species[m].species_symbol) == 0)
2873  {
2874  l = m;
2875  }
2876  }
2877 
2878  for (m = 0; m < aq.species_num; m ++)
2879  {
2880  if (strcmp (aq.niii_qua_params[k].species_s4, aq.aqueous_species[m].species_symbol) == 0)
2881  {
2882  s = m;
2883  }
2884  }
2885  Qniii_a = aq.niii_qua_params[k].Qssss[0];
2886  Qniii_b = aq.niii_qua_params[k].Qssss[1];
2887  Qniii_c = aq.niii_qua_params[k].Qssss[2];
2888  Qniii_d = aq.niii_qua_params[k].Qssss[3];
2889  Qniii_e = aq.niii_qua_params[k].Qssss[4];
2890  Qniii_f = aq.niii_qua_params[k].Qssss[5];
2891  Yniii_a = aq.niii_qua_params[k].Yssss[0];
2892  Yniii_b = aq.niii_qua_params[k].Yssss[1];
2893  Yniii_c = aq.niii_qua_params[k].Yssss[2];
2894  Yniii_d = aq.niii_qua_params[k].Yssss[3];
2895  Yniii_e = aq.niii_qua_params[k].Yssss[4];
2896  Yniii_f = aq.niii_qua_params[k].Yssss[5];
2897 
2898  double Qniii_value = Qniii_a + Qniii_b * T + Qniii_c * T * log(T) + Qniii_d * T * T + Qniii_e * T * T * T + Qniii_f/T;
2899  double Yniii_value = Yniii_a + Yniii_b * T + Yniii_c * T * log(T) + Yniii_d * T * T + Yniii_e * T * T * T + Yniii_f/T;
2900 
2901  Qniii[i][j][l][s] = Qniii_value;
2902  Qniii[i][j][s][l] = Qniii[i][j][l][s];
2903  Qniii[i][l][j][s] = Qniii[i][j][l][s];
2904  Qniii[i][l][s][j] = Qniii[i][j][l][s];
2905  Qniii[i][s][j][l] = Qniii[i][j][l][s];
2906  Qniii[i][s][l][j] = Qniii[i][j][l][s];
2907 
2908  Qniii[j][i][l][s] = Qniii[i][j][l][s];
2909  Qniii[j][i][s][l] = Qniii[i][j][l][s];
2910  Qniii[j][l][i][s] = Qniii[i][j][l][s];
2911  Qniii[j][l][s][i] = Qniii[i][j][l][s];
2912  Qniii[j][s][i][l] = Qniii[i][j][l][s];
2913  Qniii[j][s][l][i] = Qniii[i][j][l][s];
2914 
2915  Qniii[l][i][j][s] = Qniii[i][j][l][s];
2916  Qniii[l][i][s][j] = Qniii[i][j][l][s];
2917  Qniii[l][j][i][s] = Qniii[i][j][l][s];
2918  Qniii[l][j][s][i] = Qniii[i][j][l][s];
2919  Qniii[l][s][i][j] = Qniii[i][j][l][s];
2920  Qniii[l][s][j][i] = Qniii[i][j][l][s];
2921 
2922  Qniii[s][i][j][l] = Qniii[i][j][l][s];
2923  Qniii[s][i][l][j] = Qniii[i][j][l][s];
2924  Qniii[s][j][l][i] = Qniii[i][j][l][s];
2925  Qniii[s][j][i][l] = Qniii[i][j][l][s];
2926  Qniii[s][l][i][j] = Qniii[i][j][l][s];
2927  Qniii[s][l][j][i] = Qniii[i][j][l][s];
2928 
2929  Yniii[i][j][l][s] = Yniii_value;
2930  Yniii[i][j][s][l] = Yniii[i][j][l][s];
2931  Yniii[i][l][j][s] = Yniii[i][j][l][s];
2932  Yniii[i][l][s][j] = Yniii[i][j][l][s];
2933  Yniii[i][s][j][l] = Yniii[i][j][l][s];
2934  Yniii[i][s][l][j] = Yniii[i][j][l][s];
2935 
2936  Yniii[j][i][l][s] = Yniii[i][j][l][s];
2937  Yniii[j][i][s][l] = Yniii[i][j][l][s];
2938  Yniii[j][l][i][s] = Yniii[i][j][l][s];
2939  Yniii[j][l][s][i] = Yniii[i][j][l][s];
2940  Yniii[j][s][i][l] = Yniii[i][j][l][s];
2941  Yniii[j][s][l][i] = Yniii[i][j][l][s];
2942 
2943  Yniii[l][i][j][s] = Yniii[i][j][l][s];
2944  Yniii[l][i][s][j] = Yniii[i][j][l][s];
2945  Yniii[l][j][i][s] = Yniii[i][j][l][s];
2946  Yniii[l][j][s][i] = Yniii[i][j][l][s];
2947  Yniii[l][s][i][j] = Yniii[i][j][l][s];
2948  Yniii[l][s][j][i] = Yniii[i][j][l][s];
2949 
2950  Yniii[s][i][j][l] = Yniii[i][j][l][s];
2951  Yniii[s][i][l][j] = Yniii[i][j][l][s];
2952  Yniii[s][j][l][i] = Yniii[i][j][l][s];
2953  Yniii[s][j][i][l] = Yniii[i][j][l][s];
2954  Yniii[s][l][i][j] = Yniii[i][j][l][s];
2955  Yniii[s][l][j][i] = Yniii[i][j][l][s];
2956  }
2957 
2958  for (k = 0; k < aq.nnii_qua_param_num; k++)
2959  {
2960  for (m = 0; m < aq.species_num; m ++)
2961  {
2962  if (strcmp (aq.nnii_qua_params[k].species_s1, aq.aqueous_species[m].species_symbol) == 0)
2963  {
2964  i = m;
2965  }
2966  }
2967 
2968  for (m = 0; m < aq.species_num; m ++)
2969  {
2970  if (strcmp (aq.nnii_qua_params[k].species_s2, aq.aqueous_species[m].species_symbol) == 0)
2971  {
2972  j = m;
2973  }
2974  }
2975 
2976  for (m = 0; m < aq.species_num; m ++)
2977  {
2978  if (strcmp (aq.nnii_qua_params[k].species_s3, aq.aqueous_species[m].species_symbol) == 0)
2979  {
2980  l = m;
2981  }
2982  }
2983 
2984  for (m = 0; m < aq.species_num; m ++)
2985  {
2986  if (strcmp (aq.nnii_qua_params[k].species_s4, aq.aqueous_species[m].species_symbol) == 0)
2987  {
2988  s = m;
2989  }
2990  }
2991 
2992  Ynnii_a = aq.nnii_qua_params[k].Yssss[0];
2993  Ynnii_b = aq.nnii_qua_params[k].Yssss[1];
2994  Ynnii_c = aq.nnii_qua_params[k].Yssss[2];
2995  Ynnii_d = aq.nnii_qua_params[k].Yssss[3];
2996  Ynnii_e = aq.nnii_qua_params[k].Yssss[4];
2997  Ynnii_f = aq.nnii_qua_params[k].Yssss[5];
2998 
2999  double Ynnii_value = Ynnii_a + Ynnii_b * T + Ynnii_c * T * log(T) + Ynnii_d * T * T + Ynnii_e * T * T * T + Ynnii_f/T;
3000 
3001  Ynnii[i][j][l][s] = Ynnii_value;
3002  Ynnii[i][j][s][l] = Ynnii[i][j][l][s];
3003  Ynnii[i][l][j][s] = Ynnii[i][j][l][s];
3004  Ynnii[i][l][s][j] = Ynnii[i][j][l][s];
3005  Ynnii[i][s][j][l] = Ynnii[i][j][l][s];
3006  Ynnii[i][s][l][j] = Ynnii[i][j][l][s];
3007 
3008  Ynnii[j][i][l][s] = Ynnii[i][j][l][s];
3009  Ynnii[j][i][s][l] = Ynnii[i][j][l][s];
3010  Ynnii[j][l][i][s] = Ynnii[i][j][l][s];
3011  Ynnii[j][l][s][i] = Ynnii[i][j][l][s];
3012  Ynnii[j][s][i][l] = Ynnii[i][j][l][s];
3013  Ynnii[j][s][l][i] = Ynnii[i][j][l][s];
3014 
3015  Ynnii[l][i][j][s] = Ynnii[i][j][l][s];
3016  Ynnii[l][i][s][j] = Ynnii[i][j][l][s];
3017  Ynnii[l][j][i][s] = Ynnii[i][j][l][s];
3018  Ynnii[l][j][s][i] = Ynnii[i][j][l][s];
3019  Ynnii[l][s][i][j] = Ynnii[i][j][l][s];
3020  Ynnii[l][s][j][i] = Ynnii[i][j][l][s];
3021 
3022  Ynnii[s][i][j][l] = Ynnii[i][j][l][s];
3023  Ynnii[s][i][l][j] = Ynnii[i][j][l][s];
3024  Ynnii[s][j][l][i] = Ynnii[i][j][l][s];
3025  Ynnii[s][j][i][l] = Ynnii[i][j][l][s];
3026  Ynnii[s][l][i][j] = Ynnii[i][j][l][s];
3027  Ynnii[s][l][j][i] = Ynnii[i][j][l][s];
3028  }
3029 
3030  for (i = 0; i < aq.species_num; i ++)
3031  {
3032  E[i] = x[i] * fabs (Z[i]);
3033  double sum = 0;
3034 
3035  for (k = 1; k < aq.species_num; k++)
3036  {
3037  if (Z[i] * Z[k] > 0) {sum = sum + x[k] * fabs (Z[k]);}
3038  }
3039 
3040  E[i] = E[i] / sum;
3041  }
3042 
3043  double sumc = 0, suma = 0;
3044  for (i = 0; i < aq.species_num; i ++)
3045  {
3046  if (Z[i] > 0) {sumc += x[i] * fabs (Z[i]);}
3047  if (Z[i] < 0) {suma += x[i] * fabs (Z[i]);}
3048  }
3049 
3050  for (i = 0; i < aq.species_num; i ++)
3051  {
3052  for (j = 0; j < aq.species_num; j ++)
3053  {
3054  if (Z[i] * Z[j] > 0)
3055  {
3056  if (i == j && Z[i] > 0) {iEj[i][j] = fabs (Z[i]) / (sumc) * (1.0 - E[i]);}
3057  if (i != j && Z[i] > 0) {iEj[i][j] = -fabs (Z[i]) * E[j] / (sumc);}
3058  if (i == j && Z[i] < 0) {iEj[i][j] = fabs (Z[i]) / (suma) * (1.0 - E[i]);}
3059  if (i != j && Z[i] < 0) {iEj[i][j] = -fabs (Z[i]) * E[j] / (suma);}
3060  }
3061  }
3062  }
3063 
3064  // DH term
3065  lnfi = -Z[index] * Z[index] * Ax * ((2.0 / ROU) * log (1.0 + ROU * sqrt (Ix)) +
3066  sqrt (Ix) * (1.0 - 2.0 * Ix / (Z[index] * Z[index])) / (1.0 + ROU * sqrt (Ix)));
3067 
3068  // cation-anion binary interaction
3069  for (i = 0; i < aq.species_num; i ++)
3070  {
3071  if (Z[i]*Z[index] != 0)
3072  {
3073  lnfi = lnfi + (x[i] * Bii[i][index] * g (aii[i][index] * sqrt (Ix)) + x[i] * B1ii[i][index] * g (a1ii[i][index] * sqrt (Ix)));
3074  }
3075  }
3076 
3077  for (i = 0; i < aq.species_num; i ++)
3078  {
3079  for (j = 0; j < aq.species_num; j ++)
3080  {
3081  if (Z[i] * Z[j] < 0)
3082  {
3083  lnfi = lnfi - 0.5 * (x[i] * x[j] * Bii[i][j] * (Z[index] * Z[index] * g (aii[i][j] * sqrt (Ix)) / (2.0 * Ix) +
3084  (1.0 - Z[index] * Z[index] / (2.0 * Ix)) * exp (-aii[i][j] * sqrt (Ix))) +
3085  x[i] * x[j] * B1ii[i][j] * (Z[index] * Z[index] * g (a1ii[i][j] * sqrt (Ix)) / (2.0 * Ix) +
3086  (1.0 - Z[index] * Z[index] / (2.0 * Ix)) * exp (-a1ii[i][j] * sqrt (Ix))));
3087  }
3088  }
3089  }
3090 
3091  // cation-cation or anion-anion binary interaction HOE term
3092  for (i = 0; i < aq.species_num; i ++)
3093  {
3094  if (Z[i] * Z[index] > 0 && i != index && Ix != 0)
3095  {
3096  double xij = 6.0 * Z[i] * Z[index] * Ax * sqrt (Ix);
3097  double xii = 6.0 * Z[i] * Z[i] * Ax * sqrt (Ix);
3098  double xjj = 6.0 * Z[index] * Z[index] * Ax * sqrt (Ix);
3099 
3100  double vij = Z[i] * Z[index] / (4.0 * Ix) * (jj(xij) - 0.5 * jj(xii) - 0.5 * jj(xjj));
3101  double vpij = - vij/Ix + (Z[i]*Z[index]/(8.0 * Ix * Ix)) * (xij * jp(xij) - 0.5 * xii * jp(xii) - 0.5 * xjj * jp(xjj));
3102 
3103  lnfi = lnfi + (2.0 * x[i] * (vij - x[index] * (vij + vpij * (Ix - Z[index] * Z[index] / 2.0))));
3104  }
3105  }
3106 
3107  for (i = 0; i < aq.species_num; i ++)
3108  {
3109  for (j = 0; j < aq.species_num; j ++)
3110  {
3111  if (Z[i] * Z[j] > 0 && i != j && i != index && j != index && Ix != 0)
3112  {
3113  double xij = 6.0 * Z[i] * Z[j] * Ax * sqrt (Ix);
3114  double xii = 6.0 * Z[i] * Z[i] * Ax * sqrt (Ix);
3115  double xjj = 6.0 * Z[j] * Z[j] * Ax * sqrt (Ix);
3116 
3117  double vij = Z[i] * Z[j] / (4.0 * Ix) * (jj(xij) - 0.5 * jj(xii) - 0.5 * jj(xjj));
3118  double vpij = - vij/Ix + (Z[i]*Z[j]/(8.0 * Ix * Ix)) * (xij * jp(xij) - 0.5 * xii * jp(xii) - 0.5 * xjj * jp(xjj));
3119 
3120  lnfi = lnfi - 0.5 * (2.0 * x[i] * x[j] * (vij + vpij * (Ix - Z[index] * Z[index] / 2.0)));
3121  }
3122  }
3123  }
3124 
3125  // neutral-cation-anion ternary interaction
3126  for (k = 0; k < aq.species_num; k ++)
3127  {
3128  double sum1 = 0;
3129 
3130  for (i = 0; i < aq.species_num; i ++)
3131  {
3132  if (Z[i] * Z[index] < 0 && Z[k] == 0)
3133  {
3134  double sum = 0;
3135 
3136  for (j = 0; j < aq.species_num; j ++)
3137  {
3138  if (Z[j] != 0)
3139  {
3140  sum += E[j] * (fabs (Z[i]) + fabs (Z[j])) / (fabs (Z[i]) * fabs (Z[j])) * Wnii[k][i][j];
3141  }
3142  }
3143 
3144  sum1 = sum1 + E[i] * ((fabs (Z[index]) + fabs (Z[i])) / fabs (Z[i]) * Wnii[k][i][index] - (fabs (Z[index]) / 2.0 + 1.0 / F) * sum);
3145  }
3146  }
3147 
3148  if (Z[k] == 0)
3149  {
3150  lnfi = lnfi + x[k] * sum1;
3151  }
3152  }
3153 
3154  for (k = 0; k < aq.species_num; k ++)
3155  {
3156  double sum1 = 0;
3157 
3158  for (i = 0; i < aq.species_num; i ++)
3159  {
3160  if (Z[i] * Z[index] < 0 && Z[k] == 0)
3161  {
3162  double sum = 0;
3163 
3164  for (j = 0; j < aq.species_num; j ++)
3165  {
3166  if (Z[j] != 0)
3167  {
3168  sum += 2.0 * x[j] * pow(fabs (Z[i]) + fabs (Z[j]), 2) / (fabs (Z[i]) * fabs (Z[j])) * Unii[k][i][j];
3169  }
3170  }
3171 
3172  sum1 = sum1 + x[i] * (pow (fabs (Z[index]) + fabs (Z[i]), 2) / (fabs (Z[i]) * fabs (Z[index])) * Unii[k][i][index] - sum);
3173  }
3174  }
3175 
3176  if (Z[k] == 0)
3177  {
3178  lnfi = lnfi + x[k] * sum1;
3179  }
3180  }
3181 
3182  for (k = 0; k < aq.species_num; k ++)
3183  {
3184  double sum1 = 0;
3185 
3186  for (i = 0; i < aq.species_num; i ++)
3187  {
3188  if (Z[i] * Z[index] < 0 && Z[k] == 0)
3189  {
3190  double sum = 0;
3191 
3192  for (j = 0; j < aq.species_num; j ++)
3193  {
3194  if (Z[j] != 0)
3195  {
3196  sum += 3.0 * x[j] * Vnii[k][i][j];
3197  }
3198  }
3199 
3200  sum1 = sum1 + x[i] * (Vnii[k][i][index] - sum);
3201  }
3202  }
3203 
3204  if (Z[k] == 0)
3205  {
3206  lnfi = lnfi + 4.0 * x[k] * x[k] * sum1;
3207  }
3208  }
3209 
3210  for (i = 0; i < aq.species_num; i ++)
3211  {
3212  for (k = 0; k < aq.species_num; k ++)
3213  {
3214  if (Z[k] == 0 && strcmp (aq.aqueous_species[k].species_symbol, WATER_SYMBOL) == 0)
3215  {
3216  break;
3217  }
3218  }
3219 
3220  if (Z[i] * Z[index] < 0)
3221  {
3222  double sum = 0;
3223 
3224  for (j = 0; j < aq.species_num; j ++)
3225  {
3226  if (j != index)
3227  {
3228  if (Z[j] != 0)
3229  {
3230  sum += E[j] * (fabs (Z[i]) + fabs (Z[j])) / (fabs (Z[i]) * fabs (Z[j])) * Wnii[k][i][j];
3231  }
3232  }
3233  }
3234  lnfi = lnfi - E[i] * ((1.0 - E[index] / 2.0) * (fabs (Z[index]) + fabs (Z[i])) / (fabs (Z[i])) * Wnii[k][i][index] - (fabs(Z[index]) / 2.0) * sum);
3235  }
3236  }
3237 
3238  // cation-anion-anion or cation-cation-anion ternary interacton
3239  for (i = 1; i < aq.species_num; i ++)
3240  {
3241  double sum = 0;
3242 
3243  for (j = 0; j < aq.species_num; j ++)
3244  {
3245  if (Z[j] * Z[index] > 0) {sum += x[j] * Wiii[index][i][j];}
3246  }
3247 
3248  for (j = 0; j < aq.species_num; j ++)
3249  {
3250  for (l = 0; l < aq.species_num; l ++)
3251  {
3252  if ((Z[j] * Z[l] > 0) && (Z[j] * Z[i] < 0))
3253  {
3254  sum = sum - 0.5 * x[j] * x[l] * Wiii[i][j][l];
3255  }
3256  }
3257  }
3258 
3259  if (Z[i] * Z[index] < 0)
3260  {
3261  lnfi = lnfi + 2.0 * E[i] * sum;
3262  }
3263  }
3264 
3265  for (i = 0; i < aq.species_num; i ++)
3266  {
3267  double sum = 0;
3268 
3269  for (j = 0; j < aq.species_num; j ++)
3270  {
3271  for (l = 0; l < aq.species_num; l ++)
3272  {
3273  if ((Z[j] * Z[l] > 0) && (Z[j] * Z[i] < 0))
3274  {
3275  sum = sum + 0.5 * x[j] * x[l] * Wiii[i][j][l];
3276  }
3277  }
3278  }
3279 
3280  if (Z[i] * Z[index] > 0)
3281  {
3282  lnfi = lnfi - 2.0 * (E[i] - iEj[index][i]) * sum;
3283  }
3284  }
3285 
3286  // neutral-ion-ion-ion quaternary interaction
3287  for (k = 0; k < aq.species_num; k ++)
3288  {
3289  double sum1 = 0;
3290 
3291  for (i = 0; i < aq.species_num; i ++)
3292  {
3293  double sum = 0;
3294 
3295  for (j = 0; j < aq.species_num; j ++)
3296  {
3297  if (Z[j] * Z[index] > 0 && Z[i] * Z[index] < 0 && j != index)
3298  {
3299  sum += x[j] * Qniii[k][index][i][j];
3300  }
3301  }
3302 
3303  for (j = 0; j < aq.species_num; j ++)
3304  {
3305  for (l = 0; l < aq.species_num; l ++)
3306  {
3307  if ((Z[j] * Z[l] > 0) && (Z[j] * Z[i] < 0))
3308  {
3309  sum = sum - 0.5 * 2.0 * x[j] * x[l] * Qniii[k][i][j][l];
3310  }
3311  }
3312  }
3313 
3314  if (Z[i] * Z[index] < 0)
3315  {
3316  sum1 = sum1 + E[i] * sum;
3317  }
3318  }
3319 
3320  if (Z[k] == 0)
3321  {
3322  lnfi = lnfi + 4.0 * x[k] * sum1;
3323  }
3324  }
3325 
3326  for (k = 0; k < aq.species_num; k ++)
3327  {
3328  double sum1 = 0;
3329 
3330  for (i = 0; i < aq.species_num; i ++)
3331  {
3332  double sum = 0;
3333 
3334  for (j = 0; j < aq.species_num; j ++)
3335  {
3336  for (l = 0; l < aq.species_num; l ++)
3337  {
3338  if ((Z[j] * Z[l] > 0) && (Z[j] * Z[i] < 0))
3339  {
3340  sum = sum + 0.5 * x[j] * x[l] * Qniii[k][i][j][l];
3341  }
3342  }
3343  }
3344 
3345  if (Z[i] * Z[index] > 0)
3346  {
3347  sum1 = sum1 + (2.0 * E[i] - iEj[index][i]) * sum;
3348  }
3349  }
3350 
3351  if (Z[k] == 0)
3352  {
3353  lnfi = lnfi - 4.0 * x[k] * sum1;
3354  }
3355  }
3356 
3358  for (k = 0; k < aq.species_num; k ++)
3359  {
3360  double sum1 = 0;
3361 
3362  for (i = 0; i < aq.species_num; i ++)
3363  {
3364  double sum = 0;
3365 
3366  for (j = 0; j < aq.species_num; j ++)
3367  {
3368  if (Z[j] * Z[index] > 0 && j != index)
3369  {
3370  sum += x[j] * Yniii[k][index][i][j];
3371  }
3372  }
3373 
3374  for (j = 1; j < aq.species_num; j ++)
3375  {
3376  for (l = 1; l < aq.species_num; l ++)
3377  {
3378  if ((Z[j] * Z[l] > 0) && (Z[j] * Z[i] < 0))
3379  {
3380  sum = sum - 0.5 * 3.0 * x[j] * x[l] * Yniii[k][i][j][l];
3381  }
3382  }
3383  }
3384 
3385  if (Z[i] * Z[index] < 0)
3386  {
3387  sum1 = sum1 + E[i] * sum;
3388  }
3389  }
3390 
3391  if (Z[k] == 0)
3392  {
3393  lnfi = lnfi + 4.0 * x[k] * x[k] * sum1;
3394  }
3395  }
3396 
3397  for (k = 0; k < aq.species_num; k ++)
3398  {
3399  double sum1 = 0;
3400 
3401  for (i = 0; i < aq.species_num; i ++)
3402  {
3403  double sum = 0;
3404 
3405  for (j = 0; j < aq.species_num; j ++)
3406  {
3407  for (l = 0; l < aq.species_num; l ++)
3408  {
3409  if ((Z[j] * Z[l] > 0) && (Z[j] * Z[i] < 0))
3410  {
3411  sum = sum + 0.5 * x[j] * x[l] * Yniii[k][i][j][l];
3412  }
3413  }
3414  }
3415 
3416  if (Z[i] * Z[index] > 0)
3417  {
3418  sum1 = sum1 + (3.0 * E[i] - iEj[index][i]) * sum;
3419  }
3420  }
3421 
3422  if (Z[k] == 0)
3423  {
3424  lnfi = lnfi - 4.0 * x[k] * x[k] * sum1;
3425  }
3426  }
3428 
3429  // neutral-neutral-cation-anion quaterary interaction
3430  for (i = 0; i < aq.species_num; i ++)
3431  {
3432  for (j = 0; j < aq.species_num; j ++)
3433  {
3434  double sum1 = 0;
3435  for (k = 0; k < aq.species_num; k ++)
3436  {
3437  double sum = 0;
3438  for (l = 0; l < aq.species_num; l ++)
3439  {
3440  if (Z[k] * Z[l] < 0 && Z[i] == 0 && Z[j] == 0)
3441  {
3442  sum += E[k] * (fabs (Z[l]) + fabs (Z[k])) / (fabs (Z[l]) * fabs (Z[k])) * Ynnii[i][j][k][l];
3443  }
3444  }
3445 
3446  if (Z[k] * Z[index] < 0)
3447  {
3448  sum1 += E[k] * ((fabs (Z[index]) + fabs (Z[k])) / fabs (Z[k]) * Ynnii[index][i][j][k] - (fabs (Z[index])/2 + 2/F) * sum);
3449  }
3450  }
3451 
3452  if (Z[i] == 0 && Z[j] == 0 && i != j)
3453  {
3454  lnfi = lnfi + 0.5 * x[i] * x[j] * sum1;
3455  }
3456  }
3457  }
3458 
3459  // neutral-neutral binary interaction
3460  for (i = 0; i < aq.species_num; i ++)
3461  {
3462  for (j = 0; j < i; j ++)
3463  {
3464  if (Z[i] == 0 && Z[j] == 0 && i != j)
3465  {
3466  lnfi = lnfi - 2.0 * x[i] * x[j] * (Wnn[i][j] + 2.0 * (x[i] - x[j]) * Unn[i][j]);
3467  }
3468  }
3469  }
3470 
3471  for (k = 0; k < aq.species_num; k ++)
3472  {
3473  if (Z[k] == 0 &&
3474  strcmp (aq.aqueous_species[k].species_symbol, WATER_SYMBOL) == 0)
3475  {
3476  break;
3477  }
3478  }
3479 
3480  return exp(lnfi) * x[k] * (n[index] / n[k] * 1000/CONSTANT_MW);
3481 }
3482 
3483 extern double psc_r (AQUEOUS_PHASE aq, int index, double *n, double T, double P)
3484 {
3485  if (aq.aqueous_species[index].charge == 0)
3486  {
3487  return psc_rn (aq, index, n, T, P);
3488  }
3489  else
3490  {
3491  return psc_ri (aq, index, n, T, P);
3492  }
3493 }
3494 
3495 
3496 extern double psc_a (AQUEOUS_PHASE aq, int index, double *n, double T, double P)
3497 {
3498  if (aq.aqueous_species[index].charge == 0)
3499  {
3500  return psc_an (aq, index, n, T, P);
3501  }
3502  else
3503  {
3504  return psc_ai (aq, index, n, T, P);
3505  }
3506 }
3507 
NNBINPARAM::species_s1
char species_s1[64]
Definition: islec.h:97
AQUEOUS_PHASE::ii_bin_param_num
int ii_bin_param_num
Definition: islec.h:158
NIITERPARAM::species_s3
char species_s3[64]
Definition: islec.h:107
IIBINPARAM::species_s1
char species_s1[64]
Definition: islec.h:87
p_T
#define p_T
Definition: psc_model.c:48
IIITERPARAM::species_s2
char species_s2[64]
Definition: islec.h:116
C4
#define C4
Definition: islec.h:45
IIBINPARAM::B0ss
double B0ss[6]
Definition: islec.h:91
IIITERPARAM::species_s3
char species_s3[64]
Definition: islec.h:117
AQUEOUS_PHASE::ii_bin_params
IIBINPARAM * ii_bin_params
Definition: islec.h:159
NIITERPARAM::species_s1
char species_s1[64]
Definition: islec.h:105
NIIIQUAPARAM::Qssss
double Qssss[6]
Definition: islec.h:127
AQUEOUS_PHASE::species_num
int species_num
Definition: islec.h:156
CONSTANT_k
#define CONSTANT_k
Definition: psc_model.c:35
psc_a
double psc_a(AQUEOUS_PHASE aq, int index, double *n, double T, double P)
Definition: psc_model.c:3496
AQUEOUS_PHASE::aqueous_species
SPECIES * aqueous_species
Definition: islec.h:157
AQUEOUS_PHASE
Definition: islec.h:154
NIITERPARAM::Usss
double Usss[6]
Definition: islec.h:109
AQUEOUS_PHASE::niii_qua_param_num
int niii_qua_param_num
Definition: islec.h:166
IIBINPARAM::alpha
double alpha
Definition: islec.h:89
NNIIQUAPARAM::species_s1
char species_s1[64]
Definition: islec.h:133
CONSTANT_pi
#define CONSTANT_pi
Definition: psc_model.c:33
ROU
#define ROU
Definition: psc_model.c:98
AQUEOUS_PHASE::nn_bin_params
NNBINPARAM * nn_bin_params
Definition: islec.h:161
C1
#define C1
Definition: islec.h:42
SPECIES::charge
int charge
Definition: islec.h:80
WATER_SYMBOL
#define WATER_SYMBOL
Definition: islec.h:47
AQUEOUS_PHASE::nii_ter_params
NIITERPARAM * nii_ter_params
Definition: islec.h:163
IIITERPARAM::Wsss
double Wsss[6]
Definition: islec.h:118
NNIIQUAPARAM::species_s2
char species_s2[64]
Definition: islec.h:134
AQUEOUS_PHASE::nn_bin_param_num
int nn_bin_param_num
Definition: islec.h:160
NIIIQUAPARAM::species_s2
char species_s2[64]
Definition: islec.h:124
dlnD_T
#define dlnD_T
Definition: psc_model.c:81
C2
#define C2
Definition: islec.h:43
NNBINPARAM::species_s2
char species_s2[64]
Definition: islec.h:98
CONSTANT_NA
#define CONSTANT_NA
Definition: psc_model.c:34
NNIIQUAPARAM::species_s3
char species_s3[64]
Definition: islec.h:135
NIITERPARAM::Wsss
double Wsss[6]
Definition: islec.h:108
NIITERPARAM::Vsss
double Vsss[6]
Definition: islec.h:110
NNIIQUAPARAM::species_s4
char species_s4[64]
Definition: islec.h:136
NIIIQUAPARAM::species_s4
char species_s4[64]
Definition: islec.h:126
NIIIQUAPARAM::species_s3
char species_s3[64]
Definition: islec.h:125
R
#define R
Definition: islec.h:37
C3
#define C3
Definition: islec.h:44
dp_T
#define dp_T
Definition: psc_model.c:67
IIBINPARAM::B1ss
double B1ss[6]
Definition: islec.h:92
NIIIQUAPARAM::Yssss
double Yssss[6]
Definition: islec.h:128
AQUEOUS_PHASE::nii_ter_param_num
int nii_ter_param_num
Definition: islec.h:162
NNIIQUAPARAM::Yssss
double Yssss[6]
Definition: islec.h:137
NIITERPARAM::species_s2
char species_s2[64]
Definition: islec.h:106
d2D_T
#define d2D_T
Definition: psc_model.c:88
AQUEOUS_PHASE::iii_ter_params
IIITERPARAM * iii_ter_params
Definition: islec.h:165
d2p_T
#define d2p_T
Definition: psc_model.c:91
AQUEOUS_PHASE::nnii_qua_params
NNIIQUAPARAM * nnii_qua_params
Definition: islec.h:169
SPECIES::species_symbol
char species_symbol[64]
Definition: islec.h:79
CONSTANT_MW
#define CONSTANT_MW
Definition: psc_model.c:32
IIBINPARAM::species_s2
char species_s2[64]
Definition: islec.h:88
IIBINPARAM::alpha1
double alpha1
Definition: islec.h:90
psc_r
double psc_r(AQUEOUS_PHASE aq, int index, double *n, double T, double P)
Definition: psc_model.c:3483
NIIIQUAPARAM::species_s1
char species_s1[64]
Definition: islec.h:123
AQUEOUS_PHASE::niii_qua_params
NIIIQUAPARAM * niii_qua_params
Definition: islec.h:167
aw_T
#define aw_T
Definition: psc_model.c:74
AQUEOUS_PHASE::iii_ter_param_num
int iii_ter_param_num
Definition: islec.h:164
IIITERPARAM::species_s1
char species_s1[64]
Definition: islec.h:115
CONSTANT_e
#define CONSTANT_e
Definition: psc_model.c:36
AQUEOUS_PHASE::nnii_qua_param_num
int nnii_qua_param_num
Definition: islec.h:168
NNBINPARAM::Uss
double Uss[6]
Definition: islec.h:100
D_T
#define D_T
Definition: psc_model.c:65
NNBINPARAM::Wss
double Wss[6]
Definition: islec.h:99
Generated by   doxygen 1.8.13