40 m_formPitzerTemp(PITZER_TEMP_CONSTANT),
42 m_IionicMolality(0.0),
43 m_maxIionicStrength(100.0),
44 m_TempPitzerRef(298.15),
45 m_IionicMolalityStoich(0.0),
46 m_form_A_Debye(A_DEBYE_WATER),
51 m_molalitiesAreCropped(false),
54 IMS_gamma_o_min_(1.0E-5),
55 IMS_gamma_k_min_(10.0),
75 CROP_ln_gamma_o_min(-6.0),
76 CROP_ln_gamma_o_max(3.0),
77 CROP_ln_gamma_k_min(-5.0),
78 CROP_ln_gamma_k_max(15.0),
81 for (
size_t i = 0; i < 17; i++) {
97 m_formPitzerTemp(PITZER_TEMP_CONSTANT),
99 m_IionicMolality(0.0),
100 m_maxIionicStrength(100.0),
101 m_TempPitzerRef(298.15),
102 m_IionicMolalityStoich(0.0),
103 m_form_A_Debye(A_DEBYE_WATER),
106 m_densWaterSS(1000.),
108 m_molalitiesAreCropped(false),
110 IMS_X_o_cutoff_(0.2),
111 IMS_gamma_o_min_(1.0E-5),
112 IMS_gamma_k_min_(10.0),
132 CROP_ln_gamma_o_min(-6.0),
133 CROP_ln_gamma_o_max(3.0),
134 CROP_ln_gamma_k_min(-5.0),
135 CROP_ln_gamma_k_max(15.0),
138 for (
int i = 0; i < 17; i++) {
148 m_formPitzerTemp(PITZER_TEMP_CONSTANT),
150 m_IionicMolality(0.0),
151 m_maxIionicStrength(100.0),
152 m_TempPitzerRef(298.15),
153 m_IionicMolalityStoich(0.0),
154 m_form_A_Debye(A_DEBYE_WATER),
157 m_densWaterSS(1000.),
159 m_molalitiesAreCropped(false),
161 IMS_X_o_cutoff_(0.2),
162 IMS_gamma_o_min_(1.0E-5),
163 IMS_gamma_k_min_(10.0),
183 CROP_ln_gamma_o_min(-6.0),
184 CROP_ln_gamma_o_max(3.0),
185 CROP_ln_gamma_k_min(-5.0),
186 CROP_ln_gamma_k_max(15.0),
189 for (
int i = 0; i < 17; i++) {
205 m_formPitzerTemp(PITZER_TEMP_CONSTANT),
207 m_IionicMolality(0.0),
208 m_maxIionicStrength(100.0),
209 m_TempPitzerRef(298.15),
210 m_IionicMolalityStoich(0.0),
211 m_form_A_Debye(A_DEBYE_WATER),
214 m_densWaterSS(1000.),
216 m_molalitiesAreCropped(false),
218 IMS_X_o_cutoff_(0.2),
219 IMS_gamma_o_min_(1.0E-5),
220 IMS_gamma_k_min_(10.0),
240 CROP_ln_gamma_o_min(-6.0),
241 CROP_ln_gamma_o_max(3.0),
242 CROP_ln_gamma_k_min(-5.0),
243 CROP_ln_gamma_k_max(15.0),
278 throw CanteraError(
"HMWSoln::operator=()",
"Dynamic cast to PDSS_Water failed");
440 m_formPitzerTemp(PITZER_TEMP_CONSTANT),
442 m_IionicMolality(0.0),
443 m_maxIionicStrength(100.0),
444 m_TempPitzerRef(298.15),
445 m_IionicMolalityStoich(0.0),
446 m_form_A_Debye(A_DEBYE_WATER),
449 m_densWaterSS(1000.),
451 m_molalitiesAreCropped(false),
453 IMS_X_o_cutoff_(0.2),
454 IMS_gamma_o_min_(1.0E-5),
455 IMS_gamma_k_min_(10.0),
475 CROP_ln_gamma_o_min(-6.0),
476 CROP_ln_gamma_o_max(3.0),
477 CROP_ln_gamma_k_min(-5.0),
478 CROP_ln_gamma_k_max(15.0),
482 printf(
"unknown test problem\n");
490 size_t n = i *
m_kk + j;
532 double param = -0.004;
649 for (
size_t k = 0; k <
m_kk; k++) {
653 return (hbar - h0bar);
663 size_t kcation =
npos;
664 double xcation = 0.0;
665 size_t kanion =
npos;
667 for (
size_t k = 0; k <
m_kk; k++) {
668 if (charge[k] > 0.0) {
673 }
else if (charge[k] < 0.0) {
674 if (
m_tmpV[k] > xcation) {
680 if (kcation ==
npos || kanion ==
npos) {
683 double xuse = xcation;
685 if (xanion < xcation) {
687 if (charge[kcation] != 1.0) {
688 factor = charge[kcation];
691 if (charge[kanion] != 1.0) {
692 factor = charge[kanion];
695 xuse = xuse / factor;
711 double uu = hh - pres * molarV;
753 double cv = cp - beta * beta * tt * molarV / kappa_t;
783 double* vbar = &
m_pp[0];
787 doublereal vtotal = 0.0;
788 for (
size_t i = 0; i <
m_kk; i++) {
789 vtotal += vbar[i] * x[i];
807 throw CanteraError(
"HMWSoln::isothermalCompressibility",
860 if (rho != dens_old) {
862 "Density is not an independent variable");
878 "Density is not an independent variable");
946 for (
size_t k = 1; k <
m_kk; k++) {
981 return 1.0 / mvSolvent;
991 return log(c_solvent);
1018 for (
int i = 0; i < sizeUA; i++) {
1023 uA[1] = -int(
nDim());
1058 for (
size_t k = 0; k <
m_kk; k++) {
1090 for (
size_t k = 0; k <
m_kk; k++) {
1091 acMolality[k] = exp(acMolality[k]);
1116 const double xxSmall = 1.0E-150;
1131 for (
size_t k = 0; k <
m_kk; k++) {
1137 xx =
std::max(xmolSolvent, xxSmall);
1173 for (
size_t k = 0; k <
m_kk; k++) {
1182 double RTT = RT * T;
1183 for (
size_t k = 0; k <
m_kk; k++) {
1229 for (
size_t k = 0; k <
m_kk; k++) {
1242 for (
size_t k = 0; k <
m_kk; k++) {
1258 for (
size_t k = 0; k <
m_kk; k++) {
1294 for (
size_t k = 0; k <
m_kk; k++) {
1325 for (
size_t k = 0; k <
m_kk; k++) {
1337 double RTT = RT * T;
1338 for (
size_t k = 0; k <
m_kk; k++) {
1467 if (tempArg != -1.0) {
1471 if (presArg != -1.0) {
1484 printf(
"shouldn't be here\n");
1503 if (tempArg != -1.0) {
1507 if (presArg != -1.0) {
1520 printf(
"shouldn't be here\n");
1539 if (tempArg != -1.0) {
1543 if (presArg != -1.0) {
1555 printf(
"shouldn't be here\n");
1575 double dAphidT = dAdT /3.0;
1577 if (tempArg != -1.0) {
1580 double retn = dAphidT * (4.0 *
GasConstant * T * T);
1597 double dAphidP = dAdP /3.0;
1599 if (tempArg != -1.0) {
1602 double retn = - dAphidP * (4.0 *
GasConstant * T);
1627 if (tempArg != -1.0) {
1632 double d2Aphi = d2 / 3.0;
1633 double retn = 2.0 * A_L / T + 4.0 *
GasConstant * T * T *d2Aphi;
1650 if (tempArg != -1.0) {
1654 if (presArg != -1.0) {
1666 printf(
"shouldn't be here\n");
1695 "Unfinished func called: " + msg);
1727 size_t maxCounterIJlen = 1 + (
m_kk-1) * (
m_kk-2) / 2;
1733 int TCoeffLength = 1;
1807 m_BMX_IJ.resize(maxCounterIJlen, 0.0);
1819 m_Phi_IJ.resize(maxCounterIJlen, 0.0);
1828 m_CMX_IJ.resize(maxCounterIJlen, 0.0);
1865 for (
size_t k = 0; k <
m_kk; k++) {
1871 double zs_k2 = z_k - zs_k1;
1896 double lnActCoeffMolal0 = - log(xx) + (xx - 1.0)/xx;
1897 double lnxs = log(xx);
1899 for (
size_t k = 1; k <
m_kk; k++) {
1940 doublereal Imax = 0.0, Itmp;
1944 for (
size_t k = 0; k <
m_kk; k++) {
1956 if (cropMethod == 0) {
1967 for (
size_t i = 1; i < (m_kk - 1); i++) {
1969 double abs_charge_i = fabs(charge_i);
1970 if (charge_i == 0.0) {
1973 for (
size_t j = (i+1); j <
m_kk; j++) {
1975 double abs_charge_j = fabs(charge_j);
1984 if (charge_i * charge_j < 0) {
1989 if (Imax > Iac_max) {
1993 if (Imax > Iac_max) {
1998 if (Imax > Iac_max) {
2002 if (Imax > Iac_max) {
2014 for (
int times = 0; times< 10; times++) {
2015 double anion_charge = 0.0;
2016 double cation_charge = 0.0;
2017 size_t anion_contrib_max_i =
npos;
2018 double anion_contrib_max = -1.0;
2019 size_t cation_contrib_max_i =
npos;
2020 double cation_contrib_max = -1.0;
2021 for (
size_t i = 0; i <
m_kk; i++) {
2023 if (charge_i < 0.0) {
2025 anion_charge += anion_contrib ;
2026 if (anion_contrib > anion_contrib_max) {
2027 anion_contrib_max = anion_contrib;
2028 anion_contrib_max_i = i;
2030 }
else if (charge_i > 0.0) {
2032 cation_charge += cation_contrib ;
2033 if (cation_contrib > cation_contrib_max) {
2034 cation_contrib_max = cation_contrib;
2035 cation_contrib_max_i = i;
2039 double total_charge = cation_charge - anion_charge;
2040 if (total_charge > 1.0E-8) {
2041 double desiredCrop = total_charge/
m_speciesCharge[cation_contrib_max_i];
2043 if (desiredCrop < maxCrop) {
2049 }
else if (total_charge < -1.0E-8) {
2050 double desiredCrop = total_charge/
m_speciesCharge[anion_contrib_max_i];
2052 if (desiredCrop < maxCrop) {
2064 if (cropMethod == 1) {
2078 for (
size_t k = 0; k <
m_kk; k++) {
2086 for (
size_t k = 0; k <
m_kk; k++) {
2092 for (
size_t k = 0; k <
m_kk; k++) {
2094 if (charge != 0.0) {
2115 for (i = 0; i <
m_kk; i++) {
2121 for (i = 1; i < (m_kk - 1); i++) {
2124 for (j = (i+1); j <
m_kk; j++) {
2150 size_t i, j, n, counterIJ;
2151 const double* beta0MX_coeff;
2152 const double* beta1MX_coeff;
2153 const double* beta2MX_coeff;
2154 const double* CphiMX_coeff;
2155 const double* Theta_coeff;
2158 double tinv = 0.0, tln = 0.0, tlin = 0.0, tquad = 0.0;
2163 tquad = T * T - Tr * Tr;
2165 tinv = 1.0/T - 1.0/Tr;
2168 for (i = 1; i < (
m_kk - 1); i++) {
2169 for (j = (i+1); j <
m_kk; j++) {
2184 case PITZER_TEMP_CONSTANT:
2186 case PITZER_TEMP_LINEAR:
2189 + beta0MX_coeff[1]*tlin;
2194 + beta1MX_coeff[1]*tlin;
2199 + beta2MX_coeff[1]*tlin;
2204 + CphiMX_coeff[1]*tlin;
2208 m_Theta_ij[counterIJ] = Theta_coeff[0] + Theta_coeff[1]*tlin;
2214 case PITZER_TEMP_COMPLEX1:
2216 + beta0MX_coeff[1]*tlin
2217 + beta0MX_coeff[2]*tquad
2218 + beta0MX_coeff[3]*tinv
2219 + beta0MX_coeff[4]*tln;
2222 + beta1MX_coeff[1]*tlin
2223 + beta1MX_coeff[2]*tquad
2224 + beta1MX_coeff[3]*tinv
2225 + beta1MX_coeff[4]*tln;
2228 + beta2MX_coeff[1]*tlin
2229 + beta2MX_coeff[2]*tquad
2230 + beta2MX_coeff[3]*tinv
2231 + beta2MX_coeff[4]*tln;
2234 + CphiMX_coeff[1]*tlin
2235 + CphiMX_coeff[2]*tquad
2236 + CphiMX_coeff[3]*tinv
2237 + CphiMX_coeff[4]*tln;
2240 + Theta_coeff[1]*tlin
2241 + Theta_coeff[2]*tquad
2242 + Theta_coeff[3]*tinv
2243 + Theta_coeff[4]*tln;
2246 + beta0MX_coeff[2]*2.0*T
2247 - beta0MX_coeff[3]/(T*T)
2248 + beta0MX_coeff[4]/T;
2251 + beta1MX_coeff[2]*2.0*T
2252 - beta1MX_coeff[3]/(T*T)
2253 + beta1MX_coeff[4]/T;
2256 + beta2MX_coeff[2]*2.0*T
2257 - beta2MX_coeff[3]/(T*T)
2258 + beta2MX_coeff[4]/T;
2261 + CphiMX_coeff[2]*2.0*T
2262 - CphiMX_coeff[3]/(T*T)
2263 + CphiMX_coeff[4]/T;
2266 + Theta_coeff[2]*2.0*T
2267 - Theta_coeff[3]/(T*T)
2273 + beta0MX_coeff[2]*2.0
2274 + 2.0*beta0MX_coeff[3]/(T*T*T)
2275 - beta0MX_coeff[4]/(T*T);
2278 + beta1MX_coeff[2]*2.0
2279 + 2.0*beta1MX_coeff[3]/(T*T*T)
2280 - beta1MX_coeff[4]/(T*T);
2283 + beta2MX_coeff[2]*2.0
2284 + 2.0*beta2MX_coeff[3]/(T*T*T)
2285 - beta2MX_coeff[4]/(T*T);
2288 + CphiMX_coeff[2]*2.0
2289 + 2.0*CphiMX_coeff[3]/(T*T*T)
2290 - CphiMX_coeff[4]/(T*T);
2293 + Theta_coeff[2]*2.0
2294 + 2.0*Theta_coeff[3]/(T*T*T)
2295 - Theta_coeff[4]/(T*T);
2317 for (i = 1; i <
m_kk; i++) {
2319 for (j = 1; j <
m_kk; j++) {
2323 case PITZER_TEMP_CONSTANT:
2326 case PITZER_TEMP_LINEAR:
2327 m_Lambda_nj(i,j) = Lambda_coeff[0] + Lambda_coeff[1]*tlin;
2331 case PITZER_TEMP_COMPLEX1:
2333 + Lambda_coeff[1]*tlin
2334 + Lambda_coeff[2]*tquad
2335 + Lambda_coeff[3]*tinv
2336 + Lambda_coeff[4]*tln;
2339 + Lambda_coeff[2]*2.0*T
2340 - Lambda_coeff[3]/(T*T)
2341 + Lambda_coeff[4]/T;
2345 + 2.0*Lambda_coeff[3]/(T*T*T)
2346 - Lambda_coeff[4]/(T*T);
2352 case PITZER_TEMP_CONSTANT:
2355 case PITZER_TEMP_LINEAR:
2356 m_Mu_nnn[i] = Mu_coeff[0] + Mu_coeff[1]*tlin;
2360 case PITZER_TEMP_COMPLEX1:
2374 + 2.0*Mu_coeff[3]/(T*T*T)
2375 - Mu_coeff[4]/(T*T);
2383 for (i = 1; i <
m_kk; i++) {
2384 for (j = 1; j <
m_kk; j++) {
2385 for (
size_t k = 1; k <
m_kk; k++) {
2386 n = i * m_kk *m_kk + j * m_kk + k ;
2389 case PITZER_TEMP_CONSTANT:
2392 case PITZER_TEMP_LINEAR:
2393 m_Psi_ijk[n] = Psi_coeff[0] + Psi_coeff[1]*tlin;
2397 case PITZER_TEMP_COMPLEX1:
2400 + Psi_coeff[2]*tquad
2405 + Psi_coeff[2]*2.0*T
2406 - Psi_coeff[3]/(T*T)
2411 + 2.0*Psi_coeff[3]/(T*T*T)
2412 - Psi_coeff[4]/(T*T);
2436 printf(
"Wrong index solvent value!\n");
2446 std::string sni, snj, snk;
2477 double etheta[5][5], etheta_prime[5][5], sqrtIs;
2486 double molarcharge = 0.0;
2491 double molalitysumUncropped = 0.0;
2507 double Aphi, F, zsqF;
2508 double sum1, sum2, sum3, sum4, sum5, term1;
2509 double sum_m_phi_minus_1, osmotic_coef, lnwateract;
2512 size_t n, i, j, k, m, counterIJ, counterIJ2;
2516 printf(
"\n Debugging information from hmw_act \n");
2527 for (n = 1; n <
m_kk; n++) {
2529 Is += charge[n] * charge[n] * molality[n];
2531 molarcharge += fabs(charge[n]) * molality[n];
2543 printf(
" Step 1: \n");
2544 printf(
" ionic strenth = %14.7le \n total molar "
2545 "charge = %14.7le \n", Is, molarcharge);
2563 printf(
" Step 2: \n");
2566 for (z1 = 1; z1 <=4; z1++) {
2567 for (z2 =1; z2 <=4; z2++) {
2568 calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]);
2571 printf(
" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n",
2572 z1, z2, etheta[z1][z2], etheta_prime[z1][z2]);
2580 printf(
" Step 3: \n");
2581 printf(
" Species Species g(x) "
2592 for (i = 1; i < (m_kk - 1); i++) {
2593 for (j = (i+1); j <
m_kk; j++) {
2603 if (charge[i]*charge[j] < 0) {
2607 x1 = sqrtIs * alpha1MX[counterIJ];
2608 if (x1 > 1.0E-100) {
2609 gfunc[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 * x1);
2610 hfunc[counterIJ] = -2.0 *
2611 (1.0-(1.0 + x1 + 0.5 * x1 * x1) * exp(-x1)) / (x1 * x1);
2613 gfunc[counterIJ] = 0.0;
2614 hfunc[counterIJ] = 0.0;
2617 if (beta2MX[counterIJ] != 0.0) {
2618 x2 = sqrtIs * alpha2MX[counterIJ];
2619 if (x2 > 1.0E-100) {
2620 g2func[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
2621 h2func[counterIJ] = -2.0 *
2622 (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
2624 g2func[counterIJ] = 0.0;
2625 h2func[counterIJ] = 0.0;
2629 gfunc[counterIJ] = 0.0;
2630 hfunc[counterIJ] = 0.0;
2636 printf(
" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(),
2637 gfunc[counterIJ], hfunc[counterIJ]);
2649 printf(
" Step 4: \n");
2650 printf(
" Species Species BMX "
2651 "BprimeMX BphiMX \n");
2655 for (i = 1; i < m_kk - 1; i++) {
2656 for (j = i+1; j <
m_kk; j++) {
2665 if (counterIJ == 2) {
2668 printf(
"beta0MX[%d] = %g\n", counterIJ, beta0MX[counterIJ]);
2669 printf(
"beta1MX[%d] = %g\n", counterIJ, beta1MX[counterIJ]);
2670 printf(
"beta2MX[%d] = %g\n", counterIJ, beta2MX[counterIJ]);
2678 if (charge[i]*charge[j] < 0.0) {
2679 BMX[counterIJ] = beta0MX[counterIJ]
2680 + beta1MX[counterIJ] * gfunc[counterIJ]
2681 + beta2MX[counterIJ] * g2func[counterIJ];
2684 printf(
"%d %g: %g %g %g %g\n",
2685 counterIJ, BMX[counterIJ], beta0MX[counterIJ],
2686 beta1MX[counterIJ], beta2MX[counterIJ], gfunc[counterIJ]);
2689 if (Is > 1.0E-150) {
2690 BprimeMX[counterIJ] = (beta1MX[counterIJ] * hfunc[counterIJ]/Is +
2691 beta2MX[counterIJ] * h2func[counterIJ]/Is);
2693 BprimeMX[counterIJ] = 0.0;
2695 BphiMX[counterIJ] = BMX[counterIJ] + Is*BprimeMX[counterIJ];
2697 BMX[counterIJ] = 0.0;
2698 BprimeMX[counterIJ] = 0.0;
2699 BphiMX[counterIJ] = 0.0;
2705 printf(
" %-16s %-16s %11.7f %11.7f %11.7f \n",
2706 sni.c_str(), snj.c_str(),
2707 BMX[counterIJ], BprimeMX[counterIJ], BphiMX[counterIJ]);
2719 printf(
" Step 5: \n");
2720 printf(
" Species Species CMX \n");
2723 for (i = 1; i < m_kk-1; i++) {
2724 for (j = i+1; j <
m_kk; j++) {
2734 if (charge[i]*charge[j] < 0.0) {
2735 CMX[counterIJ] = CphiMX[counterIJ]/
2736 (2.0* sqrt(fabs(charge[i]*charge[j])));
2738 CMX[counterIJ] = 0.0;
2742 if (counterIJ == 2) {
2745 printf(
"CphiMX[%d] = %g\n", counterIJ, CphiMX[counterIJ]);
2753 printf(
" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(),
2766 printf(
" Step 6: \n");
2767 printf(
" Species Species Phi_ij "
2768 " Phiprime_ij Phi^phi_ij \n");
2771 for (i = 1; i < m_kk-1; i++) {
2772 for (j = i+1; j <
m_kk; j++) {
2782 if (charge[i]*charge[j] > 0) {
2783 z1 = (int) fabs(charge[i]);
2784 z2 = (int) fabs(charge[j]);
2785 Phi[counterIJ] = thetaij[counterIJ] + etheta[z1][z2];
2786 Phiprime[counterIJ] = etheta_prime[z1][z2];
2787 Phiphi[counterIJ] = Phi[counterIJ] + Is * Phiprime[counterIJ];
2789 Phi[counterIJ] = 0.0;
2790 Phiprime[counterIJ] = 0.0;
2791 Phiphi[counterIJ] = 0.0;
2797 printf(
" %-16s %-16s %10.6f %10.6f %10.6f \n",
2798 sni.c_str(), snj.c_str(),
2799 Phi[counterIJ], Phiprime[counterIJ], Phiphi[counterIJ]);
2811 printf(
" Step 7: \n");
2819 F = -Aphi * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
2820 + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
2823 printf(
"Aphi = %20.13g\n", Aphi);
2828 printf(
" initial value of F = %10.6f \n", F);
2831 for (i = 1; i < m_kk-1; i++) {
2832 for (j = i+1; j <
m_kk; j++) {
2842 if (charge[i]*charge[j] < 0) {
2843 F = F + molality[i]*molality[j] * BprimeMX[counterIJ];
2849 if (charge[i]*charge[j] > 0) {
2850 F = F + molality[i]*molality[j] * Phiprime[counterIJ];
2854 printf(
" F = %10.6f \n", F);
2861 printf(
" Step 8: Summing in All Contributions to Activity Coefficients \n");
2865 for (i = 1; i <
m_kk; i++) {
2872 if (charge[i] > 0.0) {
2877 printf(
" Contributions to ln(ActCoeff_%s):\n", sni.c_str());
2881 zsqF = charge[i]*charge[i]*F;
2884 printf(
" Unary term: z*z*F = %10.5f\n", zsqF);
2892 for (j = 1; j <
m_kk; j++) {
2899 if (charge[j] < 0.0) {
2901 sum1 = sum1 + molality[j]*
2902 (2.0*BMX[counterIJ] + molarcharge*CMX[counterIJ]);
2906 printf(
" Bin term with %-13s 2 m_j BMX = %10.5f\n", snj.c_str(),
2907 molality[j]*2.0*BMX[counterIJ]);
2908 printf(
" m_j Z CMX = %10.5f\n",
2909 molality[j]* molarcharge*CMX[counterIJ]);
2918 for (k = j+1; k <
m_kk; k++) {
2920 if (charge[k] < 0.0) {
2921 n = k + j * m_kk + i * m_kk *
m_kk;
2922 sum3 = sum3 + molality[j]*molality[k]*psi_ijk[n];
2925 if (psi_ijk[n] != 0.0) {
2927 printf(
" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(),
2928 molality[j]*molality[k]*psi_ijk[n]);
2938 if (charge[j] > 0.0) {
2941 sum2 = sum2 + molality[j]*(2.0*Phi[counterIJ]);
2944 if ((molality[j] * Phi[counterIJ])!= 0.0) {
2946 printf(
" Phi term with %-12s 2 m_j Phi_cc = %10.5f\n", snj.c_str(),
2947 molality[j]*(2.0*Phi[counterIJ]));
2952 for (k = 1; k <
m_kk; k++) {
2953 if (charge[k] < 0.0) {
2956 n = k + j * m_kk + i * m_kk *
m_kk;
2957 sum2 = sum2 + molality[j]*molality[k]*psi_ijk[n];
2960 if (psi_ijk[n] != 0.0) {
2962 printf(
" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(),
2963 molality[j]*molality[k]*psi_ijk[n]);
2972 sum4 = sum4 + (fabs(charge[i])*
2973 molality[j]*molality[k]*CMX[counterIJ2]);
2976 if ((molality[j]*molality[k]*CMX[counterIJ2]) != 0.0) {
2978 printf(
" Tern CMX term on %-16s abs(z_i) m_j m_k CMX = %10.5f\n", snj.c_str(),
2979 fabs(charge[i])* molality[j]*molality[k]*CMX[counterIJ2]);
2990 if (charge[j] == 0) {
2996 printf(
" Lambda term with %-12s 2 m_j lam_ji = %10.5f\n", snj.c_str(),
3004 for (
size_t k = 1; k <
m_kk; k++) {
3005 if (charge[k] < 0.0) {
3008 n = izeta * m_kk * m_kk + jzeta * m_kk + k;
3009 double zeta = psi_ijk[n];
3011 sum5 = sum5 + molality[j]*molality[k]*zeta;
3015 printf(
" Zeta term on %-16s m_n m_a zeta_nMa = %10.5f\n", snj.c_str(),
3016 molality[j]*molality[k]*psi_ijk[n]);
3033 printf(
" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f \n",
3044 if (charge[i] < 0) {
3049 printf(
" Contributions to ln(ActCoeff_%s):\n", sni.c_str());
3054 zsqF = charge[i]*charge[i]*F;
3057 printf(
" Unary term: z*z*F = %10.5f\n", zsqF);
3065 for (j = 1; j <
m_kk; j++) {
3075 if (charge[j] > 0) {
3076 sum1 = sum1 + molality[j]*
3077 (2.0*BMX[counterIJ]+molarcharge*CMX[counterIJ]);
3081 printf(
" Bin term with %-13s 2 m_j BMX = %10.5f\n", snj.c_str(),
3082 molality[j]*2.0*BMX[counterIJ]);
3083 printf(
" m_j Z CMX = %10.5f\n",
3084 molality[j]* molarcharge*CMX[counterIJ]);
3088 for (k = j+1; k <
m_kk; k++) {
3090 if (charge[k] > 0) {
3091 n = k + j * m_kk + i * m_kk *
m_kk;
3092 sum3 = sum3 + molality[j]*molality[k]*psi_ijk[n];
3095 if (psi_ijk[n] != 0.0) {
3097 printf(
" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(),
3098 molality[j]*molality[k]*psi_ijk[n]);
3110 if (charge[j] < 0.0) {
3113 sum2 = sum2 + molality[j]*(2.0*Phi[counterIJ]);
3116 if ((molality[j] * Phi[counterIJ])!= 0.0) {
3118 printf(
" Phi term with %-12s 2 m_j Phi_aa = %10.5f\n", snj.c_str(),
3119 molality[j]*(2.0*Phi[counterIJ]));
3124 for (k = 1; k <
m_kk; k++) {
3125 if (charge[k] > 0.0) {
3127 n = k + j * m_kk + i * m_kk *
m_kk;
3128 sum2 = sum2 + molality[j]*molality[k]*psi_ijk[n];
3131 if (psi_ijk[n] != 0.0) {
3133 printf(
" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(),
3134 molality[j]*molality[k]*psi_ijk[n]);
3145 molality[j]*molality[k]*CMX[counterIJ2]);
3148 if ((molality[j]*molality[k]*CMX[counterIJ2]) != 0.0) {
3150 printf(
" Tern CMX term on %-16s abs(z_i) m_j m_k CMX = %10.5f\n", snj.c_str(),
3151 fabs(charge[i])* molality[j]*molality[k]*CMX[counterIJ2]);
3162 if (charge[j] == 0.0) {
3168 printf(
" Lambda term with %-12s 2 m_j lam_ji = %10.5f\n", snj.c_str(),
3176 for (
size_t k = 1; k <
m_kk; k++) {
3177 if (charge[k] > 0.0) {
3181 n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
3182 double zeta = psi_ijk[n];
3184 sum5 = sum5 + molality[j]*molality[k]*zeta;
3188 printf(
" Zeta term on %-16s m_n m_c zeta_ncX = %10.5f\n", snj.c_str(),
3189 molality[j]*molality[k]*psi_ijk[n]);
3202 printf(
" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f\n",
3212 if (charge[i] == 0.0) {
3216 printf(
" Contributions to ln(ActCoeff_%s):\n", sni.c_str());
3221 for (j = 1; j <
m_kk; j++) {
3227 printf(
" Lambda_n term on %-16s 2 m_j lambda_n_j = %10.5f\n", snj.c_str(),
3236 if (charge[j] > 0.0) {
3237 for (k = 1; k <
m_kk; k++) {
3238 if (charge[k] < 0.0) {
3239 n = k + j * m_kk + i * m_kk *
m_kk;
3240 sum3 = sum3 + molality[j]*molality[k]*psi_ijk[n];
3243 if (psi_ijk[n] != 0.0) {
3245 printf(
" Zeta term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(),
3246 molality[j]*molality[k]*psi_ijk[n]);
3254 sum2 = 3.0 * molality[i]* molality[i] *
m_Mu_nnn[i];
3257 if (m_Mu_nnn[i] != 0.0) {
3258 printf(
" Mu_nnn term 3 m_n m_n Mu_n_n = %10.5f\n",
3259 3.0 * molality[i]* molality[i] * m_Mu_nnn[i]);
3269 printf(
" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f\n",
3278 printf(
" Step 9: \n");
3300 term1 = -Aphi * pow(Is,1.5) / (1.0 + 1.2 * sqrt(Is));
3302 for (j = 1; j <
m_kk; j++) {
3306 if (charge[j] > 0.0) {
3307 for (k = 1; k <
m_kk; k++) {
3308 if (charge[k] < 0.0) {
3315 sum1 = sum1 + molality[j]*molality[k]*
3316 (BphiMX[counterIJ] + molarcharge*CMX[counterIJ]);
3320 for (k = j+1; k <
m_kk; k++) {
3321 if (j == (m_kk-1)) {
3323 printf(
"logic error 1 in Step 9 of hmw_act");
3326 if (charge[k] > 0.0) {
3333 sum2 = sum2 + molality[j]*molality[k]*Phiphi[counterIJ];
3334 for (m = 1; m <
m_kk; m++) {
3335 if (charge[m] < 0.0) {
3337 n = m + k * m_kk + j * m_kk *
m_kk;
3339 molality[j]*molality[k]*molality[m]*psi_ijk[n];
3349 if (charge[j] < 0) {
3350 for (k = j+1; k <
m_kk; k++) {
3353 printf(
"logic error 2 in Step 9 of hmw_act");
3356 if (charge[k] < 0) {
3364 sum3 = sum3 + molality[j]*molality[k]*Phiphi[counterIJ];
3365 for (m = 1; m <
m_kk; m++) {
3366 if (charge[m] > 0.0) {
3367 n = m + k * m_kk + j * m_kk *
m_kk;
3369 molality[j]*molality[k]*molality[m]*psi_ijk[n];
3379 if (charge[j] == 0) {
3380 for (k = 1; k <
m_kk; k++) {
3381 if (charge[k] < 0.0) {
3382 sum4 = sum4 + molality[j]*molality[k]*
m_Lambda_nj(j,k);
3384 if (charge[k] > 0.0) {
3385 sum5 = sum5 + molality[j]*molality[k]*
m_Lambda_nj(j,k);
3387 if (charge[k] == 0.0) {
3389 sum6 = sum6 + molality[j]*molality[k]*
m_Lambda_nj(j,k);
3390 }
else if (k == j) {
3391 sum6 = sum6 + 0.5 * molality[j]*molality[k]*
m_Lambda_nj(j,k);
3394 if (charge[k] < 0.0) {
3396 for (m = 1; m <
m_kk; m++) {
3397 if (charge[m] > 0.0) {
3399 n = k + jzeta * m_kk + izeta * m_kk *
m_kk;
3400 double zeta = psi_ijk[n];
3402 sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta;
3408 sum7 += molality[j]*molality[j]*molality[j]*
m_Mu_nnn[j];
3411 sum_m_phi_minus_1 = 2.0 *
3412 (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
3417 if (molalitysumUncropped > 1.0E-150) {
3418 osmotic_coef = 1.0 + (sum_m_phi_minus_1 / molalitysumUncropped);
3424 printf(
"OsmCoef - 1 = %20.13g\n", osmotic_coef - 1.0);
3429 printf(
" term1=%10.6f sum1=%10.6f sum2=%10.6f "
3430 "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n",
3431 term1, sum1, sum2, sum3, sum4, sum5);
3432 printf(
" sum_m_phi_minus_1=%10.6f osmotic_coef=%10.6f\n",
3433 sum_m_phi_minus_1, osmotic_coef);
3437 printf(
" Step 10: \n");
3440 lnwateract = -(
m_weightSolvent/1000.0) * molalitysumUncropped * osmotic_coef;
3455 double wateract = exp(lnwateract);
3457 printf(
" molalitySumUncropped = %16.7g\n", molalitysumUncropped);
3458 printf(
" ln_a_water=%10.6f a_water=%10.6f\n\n",
3459 lnwateract, wateract);
3491 for (
size_t k = 1; k <
m_kk; k++) {
3533 printf(
"Wrong index solvent value!\n");
3537 std::string sni, snj, snk;
3553 double etheta[5][5], etheta_prime[5][5], sqrtIs;
3562 double molarcharge = 0.0;
3567 double molalitysum = 0.0;
3582 double dFdT, zsqdFdT;
3583 double sum1, sum2, sum3, sum4, sum5, term1;
3584 double sum_m_phi_minus_1, d_osmotic_coef_dT, d_lnwateract_dT;
3587 size_t n, i, j, k, m, counterIJ, counterIJ2;
3591 printf(
"\n Debugging information from "
3592 "s_Pitzer_dlnMolalityActCoeff_dT()\n");
3603 for (n = 1; n <
m_kk; n++) {
3605 Is += charge[n] * charge[n] * molality[n];
3607 molarcharge += fabs(charge[n]) * molality[n];
3608 molalitysum += molality[n];
3619 printf(
" Step 1: \n");
3620 printf(
" ionic strenth = %14.7le \n total molar "
3621 "charge = %14.7le \n", Is, molarcharge);
3639 printf(
" Step 2: \n");
3642 for (z1 = 1; z1 <=4; z1++) {
3643 for (z2 =1; z2 <=4; z2++) {
3644 calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]);
3647 printf(
" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n",
3648 z1, z2, etheta[z1][z2], etheta_prime[z1][z2]);
3656 printf(
" Step 3: \n");
3657 printf(
" Species Species g(x) "
3668 for (i = 1; i < (m_kk - 1); i++) {
3669 for (j = (i+1); j <
m_kk; j++) {
3678 if (charge[i]*charge[j] < 0) {
3682 x1 = sqrtIs * alpha1MX[counterIJ];
3683 if (x1 > 1.0E-100) {
3684 gfunc[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 * x1);
3685 hfunc[counterIJ] = -2.0 *
3686 (1.0-(1.0 + x1 + 0.5 * x1 *x1) * exp(-x1)) / (x1 * x1);
3688 gfunc[counterIJ] = 0.0;
3689 hfunc[counterIJ] = 0.0;
3692 if (beta2MX_L[counterIJ] != 0.0) {
3693 x2 = sqrtIs * alpha2MX[counterIJ];
3694 if (x2 > 1.0E-100) {
3695 g2func[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
3696 h2func[counterIJ] = -2.0 *
3697 (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
3699 g2func[counterIJ] = 0.0;
3700 h2func[counterIJ] = 0.0;
3704 gfunc[counterIJ] = 0.0;
3705 hfunc[counterIJ] = 0.0;
3711 printf(
" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(),
3712 gfunc[counterIJ], hfunc[counterIJ]);
3725 printf(
" Step 4: \n");
3726 printf(
" Species Species BMX "
3727 "BprimeMX BphiMX \n");
3731 for (i = 1; i < m_kk - 1; i++) {
3732 for (j = i+1; j <
m_kk; j++) {
3742 if (charge[i]*charge[j] < 0.0) {
3743 BMX_L[counterIJ] = beta0MX_L[counterIJ]
3744 + beta1MX_L[counterIJ] * gfunc[counterIJ]
3745 + beta2MX_L[counterIJ] * gfunc[counterIJ];
3748 printf(
"%d %g: %g %g %g %g\n",
3749 counterIJ, BMX_L[counterIJ], beta0MX_L[counterIJ],
3750 beta1MX_L[counterIJ], beta2MX_L[counterIJ], gfunc[counterIJ]);
3753 if (Is > 1.0E-150) {
3754 BprimeMX_L[counterIJ] = (beta1MX_L[counterIJ] * hfunc[counterIJ]/Is +
3755 beta2MX_L[counterIJ] * h2func[counterIJ]/Is);
3757 BprimeMX_L[counterIJ] = 0.0;
3759 BphiMX_L[counterIJ] = BMX_L[counterIJ] + Is*BprimeMX_L[counterIJ];
3761 BMX_L[counterIJ] = 0.0;
3762 BprimeMX_L[counterIJ] = 0.0;
3763 BphiMX_L[counterIJ] = 0.0;
3769 printf(
" %-16s %-16s %11.7f %11.7f %11.7f \n",
3770 sni.c_str(), snj.c_str(),
3771 BMX_L[counterIJ], BprimeMX_L[counterIJ], BphiMX_L[counterIJ]);
3783 printf(
" Step 5: \n");
3784 printf(
" Species Species CMX \n");
3787 for (i = 1; i < m_kk-1; i++) {
3788 for (j = i+1; j <
m_kk; j++) {
3798 if (charge[i]*charge[j] < 0.0) {
3799 CMX_L[counterIJ] = CphiMX_L[counterIJ]/
3800 (2.0* sqrt(fabs(charge[i]*charge[j])));
3802 CMX_L[counterIJ] = 0.0;
3808 printf(
" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(),
3821 printf(
" Step 6: \n");
3822 printf(
" Species Species Phi_ij "
3823 " Phiprime_ij Phi^phi_ij \n");
3826 for (i = 1; i < m_kk-1; i++) {
3827 for (j = i+1; j <
m_kk; j++) {
3837 if (charge[i]*charge[j] > 0) {
3838 z1 = (int) fabs(charge[i]);
3839 z2 = (int) fabs(charge[j]);
3841 Phi_L[counterIJ] = thetaij_L[counterIJ];
3843 Phiprime[counterIJ] = 0.0;
3844 Phiphi_L[counterIJ] = Phi_L[counterIJ] + Is * Phiprime[counterIJ];
3846 Phi_L[counterIJ] = 0.0;
3847 Phiprime[counterIJ] = 0.0;
3848 Phiphi_L[counterIJ] = 0.0;
3854 printf(
" %-16s %-16s %10.6f %10.6f %10.6f \n",
3855 sni.c_str(), snj.c_str(),
3856 Phi_L[counterIJ], Phiprime[counterIJ], Phiphi_L[counterIJ]);
3867 printf(
" Step 7: \n");
3876 double dAphidT = dA_DebyedT /3.0;
3883 dFdT = -dAphidT * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
3884 + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
3887 printf(
" initial value of dFdT = %10.6f \n", dFdT);
3890 for (i = 1; i < m_kk-1; i++) {
3891 for (j = i+1; j <
m_kk; j++) {
3901 if (charge[i]*charge[j] < 0) {
3902 dFdT = dFdT + molality[i]*molality[j] * BprimeMX_L[counterIJ];
3908 if (charge[i]*charge[j] > 0) {
3909 dFdT = dFdT + molality[i]*molality[j] * Phiprime[counterIJ];
3913 printf(
" dFdT = %10.6f \n", dFdT);
3920 printf(
" Step 8: \n");
3924 for (i = 1; i <
m_kk; i++) {
3930 if (charge[i] > 0) {
3932 zsqdFdT = charge[i]*charge[i]*dFdT;
3938 for (j = 1; j <
m_kk; j++) {
3945 if (charge[j] < 0.0) {
3947 sum1 = sum1 + molality[j]*
3948 (2.0*BMX_L[counterIJ] + molarcharge*CMX_L[counterIJ]);
3955 for (k = j+1; k <
m_kk; k++) {
3957 if (charge[k] < 0.0) {
3958 n = k + j * m_kk + i * m_kk *
m_kk;
3959 sum3 = sum3 + molality[j]*molality[k]*psi_ijk_L[n];
3966 if (charge[j] > 0.0) {
3969 sum2 = sum2 + molality[j]*(2.0*Phi_L[counterIJ]);
3971 for (k = 1; k <
m_kk; k++) {
3972 if (charge[k] < 0.0) {
3975 n = k + j * m_kk + i * m_kk *
m_kk;
3976 sum2 = sum2 + molality[j]*molality[k]*psi_ijk_L[n];
3982 sum4 = sum4 + (fabs(charge[i])*
3983 molality[j]*molality[k]*CMX_L[counterIJ2]);
3991 if (charge[j] == 0) {
3997 for (
size_t k = 1; k <
m_kk; k++) {
3998 if (charge[k] < 0.0) {
4001 n = izeta * m_kk * m_kk + jzeta * m_kk + k;
4002 double zeta_L = psi_ijk_L[n];
4003 if (zeta_L != 0.0) {
4004 sum5 = sum5 + molality[j]*molality[k]*zeta_L;
4014 zsqdFdT + sum1 + sum2 + sum3 + sum4 + sum5;
4019 printf(
" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f \n",
4021 printf(
" %12g %12g %12g %12g %12g %12g\n",
4022 zsqdFdT, sum1, sum2, sum3, sum4, sum5);
4031 if (charge[i] < 0) {
4033 zsqdFdT = charge[i]*charge[i]*dFdT;
4039 for (j = 1; j <
m_kk; j++) {
4049 if (charge[j] > 0) {
4050 sum1 = sum1 + molality[j]*
4051 (2.0*BMX_L[counterIJ] + molarcharge*CMX_L[counterIJ]);
4053 for (k = j+1; k <
m_kk; k++) {
4055 if (charge[k] > 0) {
4056 n = k + j * m_kk + i * m_kk *
m_kk;
4057 sum3 = sum3 + molality[j]*molality[k]*psi_ijk_L[n];
4066 if (charge[j] < 0.0) {
4069 sum2 = sum2 + molality[j]*(2.0*Phi_L[counterIJ]);
4071 for (k = 1; k <
m_kk; k++) {
4072 if (charge[k] > 0.0) {
4074 n = k + j * m_kk + i * m_kk *
m_kk;
4075 sum2 = sum2 + molality[j]*molality[k]*psi_ijk_L[n];
4083 molality[j]*molality[k]*CMX_L[counterIJ2]);
4091 if (charge[j] == 0.0) {
4093 for (
size_t k = 1; k <
m_kk; k++) {
4094 if (charge[k] > 0.0) {
4098 n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
4099 double zeta_L = psi_ijk_L[n];
4100 if (zeta_L != 0.0) {
4101 sum5 = sum5 + molality[j]*molality[k]*zeta_L;
4108 zsqdFdT + sum1 + sum2 + sum3 + sum4 + sum5;
4113 printf(
" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f\n",
4115 printf(
" %12g %12g %12g %12g %12g %12g\n",
4116 zsqdFdT, sum1, sum2, sum3, sum4, sum5);
4125 if (charge[i] == 0.0) {
4128 for (j = 1; j <
m_kk; j++) {
4133 if (charge[j] > 0.0) {
4134 for (k = 1; k <
m_kk; k++) {
4135 if (charge[k] < 0.0) {
4136 n = k + j * m_kk + i * m_kk *
m_kk;
4137 sum3 = sum3 + molality[j]*molality[k]*psi_ijk_L[n];
4142 sum2 = 3.0 * molality[i] * molality[i] *
m_Mu_nnn_L[i];
4148 printf(
" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f \n",
4157 printf(
" Step 9: \n");
4179 term1 = -dAphidT * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
4181 for (j = 1; j <
m_kk; j++) {
4185 if (charge[j] > 0.0) {
4186 for (k = 1; k <
m_kk; k++) {
4187 if (charge[k] < 0.0) {
4194 sum1 = sum1 + molality[j]*molality[k]*
4195 (BphiMX_L[counterIJ] + molarcharge*CMX_L[counterIJ]);
4199 for (k = j+1; k <
m_kk; k++) {
4200 if (j == (m_kk-1)) {
4202 printf(
"logic error 1 in Step 9 of hmw_act");
4205 if (charge[k] > 0.0) {
4212 sum2 = sum2 + molality[j]*molality[k]*Phiphi_L[counterIJ];
4213 for (m = 1; m <
m_kk; m++) {
4214 if (charge[m] < 0.0) {
4216 n = m + k * m_kk + j * m_kk *
m_kk;
4218 molality[j]*molality[k]*molality[m]*psi_ijk_L[n];
4228 if (charge[j] < 0) {
4229 for (k = j+1; k <
m_kk; k++) {
4232 printf(
"logic error 2 in Step 9 of hmw_act");
4235 if (charge[k] < 0) {
4243 sum3 = sum3 + molality[j]*molality[k]*Phiphi_L[counterIJ];
4244 for (m = 1; m <
m_kk; m++) {
4245 if (charge[m] > 0.0) {
4246 n = m + k * m_kk + j * m_kk *
m_kk;
4248 molality[j]*molality[k]*molality[m]*psi_ijk_L[n];
4258 if (charge[j] == 0) {
4259 for (k = 1; k <
m_kk; k++) {
4260 if (charge[k] < 0.0) {
4263 if (charge[k] > 0.0) {
4266 if (charge[k] == 0.0) {
4269 }
else if (k == j) {
4270 sum6 = sum6 + 0.5 * molality[j]*molality[k]*
m_Lambda_nj_L(j,k);
4273 if (charge[k] < 0.0) {
4275 for (m = 1; m <
m_kk; m++) {
4276 if (charge[m] > 0.0) {
4278 n = k + jzeta * m_kk + izeta * m_kk *
m_kk;
4279 double zeta_L = psi_ijk_L[n];
4280 if (zeta_L != 0.0) {
4281 sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_L;
4287 sum7 += molality[j]*molality[j]*molality[j]*
m_Mu_nnn_L[j];
4290 sum_m_phi_minus_1 = 2.0 *
4291 (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
4296 if (molalitysum > 1.0E-150) {
4297 d_osmotic_coef_dT = 0.0 + (sum_m_phi_minus_1 / molalitysum);
4299 d_osmotic_coef_dT = 0.0;
4303 printf(
" term1=%10.6f sum1=%10.6f sum2=%10.6f "
4304 "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n",
4305 term1, sum1, sum2, sum3, sum4, sum5);
4306 printf(
" sum_m_phi_minus_1=%10.6f d_osmotic_coef_dT =%10.6f\n",
4307 sum_m_phi_minus_1, d_osmotic_coef_dT);
4311 printf(
" Step 10: \n");
4314 d_lnwateract_dT = -(
m_weightSolvent/1000.0) * molalitysum * d_osmotic_coef_dT;
4328 double d_wateract_dT = exp(d_lnwateract_dT);
4329 printf(
" d_ln_a_water_dT = %10.6f d_a_water_dT=%10.6f\n\n",
4330 d_lnwateract_dT, d_wateract_dT);
4355 for (
size_t k = 1; k <
m_kk; k++) {
4403 printf(
"Wrong index solvent value!\n");
4407 std::string sni, snj, snk;
4423 double etheta[5][5], etheta_prime[5][5], sqrtIs;
4432 double molarcharge = 0.0;
4437 double molalitysum = 0.0;
4453 double d2FdT2, zsqd2FdT2;
4454 double sum1, sum2, sum3, sum4, sum5, term1;
4455 double sum_m_phi_minus_1, d2_osmotic_coef_dT2, d2_lnwateract_dT2;
4458 size_t n, i, j, k, m, counterIJ, counterIJ2;
4462 printf(
"\n Debugging information from "
4463 "s_Pitzer_d2lnMolalityActCoeff_dT2()\n");
4475 for (n = 1; n <
m_kk; n++) {
4477 Is += charge[n] * charge[n] * molality[n];
4479 molarcharge += fabs(charge[n]) * molality[n];
4480 molalitysum += molality[n];
4491 printf(
" Step 1: \n");
4492 printf(
" ionic strenth = %14.7le \n total molar "
4493 "charge = %14.7le \n", Is, molarcharge);
4511 printf(
" Step 2: \n");
4514 for (z1 = 1; z1 <=4; z1++) {
4515 for (z2 =1; z2 <=4; z2++) {
4516 calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]);
4519 printf(
" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n",
4520 z1, z2, etheta[z1][z2], etheta_prime[z1][z2]);
4528 printf(
" Step 3: \n");
4529 printf(
" Species Species g(x) "
4540 for (i = 1; i < (m_kk - 1); i++) {
4541 for (j = (i+1); j <
m_kk; j++) {
4550 if (charge[i]*charge[j] < 0) {
4554 x1 = sqrtIs * alpha1MX[counterIJ];
4555 if (x1 > 1.0E-100) {
4556 gfunc[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 *x1);
4557 hfunc[counterIJ] = -2.0*
4558 (1.0-(1.0 + x1 + 0.5*x1 * x1) * exp(-x1)) / (x1 * x1);
4560 gfunc[counterIJ] = 0.0;
4561 hfunc[counterIJ] = 0.0;
4564 if (beta2MX_LL[counterIJ] != 0.0) {
4565 x2 = sqrtIs * alpha2MX[counterIJ];
4566 if (x2 > 1.0E-100) {
4567 g2func[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
4568 h2func[counterIJ] = -2.0 *
4569 (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
4571 g2func[counterIJ] = 0.0;
4572 h2func[counterIJ] = 0.0;
4576 gfunc[counterIJ] = 0.0;
4577 hfunc[counterIJ] = 0.0;
4583 printf(
" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(),
4584 gfunc[counterIJ], hfunc[counterIJ]);
4596 printf(
" Step 4: \n");
4597 printf(
" Species Species BMX "
4598 "BprimeMX BphiMX \n");
4602 for (i = 1; i < m_kk - 1; i++) {
4603 for (j = i+1; j <
m_kk; j++) {
4613 if (charge[i]*charge[j] < 0.0) {
4614 BMX_LL[counterIJ] = beta0MX_LL[counterIJ]
4615 + beta1MX_LL[counterIJ] * gfunc[counterIJ]
4616 + beta2MX_LL[counterIJ] * g2func[counterIJ];
4619 printf(
"%d %g: %g %g %g %g\n",
4620 counterIJ, BMX_LL[counterIJ], beta0MX_LL[counterIJ],
4621 beta1MX_LL[counterIJ], beta2MX_LL[counterIJ], gfunc[counterIJ]);
4624 if (Is > 1.0E-150) {
4625 BprimeMX_LL[counterIJ] = (beta1MX_LL[counterIJ] * hfunc[counterIJ]/Is +
4626 beta2MX_LL[counterIJ] * h2func[counterIJ]/Is);
4628 BprimeMX_LL[counterIJ] = 0.0;
4630 BphiMX_LL[counterIJ] = BMX_LL[counterIJ] + Is*BprimeMX_LL[counterIJ];
4632 BMX_LL[counterIJ] = 0.0;
4633 BprimeMX_LL[counterIJ] = 0.0;
4634 BphiMX_LL[counterIJ] = 0.0;
4640 printf(
" %-16s %-16s %11.7f %11.7f %11.7f \n",
4641 sni.c_str(), snj.c_str(),
4642 BMX_LL[counterIJ], BprimeMX_LL[counterIJ], BphiMX_LL[counterIJ]);
4654 printf(
" Step 5: \n");
4655 printf(
" Species Species CMX \n");
4658 for (i = 1; i < m_kk-1; i++) {
4659 for (j = i+1; j <
m_kk; j++) {
4669 if (charge[i]*charge[j] < 0.0) {
4670 CMX_LL[counterIJ] = CphiMX_LL[counterIJ]/
4671 (2.0* sqrt(fabs(charge[i]*charge[j])));
4673 CMX_LL[counterIJ] = 0.0;
4679 printf(
" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(),
4692 printf(
" Step 6: \n");
4693 printf(
" Species Species Phi_ij "
4694 " Phiprime_ij Phi^phi_ij \n");
4697 for (i = 1; i < m_kk-1; i++) {
4698 for (j = i+1; j <
m_kk; j++) {
4708 if (charge[i]*charge[j] > 0) {
4709 z1 = (int) fabs(charge[i]);
4710 z2 = (int) fabs(charge[j]);
4713 Phi_LL[counterIJ] = thetaij_LL[counterIJ];
4715 Phiprime[counterIJ] = 0.0;
4718 Phiphi_LL[counterIJ] = Phi_LL[counterIJ];
4720 Phi_LL[counterIJ] = 0.0;
4721 Phiprime[counterIJ] = 0.0;
4722 Phiphi_LL[counterIJ] = 0.0;
4731 printf(
" %-16s %-16s %10.6f %10.6f %10.6f \n",
4732 sni.c_str(), snj.c_str(),
4733 Phi_LL[counterIJ], Phiprime[counterIJ], Phiphi_LL[counterIJ]);
4744 printf(
" Step 7: \n");
4764 d2FdT2 = -d2AphidT2 * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
4765 + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
4768 printf(
" initial value of d2FdT2 = %10.6f \n", d2FdT2);
4771 for (i = 1; i < m_kk-1; i++) {
4772 for (j = i+1; j <
m_kk; j++) {
4782 if (charge[i]*charge[j] < 0) {
4783 d2FdT2 = d2FdT2 + molality[i]*molality[j] * BprimeMX_LL[counterIJ];
4789 if (charge[i]*charge[j] > 0) {
4790 d2FdT2 = d2FdT2 + molality[i]*molality[j] * Phiprime[counterIJ];
4794 printf(
" d2FdT2 = %10.6f \n", d2FdT2);
4801 printf(
" Step 8: \n");
4805 for (i = 1; i <
m_kk; i++) {
4811 if (charge[i] > 0) {
4813 zsqd2FdT2 = charge[i]*charge[i]*d2FdT2;
4819 for (j = 1; j <
m_kk; j++) {
4826 if (charge[j] < 0.0) {
4828 sum1 = sum1 + molality[j]*
4829 (2.0*BMX_LL[counterIJ] + molarcharge*CMX_LL[counterIJ]);
4836 for (k = j+1; k <
m_kk; k++) {
4838 if (charge[k] < 0.0) {
4839 n = k + j * m_kk + i * m_kk *
m_kk;
4840 sum3 = sum3 + molality[j]*molality[k]*psi_ijk_LL[n];
4847 if (charge[j] > 0.0) {
4850 sum2 = sum2 + molality[j]*(2.0*Phi_LL[counterIJ]);
4852 for (k = 1; k <
m_kk; k++) {
4853 if (charge[k] < 0.0) {
4856 n = k + j * m_kk + i * m_kk *
m_kk;
4857 sum2 = sum2 + molality[j]*molality[k]*psi_ijk_LL[n];
4863 sum4 = sum4 + (fabs(charge[i])*
4864 molality[j]*molality[k]*CMX_LL[counterIJ2]);
4872 if (charge[j] == 0) {
4877 for (
size_t k = 1; k <
m_kk; k++) {
4878 if (charge[k] < 0.0) {
4881 n = izeta * m_kk * m_kk + jzeta * m_kk + k;
4882 double zeta_LL = psi_ijk_LL[n];
4883 if (zeta_LL != 0.0) {
4884 sum5 = sum5 + molality[j]*molality[k]*zeta_LL;
4895 zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5;
4899 printf(
" %-16s d2lngammadT2[i]=%10.6f \n",
4901 printf(
" %12g %12g %12g %12g %12g %12g\n",
4902 zsqd2FdT2, sum1, sum2, sum3, sum4, sum5);
4912 if (charge[i] < 0) {
4914 zsqd2FdT2 = charge[i]*charge[i]*d2FdT2;
4920 for (j = 1; j <
m_kk; j++) {
4930 if (charge[j] > 0) {
4931 sum1 = sum1 + molality[j]*
4932 (2.0*BMX_LL[counterIJ] + molarcharge*CMX_LL[counterIJ]);
4934 for (k = j+1; k <
m_kk; k++) {
4936 if (charge[k] > 0) {
4937 n = k + j * m_kk + i * m_kk *
m_kk;
4938 sum3 = sum3 + molality[j]*molality[k]*psi_ijk_LL[n];
4947 if (charge[j] < 0.0) {
4950 sum2 = sum2 + molality[j]*(2.0*Phi_LL[counterIJ]);
4952 for (k = 1; k <
m_kk; k++) {
4953 if (charge[k] > 0.0) {
4955 n = k + j * m_kk + i * m_kk *
m_kk;
4956 sum2 = sum2 + molality[j]*molality[k]*psi_ijk_LL[n];
4964 molality[j]*molality[k]*CMX_LL[counterIJ2]);
4972 if (charge[j] == 0.0) {
4977 for (
size_t k = 1; k <
m_kk; k++) {
4978 if (charge[k] > 0.0) {
4982 n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
4983 double zeta_LL = psi_ijk_LL[n];
4984 if (zeta_LL != 0.0) {
4985 sum5 = sum5 + molality[j]*molality[k]*zeta_LL;
4992 zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5;
4996 printf(
" %-16s d2lngammadT2[i]=%10.6f\n",
4998 printf(
" %12g %12g %12g %12g %12g %12g\n",
4999 zsqd2FdT2, sum1, sum2, sum3, sum4, sum5);
5008 if (charge[i] == 0.0) {
5011 for (j = 1; j <
m_kk; j++) {
5016 if (charge[j] > 0.0) {
5017 for (k = 1; k <
m_kk; k++) {
5018 if (charge[k] < 0.0) {
5019 n = k + j * m_kk + i * m_kk *
m_kk;
5020 sum3 = sum3 + molality[j]*molality[k]*psi_ijk_LL[n];
5025 sum2 = 3.0 * molality[i] * molality[i] *
m_Mu_nnn_LL[i];
5030 printf(
" %-16s d2lngammadT2[i]=%10.6f \n",
5039 printf(
" Step 9: \n");
5062 term1 = -d2AphidT2 * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
5064 for (j = 1; j <
m_kk; j++) {
5068 if (charge[j] > 0.0) {
5069 for (k = 1; k <
m_kk; k++) {
5070 if (charge[k] < 0.0) {
5077 sum1 = sum1 + molality[j]*molality[k]*
5078 (BphiMX_LL[counterIJ] + molarcharge*CMX_LL[counterIJ]);
5082 for (k = j+1; k <
m_kk; k++) {
5083 if (j == (m_kk-1)) {
5085 printf(
"logic error 1 in Step 9 of hmw_act");
5088 if (charge[k] > 0.0) {
5095 sum2 = sum2 + molality[j]*molality[k]*Phiphi_LL[counterIJ];
5096 for (m = 1; m <
m_kk; m++) {
5097 if (charge[m] < 0.0) {
5099 n = m + k * m_kk + j * m_kk *
m_kk;
5101 molality[j]*molality[k]*molality[m]*psi_ijk_LL[n];
5111 if (charge[j] < 0) {
5112 for (k = j+1; k <
m_kk; k++) {
5115 printf(
"logic error 2 in Step 9 of hmw_act");
5118 if (charge[k] < 0) {
5126 sum3 = sum3 + molality[j]*molality[k]*Phiphi_LL[counterIJ];
5127 for (m = 1; m <
m_kk; m++) {
5128 if (charge[m] > 0.0) {
5129 n = m + k * m_kk + j * m_kk *
m_kk;
5131 molality[j]*molality[k]*molality[m]*psi_ijk_LL[n];
5141 if (charge[j] == 0) {
5142 for (k = 1; k <
m_kk; k++) {
5143 if (charge[k] < 0.0) {
5146 if (charge[k] > 0.0) {
5149 if (charge[k] == 0.0) {
5152 }
else if (k == j) {
5156 if (charge[k] < 0.0) {
5158 for (m = 1; m <
m_kk; m++) {
5159 if (charge[m] > 0.0) {
5161 n = k + jzeta * m_kk + izeta * m_kk *
m_kk;
5162 double zeta_LL = psi_ijk_LL[n];
5163 if (zeta_LL != 0.0) {
5164 sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_LL;
5171 sum7 += molality[j] * molality[j] * molality[j] *
m_Mu_nnn_LL[j];
5174 sum_m_phi_minus_1 = 2.0 *
5175 (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
5180 if (molalitysum > 1.0E-150) {
5181 d2_osmotic_coef_dT2 = 0.0 + (sum_m_phi_minus_1 / molalitysum);
5183 d2_osmotic_coef_dT2 = 0.0;
5187 printf(
" term1=%10.6f sum1=%10.6f sum2=%10.6f "
5188 "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n",
5189 term1, sum1, sum2, sum3, sum4, sum5);
5190 printf(
" sum_m_phi_minus_1=%10.6f d2_osmotic_coef_dT2=%10.6f\n",
5191 sum_m_phi_minus_1, d2_osmotic_coef_dT2);
5195 printf(
" Step 10: \n");
5198 d2_lnwateract_dT2 = -(
m_weightSolvent/1000.0) * molalitysum * d2_osmotic_coef_dT2;
5212 double d2_wateract_dT2 = exp(d2_lnwateract_dT2);
5213 printf(
" d2_ln_a_water_dT2 = %10.6f d2_a_water_dT2=%10.6f\n\n",
5214 d2_lnwateract_dT2, d2_wateract_dT2);
5239 for (
size_t k = 1; k <
m_kk; k++) {
5282 printf(
"Wrong index solvent value!\n");
5286 std::string sni, snj, snk;
5302 double etheta[5][5], etheta_prime[5][5], sqrtIs;
5311 double molarcharge = 0.0;
5316 double molalitysum = 0.0;
5331 double dFdP, zsqdFdP;
5332 double sum1, sum2, sum3, sum4, sum5, term1;
5333 double sum_m_phi_minus_1, d_osmotic_coef_dP, d_lnwateract_dP;
5336 size_t n, i, j, k, m, counterIJ, counterIJ2;
5343 printf(
"\n Debugging information from "
5344 "s_Pitzer_dlnMolalityActCoeff_dP()\n");
5355 for (n = 1; n <
m_kk; n++) {
5357 Is += charge[n] * charge[n] * molality[n];
5359 molarcharge += fabs(charge[n]) * molality[n];
5360 molalitysum += molality[n];
5371 printf(
" Step 1: \n");
5372 printf(
" ionic strenth = %14.7le \n total molar "
5373 "charge = %14.7le \n", Is, molarcharge);
5392 printf(
" Step 2: \n");
5395 for (z1 = 1; z1 <=4; z1++) {
5396 for (z2 =1; z2 <=4; z2++) {
5397 calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]);
5400 printf(
" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n",
5401 z1, z2, etheta[z1][z2], etheta_prime[z1][z2]);
5409 printf(
" Step 3: \n");
5410 printf(
" Species Species g(x) "
5421 for (i = 1; i < (m_kk - 1); i++) {
5422 for (j = (i+1); j <
m_kk; j++) {
5431 if (charge[i]*charge[j] < 0) {
5435 x1 = sqrtIs * alpha1MX[counterIJ];
5436 if (x1 > 1.0E-100) {
5437 gfunc[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 * x1);
5438 hfunc[counterIJ] = -2.0*
5439 (1.0-(1.0 + x1 + 0.5 * x1 * x1) * exp(-x1)) / (x1 * x1);
5441 gfunc[counterIJ] = 0.0;
5442 hfunc[counterIJ] = 0.0;
5445 if (beta2MX_P[counterIJ] != 0.0) {
5446 x2 = sqrtIs * alpha2MX[counterIJ];
5447 if (x2 > 1.0E-100) {
5448 g2func[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
5449 h2func[counterIJ] = -2.0 *
5450 (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
5452 g2func[counterIJ] = 0.0;
5453 h2func[counterIJ] = 0.0;
5457 gfunc[counterIJ] = 0.0;
5458 hfunc[counterIJ] = 0.0;
5464 printf(
" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(),
5465 gfunc[counterIJ], hfunc[counterIJ]);
5478 printf(
" Step 4: \n");
5479 printf(
" Species Species BMX "
5480 "BprimeMX BphiMX \n");
5484 for (i = 1; i < m_kk - 1; i++) {
5485 for (j = i+1; j <
m_kk; j++) {
5495 if (charge[i]*charge[j] < 0.0) {
5496 BMX_P[counterIJ] = beta0MX_P[counterIJ]
5497 + beta1MX_P[counterIJ] * gfunc[counterIJ]
5498 + beta2MX_P[counterIJ] * g2func[counterIJ];
5501 printf(
"%d %g: %g %g %g %g\n",
5502 counterIJ, BMX_P[counterIJ], beta0MX_P[counterIJ],
5503 beta1MX_P[counterIJ], beta2MX_P[counterIJ], gfunc[counterIJ]);
5506 if (Is > 1.0E-150) {
5507 BprimeMX_P[counterIJ] = (beta1MX_P[counterIJ] * hfunc[counterIJ]/Is +
5508 beta2MX_P[counterIJ] * h2func[counterIJ]/Is);
5510 BprimeMX_P[counterIJ] = 0.0;
5512 BphiMX_P[counterIJ] = BMX_P[counterIJ] + Is*BprimeMX_P[counterIJ];
5514 BMX_P[counterIJ] = 0.0;
5515 BprimeMX_P[counterIJ] = 0.0;
5516 BphiMX_P[counterIJ] = 0.0;
5522 printf(
" %-16s %-16s %11.7f %11.7f %11.7f \n",
5523 sni.c_str(), snj.c_str(),
5524 BMX_P[counterIJ], BprimeMX_P[counterIJ], BphiMX_P[counterIJ]);
5536 printf(
" Step 5: \n");
5537 printf(
" Species Species CMX \n");
5540 for (i = 1; i < m_kk-1; i++) {
5541 for (j = i+1; j <
m_kk; j++) {
5551 if (charge[i]*charge[j] < 0.0) {
5552 CMX_P[counterIJ] = CphiMX_P[counterIJ]/
5553 (2.0* sqrt(fabs(charge[i]*charge[j])));
5555 CMX_P[counterIJ] = 0.0;
5561 printf(
" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(),
5574 printf(
" Step 6: \n");
5575 printf(
" Species Species Phi_ij "
5576 " Phiprime_ij Phi^phi_ij \n");
5579 for (i = 1; i < m_kk-1; i++) {
5580 for (j = i+1; j <
m_kk; j++) {
5590 if (charge[i]*charge[j] > 0) {
5591 z1 = (int) fabs(charge[i]);
5592 z2 = (int) fabs(charge[j]);
5594 Phi_P[counterIJ] = thetaij_P[counterIJ];
5596 Phiprime[counterIJ] = 0.0;
5597 Phiphi_P[counterIJ] = Phi_P[counterIJ] + Is * Phiprime[counterIJ];
5599 Phi_P[counterIJ] = 0.0;
5600 Phiprime[counterIJ] = 0.0;
5601 Phiphi_P[counterIJ] = 0.0;
5607 printf(
" %-16s %-16s %10.6f %10.6f %10.6f \n",
5608 sni.c_str(), snj.c_str(),
5609 Phi_P[counterIJ], Phiprime[counterIJ], Phiphi_P[counterIJ]);
5620 printf(
" Step 7: \n");
5629 double dAphidP = dA_DebyedP /3.0;
5636 dFdP = -dAphidP * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
5637 + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
5640 printf(
" initial value of dFdP = %10.6f \n", dFdP);
5643 for (i = 1; i < m_kk-1; i++) {
5644 for (j = i+1; j <
m_kk; j++) {
5654 if (charge[i]*charge[j] < 0) {
5655 dFdP = dFdP + molality[i]*molality[j] * BprimeMX_P[counterIJ];
5661 if (charge[i]*charge[j] > 0) {
5662 dFdP = dFdP + molality[i]*molality[j] * Phiprime[counterIJ];
5666 printf(
" dFdP = %10.6f \n", dFdP);
5673 printf(
" Step 8: \n");
5678 for (i = 1; i <
m_kk; i++) {
5684 if (charge[i] > 0) {
5686 zsqdFdP = charge[i]*charge[i]*dFdP;
5692 for (j = 1; j <
m_kk; j++) {
5699 if (charge[j] < 0.0) {
5701 sum1 = sum1 + molality[j]*
5702 (2.0*BMX_P[counterIJ] + molarcharge*CMX_P[counterIJ]);
5709 for (k = j+1; k <
m_kk; k++) {
5711 if (charge[k] < 0.0) {
5712 n = k + j * m_kk + i * m_kk *
m_kk;
5713 sum3 = sum3 + molality[j]*molality[k]*psi_ijk_P[n];
5720 if (charge[j] > 0.0) {
5723 sum2 = sum2 + molality[j]*(2.0*Phi_P[counterIJ]);
5725 for (k = 1; k <
m_kk; k++) {
5726 if (charge[k] < 0.0) {
5729 n = k + j * m_kk + i * m_kk *
m_kk;
5730 sum2 = sum2 + molality[j]*molality[k]*psi_ijk_P[n];
5736 sum4 = sum4 + (fabs(charge[i])*
5737 molality[j]*molality[k]*CMX_P[counterIJ2]);