Cantera  2.0
HMWSoln.cpp
Go to the documentation of this file.
1 /**
2  * @file HMWSoln.cpp
3  * Definitions for the %HMWSoln ThermoPhase object, which
4  * models concentrated electrolyte solutions
5  * (see \ref thermoprops and \link Cantera::HMWSoln HMWSoln \endlink) .
6  *
7  * Class %HMWSoln represents a concentrated liquid electrolyte phase which
8  * obeys the Pitzer formulation for nonideality using molality-based
9  * standard states.
10  *
11  * This version of the code was modified to have the binary Beta2 Pitzer
12  * parameter consistent with the temperature expansions used for Beta0,
13  * Beta1, and Cphi.(CFJC, SNL)
14  */
15 /*
16  * Copyright (2006) Sandia Corporation. Under the terms of
17  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
18  * U.S. Government retains certain rights in this software.
19  */
20 
21 #include "cantera/thermo/HMWSoln.h"
22 #include "cantera/thermo/ThermoFactory.h"
23 #include "cantera/thermo/WaterProps.h"
24 #include "cantera/thermo/PDSS_Water.h"
25 #include "cantera/base/stringUtils.h"
26 
27 #include <cmath>
28 #include <cstdlib>
29 #include <cstdio>
30 
31 namespace Cantera
32 {
33 
34 /*
35  * Default constructor
36  */
37 HMWSoln::HMWSoln() :
38  MolalityVPSSTP(),
39  m_formPitzer(PITZERFORM_BASE),
40  m_formPitzerTemp(PITZER_TEMP_CONSTANT),
41  m_formGC(2),
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),
47  m_A_Debye(1.172576), // units = sqrt(kg/gmol)
48  m_waterSS(0),
49  m_densWaterSS(1000.),
50  m_waterProps(0),
51  m_molalitiesAreCropped(false),
52  IMS_typeCutoff_(0),
53  IMS_X_o_cutoff_(0.2),
54  IMS_gamma_o_min_(1.0E-5),
55  IMS_gamma_k_min_(10.0),
56  IMS_cCut_(0.05),
57  IMS_slopefCut_(0.6),
58  IMS_dfCut_(0.0),
59  IMS_efCut_(0.0),
60  IMS_afCut_(0.0),
61  IMS_bfCut_(0.0),
62  IMS_slopegCut_(0.0),
63  IMS_dgCut_(0.0),
64  IMS_egCut_(0.0),
65  IMS_agCut_(0.0),
66  IMS_bgCut_(0.0),
67  MC_X_o_cutoff_(0.0),
68  MC_X_o_min_(0.0),
69  MC_slopepCut_(0.0),
70  MC_dpCut_(0.0),
71  MC_epCut_(0.0),
72  MC_apCut_(0.0),
73  MC_bpCut_(0.0),
74  MC_cpCut_(0.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),
79  m_debugCalc(0)
80 {
81  for (size_t i = 0; i < 17; i++) {
82  elambda[i] = 0.0;
83  elambda1[i] = 0.0;
84  }
85 }
86 /*
87  * Working constructors
88  *
89  * The two constructors below are the normal way
90  * the phase initializes itself. They are shells that call
91  * the routine initThermo(), with a reference to the
92  * XML database to get the info for the phase.
93  */
94 HMWSoln::HMWSoln(std::string inputFile, std::string id) :
95  MolalityVPSSTP(),
96  m_formPitzer(PITZERFORM_BASE),
97  m_formPitzerTemp(PITZER_TEMP_CONSTANT),
98  m_formGC(2),
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),
104  m_A_Debye(1.172576), // units = sqrt(kg/gmol)
105  m_waterSS(0),
106  m_densWaterSS(1000.),
107  m_waterProps(0),
108  m_molalitiesAreCropped(false),
109  IMS_typeCutoff_(0),
110  IMS_X_o_cutoff_(0.2),
111  IMS_gamma_o_min_(1.0E-5),
112  IMS_gamma_k_min_(10.0),
113  IMS_cCut_(0.05),
114  IMS_slopefCut_(0.6),
115  IMS_dfCut_(0.0),
116  IMS_efCut_(0.0),
117  IMS_afCut_(0.0),
118  IMS_bfCut_(0.0),
119  IMS_slopegCut_(0.0),
120  IMS_dgCut_(0.0),
121  IMS_egCut_(0.0),
122  IMS_agCut_(0.0),
123  IMS_bgCut_(0.0),
124  MC_X_o_cutoff_(0.0),
125  MC_X_o_min_(0.0),
126  MC_slopepCut_(0.0),
127  MC_dpCut_(0.0),
128  MC_epCut_(0.0),
129  MC_apCut_(0.0),
130  MC_bpCut_(0.0),
131  MC_cpCut_(0.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),
136  m_debugCalc(0)
137 {
138  for (int i = 0; i < 17; i++) {
139  elambda[i] = 0.0;
140  elambda1[i] = 0.0;
141  }
142  constructPhaseFile(inputFile, id);
143 }
144 
145 HMWSoln::HMWSoln(XML_Node& phaseRoot, std::string id) :
146  MolalityVPSSTP(),
147  m_formPitzer(PITZERFORM_BASE),
148  m_formPitzerTemp(PITZER_TEMP_CONSTANT),
149  m_formGC(2),
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),
155  m_A_Debye(1.172576), // units = sqrt(kg/gmol)
156  m_waterSS(0),
157  m_densWaterSS(1000.),
158  m_waterProps(0),
159  m_molalitiesAreCropped(false),
160  IMS_typeCutoff_(0),
161  IMS_X_o_cutoff_(0.2),
162  IMS_gamma_o_min_(1.0E-5),
163  IMS_gamma_k_min_(10.0),
164  IMS_cCut_(0.05),
165  IMS_slopefCut_(0.6),
166  IMS_dfCut_(0.0),
167  IMS_efCut_(0.0),
168  IMS_afCut_(0.0),
169  IMS_bfCut_(0.0),
170  IMS_slopegCut_(0.0),
171  IMS_dgCut_(0.0),
172  IMS_egCut_(0.0),
173  IMS_agCut_(0.0),
174  IMS_bgCut_(0.0),
175  MC_X_o_cutoff_(0.0),
176  MC_X_o_min_(0.0),
177  MC_slopepCut_(0.0),
178  MC_dpCut_(0.0),
179  MC_epCut_(0.0),
180  MC_apCut_(0.0),
181  MC_bpCut_(0.0),
182  MC_cpCut_(0.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),
187  m_debugCalc(0)
188 {
189  for (int i = 0; i < 17; i++) {
190  elambda[i] = 0.0;
191  elambda1[i] = 0.0;
192  }
193  constructPhaseXML(phaseRoot, id);
194 }
195 
196 /*
197  * Copy Constructor:
198  *
199  * Note this stuff will not work until the underlying phase
200  * has a working copy constructor
201  */
202 HMWSoln::HMWSoln(const HMWSoln& b) :
203  MolalityVPSSTP(),
204  m_formPitzer(PITZERFORM_BASE),
205  m_formPitzerTemp(PITZER_TEMP_CONSTANT),
206  m_formGC(2),
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),
212  m_A_Debye(1.172576), // units = sqrt(kg/gmol)
213  m_waterSS(0),
214  m_densWaterSS(1000.),
215  m_waterProps(0),
216  m_molalitiesAreCropped(false),
217  IMS_typeCutoff_(0),
218  IMS_X_o_cutoff_(0.2),
219  IMS_gamma_o_min_(1.0E-5),
220  IMS_gamma_k_min_(10.0),
221  IMS_cCut_(0.05),
222  IMS_slopefCut_(0.6),
223  IMS_dfCut_(0.0),
224  IMS_efCut_(0.0),
225  IMS_afCut_(0.0),
226  IMS_bfCut_(0.0),
227  IMS_slopegCut_(0.0),
228  IMS_dgCut_(0.0),
229  IMS_egCut_(0.0),
230  IMS_agCut_(0.0),
231  IMS_bgCut_(0.0),
232  MC_X_o_cutoff_(0.0),
233  MC_X_o_min_(0.0),
234  MC_slopepCut_(0.0),
235  MC_dpCut_(0.0),
236  MC_epCut_(0.0),
237  MC_apCut_(0.0),
238  MC_bpCut_(0.0),
239  MC_cpCut_(0.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),
244  m_debugCalc(0)
245 {
246  /*
247  * Use the assignment operator to do the brunt
248  * of the work for the copy constructor.
249  */
250  *this = b;
251 }
252 
253 /*
254  * operator=()
255  *
256  * Note this stuff will not work until the underlying phase
257  * has a working assignment operator
258  */
259 HMWSoln& HMWSoln::
260 operator=(const HMWSoln& b)
261 {
262  if (&b != this) {
263  MolalityVPSSTP::operator=(b);
264  m_formPitzer = b.m_formPitzer;
265  m_formPitzerTemp = b.m_formPitzerTemp;
266  m_formGC = b.m_formGC;
267  m_Aionic = b.m_Aionic;
268  m_IionicMolality = b.m_IionicMolality;
269  m_maxIionicStrength = b.m_maxIionicStrength;
270  m_TempPitzerRef = b.m_TempPitzerRef;
271  m_IionicMolalityStoich= b.m_IionicMolalityStoich;
272  m_form_A_Debye = b.m_form_A_Debye;
273  m_A_Debye = b.m_A_Debye;
274 
275  // This is an internal shallow copy of the PDSS_Water pointer
276  m_waterSS = providePDSS(0);
277  if (!m_waterSS) {
278  throw CanteraError("HMWSoln::operator=()", "Dynamic cast to PDSS_Water failed");
279  }
280 
281  m_densWaterSS = b.m_densWaterSS;
282 
283  if (m_waterProps) {
284  delete m_waterProps;
285  m_waterProps = 0;
286  }
287  if (b.m_waterProps) {
288  m_waterProps = new WaterProps(dynamic_cast<PDSS_Water*>(m_waterSS));
289  }
290 
291  m_expg0_RT = b.m_expg0_RT;
292  m_pe = b.m_pe;
293  m_pp = b.m_pp;
294  m_tmpV = b.m_tmpV;
295  m_speciesCharge_Stoich= b.m_speciesCharge_Stoich;
296  m_Beta0MX_ij = b.m_Beta0MX_ij;
297  m_Beta0MX_ij_L = b.m_Beta0MX_ij_L;
298  m_Beta0MX_ij_LL = b.m_Beta0MX_ij_LL;
299  m_Beta0MX_ij_P = b.m_Beta0MX_ij_P;
300  m_Beta0MX_ij_coeff = b.m_Beta0MX_ij_coeff;
301  m_Beta1MX_ij = b.m_Beta1MX_ij;
302  m_Beta1MX_ij_L = b.m_Beta1MX_ij_L;
303  m_Beta1MX_ij_LL = b.m_Beta1MX_ij_LL;
304  m_Beta1MX_ij_P = b.m_Beta1MX_ij_P;
305  m_Beta1MX_ij_coeff = b.m_Beta1MX_ij_coeff;
306  m_Beta2MX_ij = b.m_Beta2MX_ij;
307  m_Beta2MX_ij_L = b.m_Beta2MX_ij_L;
308  m_Beta2MX_ij_LL = b.m_Beta2MX_ij_LL;
309  m_Beta2MX_ij_P = b.m_Beta2MX_ij_P;
310  m_Beta2MX_ij_coeff = b.m_Beta2MX_ij_coeff;
311  m_Alpha1MX_ij = b.m_Alpha1MX_ij;
312  m_Alpha2MX_ij = b.m_Alpha2MX_ij;
313  m_CphiMX_ij = b.m_CphiMX_ij;
314  m_CphiMX_ij_L = b.m_CphiMX_ij_L;
315  m_CphiMX_ij_LL = b.m_CphiMX_ij_LL;
316  m_CphiMX_ij_P = b.m_CphiMX_ij_P;
317  m_CphiMX_ij_coeff = b.m_CphiMX_ij_coeff;
318  m_Theta_ij = b.m_Theta_ij;
319  m_Theta_ij_L = b.m_Theta_ij_L;
320  m_Theta_ij_LL = b.m_Theta_ij_LL;
321  m_Theta_ij_P = b.m_Theta_ij_P;
322  m_Theta_ij_coeff = b.m_Theta_ij_coeff;
323  m_Psi_ijk = b.m_Psi_ijk;
324  m_Psi_ijk_L = b.m_Psi_ijk_L;
325  m_Psi_ijk_LL = b.m_Psi_ijk_LL;
326  m_Psi_ijk_P = b.m_Psi_ijk_P;
327  m_Psi_ijk_coeff = b.m_Psi_ijk_coeff;
328  m_Lambda_nj = b.m_Lambda_nj;
329  m_Lambda_nj_L = b.m_Lambda_nj_L;
330  m_Lambda_nj_LL = b.m_Lambda_nj_LL;
331  m_Lambda_nj_P = b.m_Lambda_nj_P;
332  m_Lambda_nj_coeff = b.m_Lambda_nj_coeff;
333  m_lnActCoeffMolal_Scaled = b.m_lnActCoeffMolal_Scaled;
334  m_lnActCoeffMolal_Unscaled = b.m_lnActCoeffMolal_Unscaled;
335  m_dlnActCoeffMolaldT_Unscaled = b.m_dlnActCoeffMolaldT_Unscaled;
336  m_d2lnActCoeffMolaldT2_Unscaled= b.m_d2lnActCoeffMolaldT2_Unscaled;
337  m_dlnActCoeffMolaldP_Unscaled = b.m_dlnActCoeffMolaldP_Unscaled;
338  m_dlnActCoeffMolaldT_Scaled = b.m_dlnActCoeffMolaldT_Unscaled;
339  m_d2lnActCoeffMolaldT2_Scaled = b.m_d2lnActCoeffMolaldT2_Unscaled;
340  m_dlnActCoeffMolaldP_Scaled = b.m_dlnActCoeffMolaldP_Unscaled;
341 
342  m_gfunc_IJ = b.m_gfunc_IJ;
343  m_g2func_IJ = b.m_g2func_IJ;
344  m_hfunc_IJ = b.m_hfunc_IJ;
345  m_h2func_IJ = b.m_h2func_IJ;
346  m_BMX_IJ = b.m_BMX_IJ;
347  m_BMX_IJ_L = b.m_BMX_IJ_L;
348  m_BMX_IJ_LL = b.m_BMX_IJ_LL;
349  m_BMX_IJ_P = b.m_BMX_IJ_P;
350  m_BprimeMX_IJ = b.m_BprimeMX_IJ;
351  m_BprimeMX_IJ_L = b.m_BprimeMX_IJ_L;
352  m_BprimeMX_IJ_LL = b.m_BprimeMX_IJ_LL;
353  m_BprimeMX_IJ_P = b.m_BprimeMX_IJ_P;
354  m_BphiMX_IJ = b.m_BphiMX_IJ;
355  m_BphiMX_IJ_L = b.m_BphiMX_IJ_L;
356  m_BphiMX_IJ_LL = b.m_BphiMX_IJ_LL;
357  m_BphiMX_IJ_P = b.m_BphiMX_IJ_P;
358  m_Phi_IJ = b.m_Phi_IJ;
359  m_Phi_IJ_L = b.m_Phi_IJ_L;
360  m_Phi_IJ_LL = b.m_Phi_IJ_LL;
361  m_Phi_IJ_P = b.m_Phi_IJ_P;
362  m_Phiprime_IJ = b.m_Phiprime_IJ;
363  m_PhiPhi_IJ = b.m_PhiPhi_IJ;
364  m_PhiPhi_IJ_L = b.m_PhiPhi_IJ_L;
365  m_PhiPhi_IJ_LL = b.m_PhiPhi_IJ_LL;
366  m_PhiPhi_IJ_P = b.m_PhiPhi_IJ_P;
367  m_CMX_IJ = b.m_CMX_IJ;
368  m_CMX_IJ_L = b.m_CMX_IJ_L;
369  m_CMX_IJ_LL = b.m_CMX_IJ_LL;
370  m_CMX_IJ_P = b.m_CMX_IJ_P;
371  m_gamma_tmp = b.m_gamma_tmp;
372 
373  IMS_lnActCoeffMolal_ = b.IMS_lnActCoeffMolal_;
374  IMS_typeCutoff_ = b.IMS_typeCutoff_;
375  IMS_X_o_cutoff_ = b.IMS_X_o_cutoff_;
376  IMS_gamma_o_min_ = b.IMS_gamma_o_min_;
377  IMS_gamma_k_min_ = b.IMS_gamma_k_min_;
378  IMS_cCut_ = b.IMS_cCut_;
379  IMS_slopefCut_ = b.IMS_slopefCut_;
380  IMS_dfCut_ = b.IMS_dfCut_;
381  IMS_efCut_ = b.IMS_efCut_;
382  IMS_afCut_ = b.IMS_afCut_;
383  IMS_bfCut_ = b.IMS_bfCut_;
384  IMS_slopegCut_ = b.IMS_slopegCut_;
385  IMS_dgCut_ = b.IMS_dgCut_;
386  IMS_egCut_ = b.IMS_egCut_;
387  IMS_agCut_ = b.IMS_agCut_;
388  IMS_bgCut_ = b.IMS_bgCut_;
389  MC_X_o_cutoff_ = b.MC_X_o_cutoff_;
390  MC_X_o_min_ = b.MC_X_o_min_;
391  MC_slopepCut_ = b.MC_slopepCut_;
392  MC_dpCut_ = b.MC_dpCut_;
393  MC_epCut_ = b.MC_epCut_;
394  MC_apCut_ = b.MC_apCut_;
395  MC_bpCut_ = b.MC_bpCut_;
396  MC_cpCut_ = b.MC_cpCut_;
397  CROP_ln_gamma_o_min = b.CROP_ln_gamma_o_min;
398  CROP_ln_gamma_o_max = b.CROP_ln_gamma_o_max;
399  CROP_ln_gamma_k_min = b.CROP_ln_gamma_k_min;
400  CROP_ln_gamma_k_max = b.CROP_ln_gamma_k_max;
401  CROP_speciesCropped_ = b.CROP_speciesCropped_;
402  m_CounterIJ = b.m_CounterIJ;
403  m_molalitiesCropped = b.m_molalitiesCropped;
404  m_molalitiesAreCropped= b.m_molalitiesAreCropped;
405  m_debugCalc = b.m_debugCalc;
406  }
407  return *this;
408 }
409 
410 
411 
412 /*
413  *
414  *
415  * test problems:
416  * 1 = NaCl problem - 5 species -
417  * the thermo is read in from an XML file
418  *
419  * speci molality charge
420  * Cl- 6.0954 6.0997E+00 -1
421  * H+ 1.0000E-08 2.1628E-09 1
422  * Na+ 6.0954E+00 6.0997E+00 1
423  * OH- 7.5982E-07 1.3977E-06 -1
424  * HMW_params____beta0MX__beta1MX__beta2MX__CphiMX_____alphaMX__thetaij
425  * 10
426  * 1 2 0.1775 0.2945 0.0 0.00080 2.0 0.0
427  * 1 3 0.0765 0.2664 0.0 0.00127 2.0 0.0
428  * 1 4 0.0 0.0 0.0 0.0 0.0 -0.050
429  * 2 3 0.0 0.0 0.0 0.0 0.0 0.036
430  * 2 4 0.0 0.0 0.0 0.0 0.0 0.0
431  * 3 4 0.0864 0.253 0.0 0.0044 2.0 0.0
432  * Triplet_interaction_parameters_psiaa'_or_psicc'
433  * 2
434  * 1 2 3 -0.004
435  * 1 3 4 -0.006
436  */
437 HMWSoln::HMWSoln(int testProb) :
438  MolalityVPSSTP(),
439  m_formPitzer(PITZERFORM_BASE),
440  m_formPitzerTemp(PITZER_TEMP_CONSTANT),
441  m_formGC(2),
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),
447  m_A_Debye(1.172576), // units = sqrt(kg/gmol)
448  m_waterSS(0),
449  m_densWaterSS(1000.),
450  m_waterProps(0),
451  m_molalitiesAreCropped(false),
452  IMS_typeCutoff_(0),
453  IMS_X_o_cutoff_(0.2),
454  IMS_gamma_o_min_(1.0E-5),
455  IMS_gamma_k_min_(10.0),
456  IMS_cCut_(0.05),
457  IMS_slopefCut_(0.6),
458  IMS_dfCut_(0.0),
459  IMS_efCut_(0.0),
460  IMS_afCut_(0.0),
461  IMS_bfCut_(0.0),
462  IMS_slopegCut_(0.0),
463  IMS_dgCut_(0.0),
464  IMS_egCut_(0.0),
465  IMS_agCut_(0.0),
466  IMS_bgCut_(0.0),
467  MC_X_o_cutoff_(0.0),
468  MC_X_o_min_(0.0),
469  MC_slopepCut_(0.0),
470  MC_dpCut_(0.0),
471  MC_epCut_(0.0),
472  MC_apCut_(0.0),
473  MC_bpCut_(0.0),
474  MC_cpCut_(0.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),
479  m_debugCalc(0)
480 {
481  if (testProb != 1) {
482  printf("unknown test problem\n");
483  exit(EXIT_FAILURE);
484  }
485 
486  constructPhaseFile("HMW_NaCl.xml", "");
487 
488  size_t i = speciesIndex("Cl-");
489  size_t j = speciesIndex("H+");
490  size_t n = i * m_kk + j;
491  size_t ct = m_CounterIJ[n];
492  m_Beta0MX_ij[ct] = 0.1775;
493  m_Beta1MX_ij[ct] = 0.2945;
494  m_CphiMX_ij[ct] = 0.0008;
495  m_Alpha1MX_ij[ct]= 2.000;
496 
497 
498  i = speciesIndex("Cl-");
499  j = speciesIndex("Na+");
500  n = i * m_kk + j;
501  ct = m_CounterIJ[n];
502  m_Beta0MX_ij[ct] = 0.0765;
503  m_Beta1MX_ij[ct] = 0.2664;
504  m_CphiMX_ij[ct] = 0.00127;
505  m_Alpha1MX_ij[ct]= 2.000;
506 
507 
508  i = speciesIndex("Cl-");
509  j = speciesIndex("OH-");
510  n = i * m_kk + j;
511  ct = m_CounterIJ[n];
512  m_Theta_ij[ct] = -0.05;
513 
514  i = speciesIndex("H+");
515  j = speciesIndex("Na+");
516  n = i * m_kk + j;
517  ct = m_CounterIJ[n];
518  m_Theta_ij[ct] = 0.036;
519 
520  i = speciesIndex("Na+");
521  j = speciesIndex("OH-");
522  n = i * m_kk + j;
523  ct = m_CounterIJ[n];
524  m_Beta0MX_ij[ct] = 0.0864;
525  m_Beta1MX_ij[ct] = 0.253;
526  m_CphiMX_ij[ct] = 0.0044;
527  m_Alpha1MX_ij[ct]= 2.000;
528 
529  i = speciesIndex("Cl-");
530  j = speciesIndex("H+");
531  size_t k = speciesIndex("Na+");
532  double param = -0.004;
533  n = i * m_kk *m_kk + j * m_kk + k ;
534  m_Psi_ijk[n] = param;
535  m_Psi_ijk_coeff(0,n) = param;
536  n = i * m_kk *m_kk + k * m_kk + j ;
537  m_Psi_ijk[n] = param;
538  m_Psi_ijk_coeff(0,n) = param;
539  n = j * m_kk *m_kk + i * m_kk + k ;
540  m_Psi_ijk[n] = param;
541  m_Psi_ijk_coeff(0,n) = param;
542  n = j * m_kk *m_kk + k * m_kk + i ;
543  m_Psi_ijk[n] = param;
544  m_Psi_ijk_coeff(0,n) = param;
545  n = k * m_kk *m_kk + j * m_kk + i ;
546  m_Psi_ijk[n] = param;
547  m_Psi_ijk_coeff(0,n) = param;
548  n = k * m_kk *m_kk + i * m_kk + j ;
549  m_Psi_ijk[n] = param;
550  m_Psi_ijk_coeff(0,n) = param;
551 
552  i = speciesIndex("Cl-");
553  j = speciesIndex("Na+");
554  k = speciesIndex("OH-");
555  param = -0.006;
556  n = i * m_kk *m_kk + j * m_kk + k ;
557  m_Psi_ijk[n] = param;
558  m_Psi_ijk_coeff(0,n) = param;
559  n = i * m_kk *m_kk + k * m_kk + j ;
560  m_Psi_ijk[n] = param;
561  m_Psi_ijk_coeff(0,n) = param;
562  n = j * m_kk *m_kk + i * m_kk + k ;
563  m_Psi_ijk[n] = param;
564  m_Psi_ijk_coeff(0,n) = param;
565  n = j * m_kk *m_kk + k * m_kk + i ;
566  m_Psi_ijk[n] = param;
567  m_Psi_ijk_coeff(0,n) = param;
568  n = k * m_kk *m_kk + j * m_kk + i ;
569  m_Psi_ijk[n] = param;
570  m_Psi_ijk_coeff(0,n) = param;
571  n = k * m_kk *m_kk + i * m_kk + j ;
572  m_Psi_ijk[n] = param;
573  m_Psi_ijk_coeff(0,n) = param;
574 
575  printCoeffs();
576 }
577 
578 /*
579  * ~HMWSoln(): (virtual)
580  *
581  * Destructor: does nothing:
582  */
583 HMWSoln::~HMWSoln()
584 {
585  if (m_waterProps) {
586  delete m_waterProps;
587  m_waterProps = 0;
588  }
589 }
590 
591 /*
592  * duplMyselfAsThermoPhase():
593  *
594  * This routine operates at the ThermoPhase level to
595  * duplicate the current object. It uses the copy constructor
596  * defined above.
597  */
598 ThermoPhase* HMWSoln::duplMyselfAsThermoPhase() const
599 {
600  HMWSoln* mtp = new HMWSoln(*this);
601  return (ThermoPhase*) mtp;
602 }
603 
604 /*
605  * Equation of state type flag. The base class returns
606  * zero. Subclasses should define this to return a unique
607  * non-zero value. Constants defined for this purpose are
608  * listed in mix_defs.h.
609  */
610 int HMWSoln::eosType() const
611 {
612  int res;
613  switch (m_formGC) {
614  case 0:
615  res = cHMWSoln0;
616  break;
617  case 1:
618  res = cHMWSoln1;
619  break;
620  case 2:
621  res = cHMWSoln2;
622  break;
623  default:
624  throw CanteraError("eosType", "Unknown type");
625  }
626  return res;
627 }
628 
629 //
630 // -------- Molar Thermodynamic Properties of the Solution ---------------
631 //
632 /*
633  * Molar enthalpy of the solution. Units: J/kmol.
634  */
635 doublereal HMWSoln::enthalpy_mole() const
636 {
637  getPartialMolarEnthalpies(DATA_PTR(m_tmpV));
638  getMoleFractions(DATA_PTR(m_pp));
639  double val = mean_X(DATA_PTR(m_tmpV));
640  return val;
641 }
642 
643 doublereal HMWSoln::relative_enthalpy() const
644 {
645  getPartialMolarEnthalpies(DATA_PTR(m_tmpV));
646  double hbar = mean_X(DATA_PTR(m_tmpV));
647  getEnthalpy_RT(DATA_PTR(m_gamma_tmp));
648  double RT = GasConstant * temperature();
649  for (size_t k = 0; k < m_kk; k++) {
650  m_gamma_tmp[k] *= RT;
651  }
652  double h0bar = mean_X(DATA_PTR(m_gamma_tmp));
653  return (hbar - h0bar);
654 }
655 
656 
657 
658 doublereal HMWSoln::relative_molal_enthalpy() const
659 {
660  double L = relative_enthalpy();
661  getMoleFractions(DATA_PTR(m_tmpV));
662  double xanion = 0.0;
663  size_t kcation = npos;
664  double xcation = 0.0;
665  size_t kanion = npos;
666  const double* charge = DATA_PTR(m_speciesCharge);
667  for (size_t k = 0; k < m_kk; k++) {
668  if (charge[k] > 0.0) {
669  if (m_tmpV[k] > xanion) {
670  xanion = m_tmpV[k];
671  kanion = k;
672  }
673  } else if (charge[k] < 0.0) {
674  if (m_tmpV[k] > xcation) {
675  xcation = m_tmpV[k];
676  kcation = k;
677  }
678  }
679  }
680  if (kcation == npos || kanion == npos) {
681  return L;
682  }
683  double xuse = xcation;
684  double factor = 1;
685  if (xanion < xcation) {
686  xuse = xanion;
687  if (charge[kcation] != 1.0) {
688  factor = charge[kcation];
689  }
690  } else {
691  if (charge[kanion] != 1.0) {
692  factor = charge[kanion];
693  }
694  }
695  xuse = xuse / factor;
696  L = L / xuse;
697  return L;
698 }
699 
700 /*
701  * Molar internal energy of the solution. Units: J/kmol.
702  *
703  * This is calculated from the soln enthalpy and then
704  * subtracting pV.
705  */
706 doublereal HMWSoln::intEnergy_mole() const
707 {
708  double hh = enthalpy_mole();
709  double pres = pressure();
710  double molarV = 1.0/molarDensity();
711  double uu = hh - pres * molarV;
712  return uu;
713 }
714 
715 /*
716  * Molar soln entropy at constant pressure. Units: J/kmol/K.
717  *
718  * This is calculated from the partial molar entropies.
719  */
720 doublereal HMWSoln::entropy_mole() const
721 {
722  getPartialMolarEntropies(DATA_PTR(m_tmpV));
723  return mean_X(DATA_PTR(m_tmpV));
724 }
725 
726 /// Molar Gibbs function. Units: J/kmol.
727 doublereal HMWSoln::gibbs_mole() const
728 {
729  getChemPotentials(DATA_PTR(m_tmpV));
730  return mean_X(DATA_PTR(m_tmpV));
731 }
732 
733 /* Molar heat capacity at constant pressure. Units: J/kmol/K.
734  *
735  * Returns the solution heat capacition at constant pressure.
736  * This is calculated from the partial molar heat capacities.
737  */
738 doublereal HMWSoln::cp_mole() const
739 {
740  getPartialMolarCp(DATA_PTR(m_tmpV));
741  double val = mean_X(DATA_PTR(m_tmpV));
742  return val;
743 }
744 
745 // Molar heat capacity at constant volume. Units: J/kmol/K.
746 doublereal HMWSoln::cv_mole() const
747 {
748  double kappa_t = isothermalCompressibility();
749  double beta = thermalExpansionCoeff();
750  double cp = cp_mole();
751  double tt = temperature();
752  double molarV = molarVolume();
753  double cv = cp - beta * beta * tt * molarV / kappa_t;
754  return cv;
755 }
756 
757 //
758 // ------- Mechanical Equation of State Properties ------------------------
759 //
760 
761 /**
762  * Pressure. Units: Pa.
763  * For this incompressible system, we return the internally stored
764  * independent value of the pressure.
765  */
766 doublereal HMWSoln::pressure() const
767 {
768  return m_Pcurrent;
769 }
770 
771 /*
772  * Set the pressure at constant temperature. Units: Pa.
773  * This method sets a constant within the object.
774  * The mass density is not a function of pressure.
775  */
776 void HMWSoln::setPressure(doublereal p)
777 {
778  setState_TP(temperature(), p);
779 }
780 
781 void HMWSoln::calcDensity()
782 {
783  double* vbar = &m_pp[0];
784  getPartialMolarVolumes(vbar);
785  double* x = &m_tmpV[0];
786  getMoleFractions(x);
787  doublereal vtotal = 0.0;
788  for (size_t i = 0; i < m_kk; i++) {
789  vtotal += vbar[i] * x[i];
790  }
791  doublereal dd = meanMolecularWeight() / vtotal;
792  Phase::setDensity(dd);
793 }
794 
795 /*
796  * The isothermal compressibility. Units: 1/Pa.
797  * The isothermal compressibility is defined as
798  * \f[
799  * \kappa_T = -\frac{1}{v}\left(\frac{\partial v}{\partial P}\right)_T
800  * \f]
801  *
802  * It's equal to zero for this model, since the molar volume
803  * doesn't change with pressure or temperature.
804  */
805 doublereal HMWSoln::isothermalCompressibility() const
806 {
807  throw CanteraError("HMWSoln::isothermalCompressibility",
808  "unimplemented");
809  return 0.0;
810 }
811 
812 /*
813  * The thermal expansion coefficient. Units: 1/K.
814  * The thermal expansion coefficient is defined as
815  *
816  * \f[
817  * \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P
818  * \f]
819  *
820  * It's equal to zero for this model, since the molar volume
821  * doesn't change with pressure or temperature.
822  */
823 doublereal HMWSoln::thermalExpansionCoeff() const
824 {
825  throw CanteraError("HMWSoln::thermalExpansionCoeff",
826  "unimplemented");
827  return 0.0;
828 }
829 
830 double HMWSoln::density() const
831 {
832  // calcDensity();
833  return Phase::density();
834 }
835 
836 /*
837  * Overwritten setDensity() function is necessary because the
838  * density is not an independent variable.
839  *
840  * This function will now throw an error condition
841  *
842  * Note, in general, setting the phase density is now a nonlinear
843  * calculation. P and T are the fundamental variables. This
844  * routine should be revamped to do the nonlinear problem
845  *
846  * @internal May have to adjust the strategy here to make
847  * the eos for these materials slightly compressible, in order
848  * to create a condition where the density is a function of
849  * the pressure.
850  *
851  * This function will now throw an error condition.
852  *
853  * NOTE: This is an overwritten function from the State.h
854  * class
855  */
856 void HMWSoln::setDensity(const doublereal rho)
857 {
858  double dens_old = density();
859 
860  if (rho != dens_old) {
861  throw CanteraError("HMWSoln::setDensity",
862  "Density is not an independent variable");
863  }
864 }
865 
866 /*
867  * Overwritten setMolarDensity() function is necessary because the
868  * density is not an independent variable.
869  *
870  * This function will now throw an error condition.
871  *
872  * NOTE: This is an overwritten function from the State.h
873  * class
874  */
875 void HMWSoln::setMolarDensity(const doublereal rho)
876 {
877  throw CanteraError("HMWSoln::setMolarDensity",
878  "Density is not an independent variable");
879 }
880 
881 /*
882  * Overwritten setTemperature(double) from State.h. This
883  * function sets the temperature, and makes sure that
884  * the value propagates to underlying objects.
885  */
886 void HMWSoln::setTemperature(const doublereal temp)
887 {
888  setState_TP(temp, m_Pcurrent);
889 }
890 
891 /*
892  * Overwritten setTemperature(double) from State.h. This
893  * function sets the temperature, and makes sure that
894  * the value propagates to underlying objects.
895  */
896 void HMWSoln::setState_TP(doublereal temp, doublereal pres)
897 {
898  Phase::setTemperature(temp);
899  /*
900  * Store the current pressure
901  */
902  m_Pcurrent = pres;
903 
904  /*
905  * update the standard state thermo
906  * -> This involves calling the water function and setting the pressure
907  */
908  updateStandardStateThermo();
909  /*
910  * Store the internal density of the water SS.
911  * Note, we would have to do this for all other
912  * species if they had pressure dependent properties.
913  */
914  m_densWaterSS = m_waterSS->density();
915  /*
916  * Calculate all of the other standard volumes
917  * -> note these are constant for now
918  */
919  calcDensity();
920 }
921 
922 //
923 // ------- Activities and Activity Concentrations
924 //
925 
926 /*
927  * This method returns an array of generalized concentrations
928  * \f$ C_k\f$ that are defined such that
929  * \f$ a_k = C_k / C^0_k, \f$ where \f$ C^0_k \f$
930  * is a standard concentration
931  * defined below. These generalized concentrations are used
932  * by kinetics manager classes to compute the forward and
933  * reverse rates of elementary reactions.
934  *
935  * @param c Array of generalized concentrations. The
936  * units depend upon the implementation of the
937  * reaction rate expressions within the phase.
938  */
939 void HMWSoln::getActivityConcentrations(doublereal* c) const
940 {
941  double cs_solvent = standardConcentration();
942  getActivities(c);
943  c[0] *= cs_solvent;
944  if (m_kk > 1) {
945  double cs_solute = standardConcentration(1);
946  for (size_t k = 1; k < m_kk; k++) {
947  c[k] *= cs_solute;
948  }
949  }
950 }
951 
952 /*
953  * The standard concentration \f$ C^0_k \f$ used to normalize
954  * the generalized concentration. In many cases, this quantity
955  * will be the same for all species in a phase - for example,
956  * for an ideal gas \f$ C^0_k = P/\hat R T \f$. For this
957  * reason, this method returns a single value, instead of an
958  * array. However, for phases in which the standard
959  * concentration is species-specific (e.g. surface species of
960  * different sizes), this method may be called with an
961  * optional parameter indicating the species.
962  *
963  * For the time being we will use the concentration of pure
964  * solvent for the the standard concentration of the solvent.
965  * We will use the concentration of the pure solvent
966  * multipled by Mnaught (kg solvent / gmol solvent) for
967  * the standard concentration of all solute species.
968  * This has the effect of making reaction rates
969  * based on the molality of species proportional to the
970  * molality of the species, but have units based on assuming
971  * all species concentrations have units of kmol/m3.
972  *
973  */
974 doublereal HMWSoln::standardConcentration(size_t k) const
975 {
976  getStandardVolumes(DATA_PTR(m_tmpV));
977  double mvSolvent = m_tmpV[m_indexSolvent];
978  if (k > 0) {
979  return m_Mnaught / mvSolvent;
980  }
981  return 1.0 / mvSolvent;
982 }
983 
984 /*
985  * Returns the natural logarithm of the standard
986  * concentration of the kth species
987  */
988 doublereal HMWSoln::logStandardConc(size_t k) const
989 {
990  double c_solvent = standardConcentration(k);
991  return log(c_solvent);
992 }
993 
994 /*
995  * Returns the units of the standard and general concentrations
996  * Note they have the same units, as their divisor is
997  * defined to be equal to the activity of the kth species
998  * in the solution, which is unitless.
999  *
1000  * This routine is used in print out applications where the
1001  * units are needed. Usually, MKS units are assumed throughout
1002  * the program and in the XML input files.
1003  *
1004  * On return uA contains the powers of the units (MKS assumed)
1005  * of the standard concentrations and generalized concentrations
1006  * for the kth species.
1007  *
1008  * uA[0] = kmol units - default = 1
1009  * uA[1] = m units - default = -nDim(), the number of spatial
1010  * dimensions in the Phase class.
1011  * uA[2] = kg units - default = 0;
1012  * uA[3] = Pa(pressure) units - default = 0;
1013  * uA[4] = Temperature units - default = 0;
1014  * uA[5] = time units - default = 0
1015  */
1016 void HMWSoln::getUnitsStandardConc(double* uA, int k, int sizeUA) const
1017 {
1018  for (int i = 0; i < sizeUA; i++) {
1019  if (i == 0) {
1020  uA[0] = 1.0;
1021  }
1022  if (i == 1) {
1023  uA[1] = -int(nDim());
1024  }
1025  if (i == 2) {
1026  uA[2] = 0.0;
1027  }
1028  if (i == 3) {
1029  uA[3] = 0.0;
1030  }
1031  if (i == 4) {
1032  uA[4] = 0.0;
1033  }
1034  if (i == 5) {
1035  uA[5] = 0.0;
1036  }
1037  }
1038 }
1039 
1040 /*
1041  * Get the array of non-dimensional activities at
1042  * the current solution temperature, pressure, and
1043  * solution concentration.
1044  * (note solvent activity coefficient is on the molar scale).
1045  *
1046  */
1047 void HMWSoln::getActivities(doublereal* ac) const
1048 {
1049  updateStandardStateThermo();
1050  /*
1051  * Update the molality array, m_molalities()
1052  * This requires an update due to mole fractions
1053  */
1054  s_update_lnMolalityActCoeff();
1055  /*
1056  * Now calculate the array of activities.
1057  */
1058  for (size_t k = 0; k < m_kk; k++) {
1059  if (k != m_indexSolvent) {
1060  ac[k] = m_molalities[k] * exp(m_lnActCoeffMolal_Scaled[k]);
1061  }
1062  }
1063  double xmolSolvent = moleFraction(m_indexSolvent);
1064  ac[m_indexSolvent] =
1065  exp(m_lnActCoeffMolal_Scaled[m_indexSolvent]) * xmolSolvent;
1066  /*
1067  * Apply the pH scale
1068  */
1069  //applyphScale(ac);
1070 }
1071 
1072 /*
1073  * getUnscaledMolalityActivityCoefficients() (virtual, const)
1074  *
1075  * Get the array of non-dimensional Molality based
1076  * activity coefficients at
1077  * the current solution temperature, pressure, and
1078  * solution concentration.
1079  * (note solvent activity coefficient is on the molar scale).
1080  *
1081  * Note, most of the work is done in an internal private routine
1082  */
1083 void HMWSoln::
1084 getUnscaledMolalityActivityCoefficients(doublereal* acMolality) const
1085 {
1086  updateStandardStateThermo();
1087  A_Debye_TP(-1.0, -1.0);
1088  s_update_lnMolalityActCoeff();
1089  std::copy(m_lnActCoeffMolal_Unscaled.begin(), m_lnActCoeffMolal_Unscaled.end(), acMolality);
1090  for (size_t k = 0; k < m_kk; k++) {
1091  acMolality[k] = exp(acMolality[k]);
1092  }
1093 }
1094 
1095 //
1096 // ------ Partial Molar Properties of the Solution -----------------
1097 //
1098 /*
1099  * Get the species chemical potentials. Units: J/kmol.
1100  *
1101  * This function returns a vector of chemical potentials of the
1102  * species in solution.
1103  *
1104  * \f[
1105  * \mu_k = \mu^{o}_k(T,P) + R T ln(m_k)
1106  * \f]
1107  *
1108  * \f[
1109  * \mu_solvent = \mu^{o}_solvent(T,P) +
1110  * R T ((X_solvent - 1.0) / X_solvent)
1111  * \f]
1112  */
1113 void HMWSoln::getChemPotentials(doublereal* mu) const
1114 {
1115  double xx;
1116  const double xxSmall = 1.0E-150;
1117  /*
1118  * First get the standard chemical potentials in
1119  * molar form.
1120  * -> this requires updates of standard state as a function
1121  * of T and P
1122  */
1123  getStandardChemPotentials(mu);
1124  /*
1125  * Update the activity coefficients
1126  * This also updates the internal molality array.
1127  */
1128  s_update_lnMolalityActCoeff();
1129  doublereal RT = GasConstant * temperature();
1130  double xmolSolvent = moleFraction(m_indexSolvent);
1131  for (size_t k = 0; k < m_kk; k++) {
1132  if (m_indexSolvent != k) {
1133  xx = std::max(m_molalities[k], xxSmall);
1134  mu[k] += RT * (log(xx) + m_lnActCoeffMolal_Scaled[k]);
1135  }
1136  }
1137  xx = std::max(xmolSolvent, xxSmall);
1138  mu[m_indexSolvent] +=
1139  RT * (log(xx) + m_lnActCoeffMolal_Scaled[m_indexSolvent]);
1140 }
1141 
1142 
1143 /*
1144  * Returns an array of partial molar enthalpies for the species
1145  * in the mixture.
1146  * Units (J/kmol)
1147  *
1148  * For this phase, the partial molar enthalpies are equal to the
1149  * standard state enthalpies modified by the derivative of the
1150  * molality-based activity coefficient wrt temperature
1151  *
1152  * \f[
1153  * \bar h_k(T,P) = h^{\triangle}_k(T,P) - R T^2 \frac{d \ln(\gamma_k^\triangle)}{dT}
1154  * \f]
1155  * The solvent partial molar enthalpy is equal to
1156  * \f[
1157  * \bar h_o(T,P) = h^{o}_o(T,P) - R T^2 \frac{d \ln(a_o)}{dT}
1158  * \f]
1159  *
1160  *
1161  */
1162 void HMWSoln::getPartialMolarEnthalpies(doublereal* hbar) const
1163 {
1164  /*
1165  * Get the nondimensional standard state enthalpies
1166  */
1167  getEnthalpy_RT(hbar);
1168  /*
1169  * dimensionalize it.
1170  */
1171  double T = temperature();
1172  double RT = GasConstant * T;
1173  for (size_t k = 0; k < m_kk; k++) {
1174  hbar[k] *= RT;
1175  }
1176  /*
1177  * Update the activity coefficients, This also update the
1178  * internally stored molalities.
1179  */
1180  s_update_lnMolalityActCoeff();
1181  s_update_dlnMolalityActCoeff_dT();
1182  double RTT = RT * T;
1183  for (size_t k = 0; k < m_kk; k++) {
1184  hbar[k] -= RTT * m_dlnActCoeffMolaldT_Scaled[k];
1185  }
1186 }
1187 
1188 /*
1189  * getPartialMolarEntropies() (virtual, const)
1190  *
1191  * Returns an array of partial molar entropies of the species in the
1192  * solution. Units: J/kmol.
1193  *
1194  * Maxwell's equations provide an insight in how to calculate this
1195  * (p.215 Smith and Van Ness)
1196  *
1197  * d(chemPot_i)/dT = -sbar_i
1198  *
1199  * Combining this with the expression H = G + TS yields:
1200  *
1201  * \f[
1202  * \bar s_k(T,P) = s^{\triangle}_k(T,P)
1203  * - R \ln( \gamma^{\triangle}_k \frac{m_k}{m^{\triangle}}))
1204  * - R T \frac{d \ln(\gamma^{\triangle}_k) }{dT}
1205  * \f]
1206  * \f[
1207  * \bar s_o(T,P) = s^o_o(T,P) - R \ln(a_o)
1208  * - R T \frac{d \ln(a_o)}{dT}
1209  * \f]
1210  *
1211  * The reference-state pure-species entropies,\f$ \hat s^0_k(T) \f$,
1212  * at the reference pressure, \f$ P_{ref} \f$, are computed by the
1213  * species thermodynamic
1214  * property manager. They are polynomial functions of temperature.
1215  * @see SpeciesThermo
1216  */
1217 void HMWSoln::
1218 getPartialMolarEntropies(doublereal* sbar) const
1219 {
1220  /*
1221  * Get the standard state entropies at the temperature
1222  * and pressure of the solution.
1223  */
1224  getEntropy_R(sbar);
1225  /*
1226  * Dimensionalize the entropies
1227  */
1228  doublereal R = GasConstant;
1229  for (size_t k = 0; k < m_kk; k++) {
1230  sbar[k] *= R;
1231  }
1232  /*
1233  * Update the activity coefficients, This also update the
1234  * internally stored molalities.
1235  */
1236  s_update_lnMolalityActCoeff();
1237  /*
1238  * First we will add in the obvious dependence on the T
1239  * term out front of the log activity term
1240  */
1241  doublereal mm;
1242  for (size_t k = 0; k < m_kk; k++) {
1243  if (k != m_indexSolvent) {
1244  mm = std::max(SmallNumber, m_molalities[k]);
1245  sbar[k] -= R * (log(mm) + m_lnActCoeffMolal_Scaled[k]);
1246  }
1247  }
1248  double xmolSolvent = moleFraction(m_indexSolvent);
1249  mm = std::max(SmallNumber, xmolSolvent);
1250  sbar[m_indexSolvent] -= R *(log(mm) + m_lnActCoeffMolal_Scaled[m_indexSolvent]);
1251  /*
1252  * Check to see whether activity coefficients are temperature
1253  * dependent. If they are, then calculate the their temperature
1254  * derivatives and add them into the result.
1255  */
1256  s_update_dlnMolalityActCoeff_dT();
1257  double RT = R * temperature();
1258  for (size_t k = 0; k < m_kk; k++) {
1259  sbar[k] -= RT * m_dlnActCoeffMolaldT_Scaled[k];
1260  }
1261 }
1262 
1263 /*
1264  * getPartialMolarVolumes() (virtual, const)
1265  *
1266  * Returns an array of partial molar volumes of the species
1267  * in the solution. Units: m^3 kmol-1.
1268  *
1269  * For this solution, the partial molar volumes are a
1270  * complex function of pressure.
1271  *
1272  * The general relation is
1273  *
1274  * vbar_i = d(chemPot_i)/dP at const T, n
1275  *
1276  * = V0_i + d(Gex)/dP)_T,M
1277  *
1278  * = V0_i + RT d(lnActCoeffi)dP _T,M
1279  *
1280  */
1281 void HMWSoln::getPartialMolarVolumes(doublereal* vbar) const
1282 {
1283  /*
1284  * Get the standard state values in m^3 kmol-1
1285  */
1286  getStandardVolumes(vbar);
1287  /*
1288  * Update the derivatives wrt the activity coefficients.
1289  */
1290  s_update_lnMolalityActCoeff();
1291  s_update_dlnMolalityActCoeff_dP();
1292  double T = temperature();
1293  double RT = GasConstant * T;
1294  for (size_t k = 0; k < m_kk; k++) {
1295  vbar[k] += RT * m_dlnActCoeffMolaldP_Scaled[k];
1296  }
1297 }
1298 
1299 /*
1300  * Partial molar heat capacity of the solution:
1301  * The kth partial molar heat capacity is equal to
1302  * the temperature derivative of the partial molar
1303  * enthalpy of the kth species in the solution at constant
1304  * P and composition (p. 220 Smith and Van Ness).
1305  *
1306  * \f[
1307  * \bar C_{p,k}(T,P) = C^{\triangle}_{p,k}(T,P)
1308  * - 2 R T \frac{d \ln( \gamma^{\triangle}_k)}{dT}
1309  * - R T^2 \frac{d^2 \ln(\gamma^{\triangle}_k) }{{dT}^2}
1310  * \f]
1311  * \f[
1312  * \bar C_{p,o}(T,P) = C^o_{p,o}(T,P)
1313  * - 2 R T \frac{d \ln(a_o)}{dT}
1314  * - R T^2 \frac{d^2 \ln(a_o)}{{dT}^2}
1315  * \f]
1316  */
1317 void HMWSoln::getPartialMolarCp(doublereal* cpbar) const
1318 {
1319  /*
1320  * Get the nondimensional gibbs standard state of the
1321  * species at the T and P of the solution.
1322  */
1323  getCp_R(cpbar);
1324 
1325  for (size_t k = 0; k < m_kk; k++) {
1326  cpbar[k] *= GasConstant;
1327  }
1328  /*
1329  * Update the activity coefficients, This also update the
1330  * internally stored molalities.
1331  */
1332  s_update_lnMolalityActCoeff();
1333  s_update_dlnMolalityActCoeff_dT();
1334  s_update_d2lnMolalityActCoeff_dT2();
1335  double T = temperature();
1336  double RT = GasConstant * T;
1337  double RTT = RT * T;
1338  for (size_t k = 0; k < m_kk; k++) {
1339  cpbar[k] -= (2.0 * RT * m_dlnActCoeffMolaldT_Scaled[k] +
1340  RTT * m_d2lnActCoeffMolaldT2_Scaled[k]);
1341  }
1342 }
1343 
1344 /*
1345  * Updates the standard state thermodynamic functions at the current T and
1346  * P of the solution.
1347  *
1348  * @internal
1349  *
1350  * This function gets called for every call to functions in this
1351  * class. It checks to see whether the temperature or pressure has changed and
1352  * thus the ss thermodynamics functions for all of the species
1353  * must be recalculated.
1354  */
1355 // void HMWSoln::_updateStandardStateThermo() const {
1356 //doublereal tnow = temperature();
1357 // doublereal pnow = m_Pcurrent;
1358 // if (m_waterSS) {
1359 // m_waterSS->setTempPressure(tnow, pnow);
1360 // }
1361 // m_VPSS_ptr->setState_TP(tnow, pnow);
1362 // VPStandardStateTP::updateStandardStateThermo();
1363 //}
1364 
1365 /*
1366  * ------ Thermodynamic Values for the Species Reference States ---
1367  */
1368 
1369 // -> This is handled by VPStandardStatesTP
1370 
1371 /*
1372  * -------------- Utilities -------------------------------
1373  */
1374 
1375 /*
1376  * @internal
1377  * Set equation of state parameters. The number and meaning of
1378  * these depends on the subclass.
1379  * @param n number of parameters
1380  * @param c array of <I>n</I> coefficients
1381  *
1382  */
1383 void HMWSoln::setParameters(int n, doublereal* const c)
1384 {
1385 }
1386 
1387 void HMWSoln::getParameters(int& n, doublereal* const c) const
1388 {
1389 }
1390 /*
1391  * Set equation of state parameter values from XML
1392  * entries. This method is called by function importPhase in
1393  * file importCTML.cpp when processing a phase definition in
1394  * an input file. It should be overloaded in subclasses to set
1395  * any parameters that are specific to that particular phase
1396  * model.
1397  *
1398  * @param eosdata An XML_Node object corresponding to
1399  * the "thermo" entry for this phase in the input file.
1400  *
1401  * HKM -> Right now, the parameters are set elsewhere (initThermoXML)
1402  * It just didn't seem to fit.
1403  */
1404 void HMWSoln::setParametersFromXML(const XML_Node& eosdata)
1405 {
1406 }
1407 
1408 /*
1409  * Get the saturation pressure for a given temperature.
1410  * Note the limitations of this function. Stability considerations
1411  * concerning multiphase equilibrium are ignored in this
1412  * calculation. Therefore, the call is made directly to the SS of
1413  * water underneath. The object is put back into its original
1414  * state at the end of the call.
1415  */
1416 doublereal HMWSoln::satPressure(doublereal t) const
1417 {
1418  double p_old = pressure();
1419  double t_old = temperature();
1420  double pres = m_waterSS->satPressure(t);
1421  /*
1422  * Set the underlying object back to its original state.
1423  */
1424  m_waterSS->setState_TP(t_old, p_old);
1425  return pres;
1426 }
1427 
1428 /*
1429  * Report the molar volume of species k
1430  *
1431  * units - \f$ m^3 kmol^-1 \f$
1432  */
1433 double HMWSoln::speciesMolarVolume(int k) const
1434 {
1435  double vol = m_speciesSize[k];
1436  if (k == 0) {
1437  double dd = m_waterSS->density();
1438  vol = molecularWeight(0)/dd;
1439  }
1440  return vol;
1441 }
1442 
1443 /*
1444  * A_Debye_TP() (virtual)
1445  *
1446  * Returns the A_Debye parameter as a function of temperature
1447  * and pressure. This function also sets the internal value
1448  * of the parameter within the object, if it is changeable.
1449  *
1450  * The default is to assume that it is constant, given
1451  * in the initialization process and stored in the
1452  * member double, m_A_Debye
1453  *
1454  * A_Debye = (1/(8 Pi)) sqrt(2 Na dw /1000)
1455  * (e e/(epsilon R T))^3/2
1456  *
1457  * where epsilon = e_rel * e_naught
1458  *
1459  * Note, this is si units. Frequently, gaussian units are
1460  * used in Pitzer's papers where D is used, D = epsilon/(4 Pi)
1461  * units = A_Debye has units of sqrt(gmol kg-1).
1462  */
1463 double HMWSoln::A_Debye_TP(double tempArg, double presArg) const
1464 {
1465  double T = temperature();
1466  double A;
1467  if (tempArg != -1.0) {
1468  T = tempArg;
1469  }
1470  double P = pressure();
1471  if (presArg != -1.0) {
1472  P = presArg;
1473  }
1474 
1475  switch (m_form_A_Debye) {
1476  case A_DEBYE_CONST:
1477  A = m_A_Debye;
1478  break;
1479  case A_DEBYE_WATER:
1480  A = m_waterProps->ADebye(T, P, 0);
1481  m_A_Debye = A;
1482  break;
1483  default:
1484  printf("shouldn't be here\n");
1485  exit(EXIT_FAILURE);
1486  }
1487  return A;
1488 }
1489 
1490 /*
1491  * dA_DebyedT_TP() (virtual)
1492  *
1493  * Returns the derivative of the A_Debye parameter with
1494  * respect to temperature as a function of temperature
1495  * and pressure.
1496  *
1497  * units = A_Debye has units of sqrt(gmol kg-1).
1498  * Temp has units of Kelvin.
1499  */
1500 double HMWSoln::dA_DebyedT_TP(double tempArg, double presArg) const
1501 {
1502  doublereal T = temperature();
1503  if (tempArg != -1.0) {
1504  T = tempArg;
1505  }
1506  doublereal P = pressure();
1507  if (presArg != -1.0) {
1508  P = presArg;
1509  }
1510  doublereal dAdT;
1511  switch (m_form_A_Debye) {
1512  case A_DEBYE_CONST:
1513  dAdT = 0.0;
1514  break;
1515  case A_DEBYE_WATER:
1516  dAdT = m_waterProps->ADebye(T, P, 1);
1517  //dAdT = WaterProps::ADebye(T, P, 1);
1518  break;
1519  default:
1520  printf("shouldn't be here\n");
1521  exit(EXIT_FAILURE);
1522  }
1523  return dAdT;
1524 }
1525 
1526 /*
1527  * dA_DebyedP_TP() (virtual)
1528  *
1529  * Returns the derivative of the A_Debye parameter with
1530  * respect to pressure, as a function of temperature
1531  * and pressure.
1532  *
1533  * units = A_Debye has units of sqrt(gmol kg-1).
1534  * Pressure has units of pascals.
1535  */
1536 double HMWSoln::dA_DebyedP_TP(double tempArg, double presArg) const
1537 {
1538  double T = temperature();
1539  if (tempArg != -1.0) {
1540  T = tempArg;
1541  }
1542  double P = pressure();
1543  if (presArg != -1.0) {
1544  P = presArg;
1545  }
1546  double dAdP;
1547  switch (m_form_A_Debye) {
1548  case A_DEBYE_CONST:
1549  dAdP = 0.0;
1550  break;
1551  case A_DEBYE_WATER:
1552  dAdP = m_waterProps->ADebye(T, P, 3);
1553  break;
1554  default:
1555  printf("shouldn't be here\n");
1556  exit(EXIT_FAILURE);
1557  }
1558  return dAdP;
1559 }
1560 
1561 
1562 /*
1563  * Calculate the DH Parameter used for the Enthalpy calcalations
1564  *
1565  * ADebye_L = 4 R T**2 d(Aphi) / dT
1566  *
1567  * where Aphi = A_Debye/3
1568  *
1569  * units -> J / (kmolK) * sqrt( kg/gmol)
1570  *
1571  */
1572 double HMWSoln::ADebye_L(double tempArg, double presArg) const
1573 {
1574  double dAdT = dA_DebyedT_TP();
1575  double dAphidT = dAdT /3.0;
1576  double T = temperature();
1577  if (tempArg != -1.0) {
1578  T = tempArg;
1579  }
1580  double retn = dAphidT * (4.0 * GasConstant * T * T);
1581  return retn;
1582 }
1583 
1584 /*
1585  * Calculate the DH Parameter used for the Volume calcalations
1586  *
1587  * ADebye_V = - 4 R T d(Aphi) / dP
1588  *
1589  * where Aphi = A_Debye/3
1590  *
1591  * units -> J / (kmolK) * sqrt( kg/gmol)
1592  *
1593  */
1594 double HMWSoln::ADebye_V(double tempArg, double presArg) const
1595 {
1596  double dAdP = dA_DebyedP_TP();
1597  double dAphidP = dAdP /3.0;
1598  double T = temperature();
1599  if (tempArg != -1.0) {
1600  T = tempArg;
1601  }
1602  double retn = - dAphidP * (4.0 * GasConstant * T);
1603  return retn;
1604 }
1605 
1606 /*
1607  * Return Pitzer's definition of A_J. This is basically the
1608  * temperature derivative of A_L, and the second derivative
1609  * of Aphi
1610  * It's the DH parameter used in heat capacity calculations
1611  *
1612  * A_J = 2 A_L/T + 4 * R * T * T * d2(A_phi)/dT2
1613  *
1614  * Units = sqrt(kg/gmol) (R)
1615  *
1616  * where
1617  * ADebye_L = 4 R T**2 d(Aphi) / dT
1618  *
1619  * where Aphi = A_Debye/3
1620  *
1621  * units -> J / (kmolK) * sqrt( kg/gmol)
1622  *
1623  */
1624 double HMWSoln::ADebye_J(double tempArg, double presArg) const
1625 {
1626  double T = temperature();
1627  if (tempArg != -1.0) {
1628  T = tempArg;
1629  }
1630  double A_L = ADebye_L(T, presArg);
1631  double d2 = d2A_DebyedT2_TP(T, presArg);
1632  double d2Aphi = d2 / 3.0;
1633  double retn = 2.0 * A_L / T + 4.0 * GasConstant * T * T *d2Aphi;
1634  return retn;
1635 }
1636 
1637 /*
1638  * d2A_DebyedT2_TP() (virtual)
1639  *
1640  * Returns the 2nd derivative of the A_Debye parameter with
1641  * respect to temperature as a function of temperature
1642  * and pressure.
1643  *
1644  * units = A_Debye has units of sqrt(gmol kg-1).
1645  * Temp has units of Kelvin.
1646  */
1647 double HMWSoln::d2A_DebyedT2_TP(double tempArg, double presArg) const
1648 {
1649  double T = temperature();
1650  if (tempArg != -1.0) {
1651  T = tempArg;
1652  }
1653  double P = pressure();
1654  if (presArg != -1.0) {
1655  P = presArg;
1656  }
1657  double d2AdT2;
1658  switch (m_form_A_Debye) {
1659  case A_DEBYE_CONST:
1660  d2AdT2 = 0.0;
1661  break;
1662  case A_DEBYE_WATER:
1663  d2AdT2 = m_waterProps->ADebye(T, P, 2);
1664  break;
1665  default:
1666  printf("shouldn't be here\n");
1667  exit(EXIT_FAILURE);
1668  }
1669  return d2AdT2;
1670 }
1671 
1672 /*
1673  * ----------- Critical State Properties --------------------------
1674  */
1675 
1676 /*
1677  * ---------- Other Property Functions
1678  */
1679 double HMWSoln::AionicRadius(int k) const
1680 {
1681  return m_Aionic[k];
1682 }
1683 
1684 /*
1685  * ------------ Private and Restricted Functions ------------------
1686  */
1687 
1688 /**
1689  * Bail out of functions with an error exit if they are not
1690  * implemented.
1691  */
1692 doublereal HMWSoln::err(std::string msg) const
1693 {
1694  throw CanteraError("HMWSoln",
1695  "Unfinished func called: " + msg);
1696  return 0.0;
1697 }
1698 
1699 
1700 
1701 /*
1702  * initLengths():
1703  *
1704  * This internal function adjusts the lengths of arrays based on
1705  * the number of species. This is done before these arrays are
1706  * populated with parameter values.
1707  */
1708 void HMWSoln::initLengths()
1709 {
1710  m_kk = nSpecies();
1711 
1712  /*
1713  * Resize lengths equal to the number of species in
1714  * the phase.
1715  */
1716  m_electrolyteSpeciesType.resize(m_kk, cEST_polarNeutral);
1717  m_speciesSize.resize(m_kk);
1718  m_speciesCharge_Stoich.resize(m_kk, 0.0);
1719  m_Aionic.resize(m_kk, 0.0);
1720 
1721  m_expg0_RT.resize(m_kk, 0.0);
1722  m_pe.resize(m_kk, 0.0);
1723  m_pp.resize(m_kk, 0.0);
1724  m_tmpV.resize(m_kk, 0.0);
1725  m_molalitiesCropped.resize(m_kk, 0.0);
1726 
1727  size_t maxCounterIJlen = 1 + (m_kk-1) * (m_kk-2) / 2;
1728 
1729  /*
1730  * Figure out the size of the temperature coefficient
1731  * arrays
1732  */
1733  int TCoeffLength = 1;
1734  if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
1735  TCoeffLength = 2;
1736  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
1737  TCoeffLength = 5;
1738  }
1739 
1740  m_Beta0MX_ij.resize(maxCounterIJlen, 0.0);
1741  m_Beta0MX_ij_L.resize(maxCounterIJlen, 0.0);
1742  m_Beta0MX_ij_LL.resize(maxCounterIJlen, 0.0);
1743  m_Beta0MX_ij_P.resize(maxCounterIJlen, 0.0);
1744  m_Beta0MX_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
1745 
1746  m_Beta1MX_ij.resize(maxCounterIJlen, 0.0);
1747  m_Beta1MX_ij_L.resize(maxCounterIJlen, 0.0);
1748  m_Beta1MX_ij_LL.resize(maxCounterIJlen, 0.0);
1749  m_Beta1MX_ij_P.resize(maxCounterIJlen, 0.0);
1750  m_Beta1MX_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
1751 
1752  m_Beta2MX_ij.resize(maxCounterIJlen, 0.0);
1753  m_Beta2MX_ij_L.resize(maxCounterIJlen, 0.0);
1754  m_Beta2MX_ij_LL.resize(maxCounterIJlen, 0.0);
1755  m_Beta2MX_ij_P.resize(maxCounterIJlen, 0.0);
1756  m_Beta2MX_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
1757 
1758  m_CphiMX_ij.resize(maxCounterIJlen, 0.0);
1759  m_CphiMX_ij_L.resize(maxCounterIJlen, 0.0);
1760  m_CphiMX_ij_LL.resize(maxCounterIJlen, 0.0);
1761  m_CphiMX_ij_P.resize(maxCounterIJlen, 0.0);
1762  m_CphiMX_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
1763 
1764  m_Alpha1MX_ij.resize(maxCounterIJlen, 2.0);
1765  m_Alpha2MX_ij.resize(maxCounterIJlen, 12.0);
1766  m_Theta_ij.resize(maxCounterIJlen, 0.0);
1767  m_Theta_ij_L.resize(maxCounterIJlen, 0.0);
1768  m_Theta_ij_LL.resize(maxCounterIJlen, 0.0);
1769  m_Theta_ij_P.resize(maxCounterIJlen, 0.0);
1770  m_Theta_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
1771 
1772  size_t n = m_kk*m_kk*m_kk;
1773  m_Psi_ijk.resize(n, 0.0);
1774  m_Psi_ijk_L.resize(n, 0.0);
1775  m_Psi_ijk_LL.resize(n, 0.0);
1776  m_Psi_ijk_P.resize(n, 0.0);
1777  m_Psi_ijk_coeff.resize(TCoeffLength, n, 0.0);
1778 
1779  m_Lambda_nj.resize(m_kk, m_kk, 0.0);
1780  m_Lambda_nj_L.resize(m_kk, m_kk, 0.0);
1781  m_Lambda_nj_LL.resize(m_kk, m_kk, 0.0);
1782  m_Lambda_nj_P.resize(m_kk, m_kk, 0.0);
1783  m_Lambda_nj_coeff.resize(TCoeffLength, m_kk * m_kk, 0.0);
1784 
1785  m_Mu_nnn.resize(m_kk, 0.0);
1786  m_Mu_nnn_L.resize(m_kk, 0.0);
1787  m_Mu_nnn_LL.resize(m_kk, 0.0);
1788  m_Mu_nnn_P.resize(m_kk, 0.0);
1789  m_Mu_nnn_coeff.resize(TCoeffLength, m_kk, 0.0);
1790 
1791  m_lnActCoeffMolal_Scaled.resize(m_kk, 0.0);
1792  m_dlnActCoeffMolaldT_Scaled.resize(m_kk, 0.0);
1793  m_d2lnActCoeffMolaldT2_Scaled.resize(m_kk, 0.0);
1794  m_dlnActCoeffMolaldP_Scaled.resize(m_kk, 0.0);
1795 
1796  m_lnActCoeffMolal_Unscaled.resize(m_kk, 0.0);
1797  m_dlnActCoeffMolaldT_Unscaled.resize(m_kk, 0.0);
1798  m_d2lnActCoeffMolaldT2_Unscaled.resize(m_kk, 0.0);
1799  m_dlnActCoeffMolaldP_Unscaled.resize(m_kk, 0.0);
1800 
1801  m_CounterIJ.resize(m_kk*m_kk, 0);
1802 
1803  m_gfunc_IJ.resize(maxCounterIJlen, 0.0);
1804  m_g2func_IJ.resize(maxCounterIJlen, 0.0);
1805  m_hfunc_IJ.resize(maxCounterIJlen, 0.0);
1806  m_h2func_IJ.resize(maxCounterIJlen, 0.0);
1807  m_BMX_IJ.resize(maxCounterIJlen, 0.0);
1808  m_BMX_IJ_L.resize(maxCounterIJlen, 0.0);
1809  m_BMX_IJ_LL.resize(maxCounterIJlen, 0.0);
1810  m_BMX_IJ_P.resize(maxCounterIJlen, 0.0);
1811  m_BprimeMX_IJ.resize(maxCounterIJlen, 0.0);
1812  m_BprimeMX_IJ_L.resize(maxCounterIJlen, 0.0);
1813  m_BprimeMX_IJ_LL.resize(maxCounterIJlen, 0.0);
1814  m_BprimeMX_IJ_P.resize(maxCounterIJlen, 0.0);
1815  m_BphiMX_IJ.resize(maxCounterIJlen, 0.0);
1816  m_BphiMX_IJ_L.resize(maxCounterIJlen, 0.0);
1817  m_BphiMX_IJ_LL.resize(maxCounterIJlen, 0.0);
1818  m_BphiMX_IJ_P.resize(maxCounterIJlen, 0.0);
1819  m_Phi_IJ.resize(maxCounterIJlen, 0.0);
1820  m_Phi_IJ_L.resize(maxCounterIJlen, 0.0);
1821  m_Phi_IJ_LL.resize(maxCounterIJlen, 0.0);
1822  m_Phi_IJ_P.resize(maxCounterIJlen, 0.0);
1823  m_Phiprime_IJ.resize(maxCounterIJlen, 0.0);
1824  m_PhiPhi_IJ.resize(maxCounterIJlen, 0.0);
1825  m_PhiPhi_IJ_L.resize(maxCounterIJlen, 0.0);
1826  m_PhiPhi_IJ_LL.resize(maxCounterIJlen, 0.0);
1827  m_PhiPhi_IJ_P.resize(maxCounterIJlen, 0.0);
1828  m_CMX_IJ.resize(maxCounterIJlen, 0.0);
1829  m_CMX_IJ_L.resize(maxCounterIJlen, 0.0);
1830  m_CMX_IJ_LL.resize(maxCounterIJlen, 0.0);
1831  m_CMX_IJ_P.resize(maxCounterIJlen, 0.0);
1832 
1833  m_gamma_tmp.resize(m_kk, 0.0);
1834 
1835  IMS_lnActCoeffMolal_.resize(m_kk, 0.0);
1836  CROP_speciesCropped_.resize(m_kk, 0);
1837 
1838  counterIJ_setup();
1839 }
1840 
1841 /**
1842  * Calculate the natural log of the molality-based
1843  * activity coefficients.
1844  *
1845  */
1846 void HMWSoln::s_update_lnMolalityActCoeff() const
1847 {
1848 
1849  /*
1850  * Calculate the molalities. Currently, the molalities
1851  * may not be current with respect to the contents of the
1852  * State objects' data.
1853  */
1854  calcMolalities();
1855  /*
1856  * Calculate a cropped set of molalities that will be used
1857  * in all activity coefficient calculations.
1858  */
1859  calcMolalitiesCropped();
1860  /*
1861  * Calculate the stoichiometric ionic charge. This isn't used in the
1862  * Pitzer formulation.
1863  */
1864  m_IionicMolalityStoich = 0.0;
1865  for (size_t k = 0; k < m_kk; k++) {
1866  double z_k = m_speciesCharge[k];
1867  double zs_k1 = m_speciesCharge_Stoich[k];
1868  if (z_k == zs_k1) {
1869  m_IionicMolalityStoich += m_molalities[k] * z_k * z_k;
1870  } else {
1871  double zs_k2 = z_k - zs_k1;
1872  m_IionicMolalityStoich
1873  += m_molalities[k] * (zs_k1 * zs_k1 + zs_k2 * zs_k2);
1874  }
1875  }
1876 
1877 
1878  /*
1879  * Update the temperature dependence of the pitzer coefficients
1880  * and their derivatives
1881  */
1882  s_updatePitzer_CoeffWRTemp();
1883 
1884  /*
1885  * Calculate the IMS cutoff factors
1886  */
1887  s_updateIMS_lnMolalityActCoeff();
1888 
1889  /*
1890  * Now do the main calculation.
1891  */
1892  s_updatePitzer_lnMolalityActCoeff();
1893 
1894  double xmolSolvent = moleFraction(m_indexSolvent);
1895  double xx = std::max(m_xmolSolventMIN, xmolSolvent);
1896  double lnActCoeffMolal0 = - log(xx) + (xx - 1.0)/xx;
1897  double lnxs = log(xx);
1898 
1899  for (size_t k = 1; k < m_kk; k++) {
1900  CROP_speciesCropped_[k] = 0;
1901  m_lnActCoeffMolal_Unscaled[k] += IMS_lnActCoeffMolal_[k];
1902  if (m_lnActCoeffMolal_Unscaled[k] > (CROP_ln_gamma_k_max- 2.5 *lnxs)) {
1903  CROP_speciesCropped_[k] = 2;
1904  m_lnActCoeffMolal_Unscaled[k] = CROP_ln_gamma_k_max - 2.5 * lnxs;
1905  }
1906  if (m_lnActCoeffMolal_Unscaled[k] < (CROP_ln_gamma_k_min - 2.5 *lnxs)) {
1907  // -1.0 and -1.5 caused multiple solutions
1908  CROP_speciesCropped_[k] = 2;
1909  m_lnActCoeffMolal_Unscaled[k] = CROP_ln_gamma_k_min - 2.5 * lnxs;
1910  }
1911  }
1912  CROP_speciesCropped_[0] = 0;
1913  m_lnActCoeffMolal_Unscaled[0] += (IMS_lnActCoeffMolal_[0] - lnActCoeffMolal0);
1914  if (m_lnActCoeffMolal_Unscaled[0] < CROP_ln_gamma_o_min) {
1915  CROP_speciesCropped_[0] = 2;
1916  m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_min;
1917  }
1918  if (m_lnActCoeffMolal_Unscaled[0] > CROP_ln_gamma_o_max) {
1919  CROP_speciesCropped_[0] = 2;
1920  // -0.5 caused multiple solutions
1921  m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_max;
1922  }
1923  if (m_lnActCoeffMolal_Unscaled[0] > CROP_ln_gamma_o_max - 0.5 * lnxs) {
1924  CROP_speciesCropped_[0] = 2;
1925  m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_max - 0.5 * lnxs;
1926  }
1927 
1928  /*
1929  * Now do the pH Scaling
1930  */
1931  s_updateScaling_pHScaling();
1932 }
1933 
1934 
1935 /*
1936  * Calculate cropped molalities
1937  */
1938 void HMWSoln::calcMolalitiesCropped() const
1939 {
1940  doublereal Imax = 0.0, Itmp;
1941  doublereal Iac_max;
1942  m_molalitiesAreCropped = false;
1943 
1944  for (size_t k = 0; k < m_kk; k++) {
1945  m_molalitiesCropped[k] = m_molalities[k];
1946  double charge = m_speciesCharge[k];
1947  Itmp = m_molalities[k] * charge * charge;
1948  if (Itmp > Imax) {
1949  Imax = Itmp;
1950  }
1951  }
1952 
1953  int cropMethod = 1;
1954 
1955 
1956  if (cropMethod == 0) {
1957 
1958  /*
1959  * Quick return
1960  */
1961  if (Imax < m_maxIionicStrength) {
1962  return;
1963  }
1964 
1965  m_molalitiesAreCropped = true;
1966 
1967  for (size_t i = 1; i < (m_kk - 1); i++) {
1968  double charge_i = m_speciesCharge[i];
1969  double abs_charge_i = fabs(charge_i);
1970  if (charge_i == 0.0) {
1971  continue;
1972  }
1973  for (size_t j = (i+1); j < m_kk; j++) {
1974  double charge_j = m_speciesCharge[j];
1975  double abs_charge_j = fabs(charge_j);
1976  /*
1977  * Find the counterIJ for the symmetric binary interaction
1978  */
1979  // n = m_kk*i + j;
1980  // counterIJ = m_CounterIJ[n];
1981  /*
1982  * Only loop over oppositely charge species
1983  */
1984  if (charge_i * charge_j < 0) {
1985  Iac_max = m_maxIionicStrength;
1986 
1987  if (m_molalitiesCropped[i] > m_molalitiesCropped[j]) {
1988  Imax = m_molalitiesCropped[i] * abs_charge_i * abs_charge_i;
1989  if (Imax > Iac_max) {
1990  m_molalitiesCropped[i] = Iac_max / (abs_charge_i * abs_charge_i);
1991  }
1992  Imax = m_molalitiesCropped[j] * fabs(abs_charge_j * abs_charge_i);
1993  if (Imax > Iac_max) {
1994  m_molalitiesCropped[j] = Iac_max / (abs_charge_j * abs_charge_i);
1995  }
1996  } else {
1997  Imax = m_molalitiesCropped[j] * abs_charge_j * abs_charge_j;
1998  if (Imax > Iac_max) {
1999  m_molalitiesCropped[j] = Iac_max / (abs_charge_j * abs_charge_j);
2000  }
2001  Imax = m_molalitiesCropped[i] * abs_charge_j * abs_charge_i;
2002  if (Imax > Iac_max) {
2003  m_molalitiesCropped[i] = Iac_max / (abs_charge_j * abs_charge_i);
2004  }
2005  }
2006  }
2007  }
2008  }
2009 
2010  /*
2011  * Do this loop 10 times until we have achieved charge neutrality
2012  * in the cropped molalities
2013  */
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++) {
2022  double charge_i = m_speciesCharge[i];
2023  if (charge_i < 0.0) {
2024  double anion_contrib = - m_molalitiesCropped[i] * charge_i;
2025  anion_charge += anion_contrib ;
2026  if (anion_contrib > anion_contrib_max) {
2027  anion_contrib_max = anion_contrib;
2028  anion_contrib_max_i = i;
2029  }
2030  } else if (charge_i > 0.0) {
2031  double cation_contrib = m_molalitiesCropped[i] * charge_i;
2032  cation_charge += cation_contrib ;
2033  if (cation_contrib > cation_contrib_max) {
2034  cation_contrib_max = cation_contrib;
2035  cation_contrib_max_i = i;
2036  }
2037  }
2038  }
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];
2042  double maxCrop = 0.66 * m_molalitiesCropped[cation_contrib_max_i];
2043  if (desiredCrop < maxCrop) {
2044  m_molalitiesCropped[cation_contrib_max_i] -= desiredCrop;
2045  break;
2046  } else {
2047  m_molalitiesCropped[cation_contrib_max_i] -= maxCrop;
2048  }
2049  } else if (total_charge < -1.0E-8) {
2050  double desiredCrop = total_charge/m_speciesCharge[anion_contrib_max_i];
2051  double maxCrop = 0.66 * m_molalitiesCropped[anion_contrib_max_i];
2052  if (desiredCrop < maxCrop) {
2053  m_molalitiesCropped[anion_contrib_max_i] -= desiredCrop;
2054  break;
2055  } else {
2056  m_molalitiesCropped[anion_contrib_max_i] -= maxCrop;
2057  }
2058  } else {
2059  break;
2060  }
2061  }
2062  }
2063 
2064  if (cropMethod == 1) {
2065  double* molF = DATA_PTR(m_gamma_tmp);
2066  getMoleFractions(molF);
2067  double xmolSolvent = molF[m_indexSolvent];
2068  if (xmolSolvent >= MC_X_o_cutoff_) {
2069  return;
2070  }
2071 
2072  m_molalitiesAreCropped = true;
2073 
2074  double poly = MC_apCut_ + MC_bpCut_ * xmolSolvent + MC_dpCut_* xmolSolvent * xmolSolvent;
2075  double p = xmolSolvent + MC_epCut_ + exp(- xmolSolvent/ MC_cpCut_) * poly;
2076  double denomInv = 1.0/ (m_Mnaught * p);
2077 
2078  for (size_t k = 0; k < m_kk; k++) {
2079  m_molalitiesCropped[k] = molF[k] * denomInv ;
2080  }
2081 
2082  // Do a further check to see if the Ionic strength is below a max value
2083  // Reduce the molalities to enforce this. Note, this algorithm preserves
2084  // the charge neutrality of the solution after cropping.
2085  Itmp = 0.0;
2086  for (size_t k = 0; k < m_kk; k++) {
2087  double charge = m_speciesCharge[k];
2088  Itmp += m_molalitiesCropped[k] * charge * charge;
2089  }
2090  if (Itmp > m_maxIionicStrength) {
2091  double ratio = Itmp / m_maxIionicStrength;
2092  for (size_t k = 0; k < m_kk; k++) {
2093  double charge = m_speciesCharge[k];
2094  if (charge != 0.0) {
2095  m_molalitiesCropped[k] *= ratio;
2096  }
2097  }
2098  }
2099  }
2100 
2101 }
2102 
2103 /*
2104  * Set up a counter variable for keeping track of symmetric binary
2105  * interactions amongst the solute species.
2106  *
2107  * n = m_kk*i + j
2108  * m_Counter[n] = counter
2109  */
2110 void HMWSoln::counterIJ_setup(void) const
2111 {
2112  size_t n, nc, i, j;
2113  m_CounterIJ.resize(m_kk * m_kk);
2114  int counter = 0;
2115  for (i = 0; i < m_kk; i++) {
2116  n = i;
2117  nc = m_kk * i;
2118  m_CounterIJ[n] = 0;
2119  m_CounterIJ[nc] = 0;
2120  }
2121  for (i = 1; i < (m_kk - 1); i++) {
2122  n = m_kk * i + i;
2123  m_CounterIJ[n] = 0;
2124  for (j = (i+1); j < m_kk; j++) {
2125  n = m_kk * j + i;
2126  nc = m_kk * i + j;
2127  counter++;
2128  m_CounterIJ[n] = counter;
2129  m_CounterIJ[nc] = counter;
2130  }
2131  }
2132 }
2133 
2134 /*
2135  * Calculates the Pitzer coefficients' dependence on the
2136  * temperature. It will also calculate the temperature
2137  * derivatives of the coefficients, as they are important
2138  * in the calculation of the latent heats and the
2139  * heat capacities of the mixtures.
2140  *
2141  * @param doDerivs If >= 1, then the routine will calculate
2142  * the first derivative. If >= 2, the
2143  * routine will calculate the first and second
2144  * temperature derivative.
2145  * default = 2
2146  */
2147 void HMWSoln::s_updatePitzer_CoeffWRTemp(int doDerivs) const
2148 {
2149 
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;
2156  double T = temperature();
2157  double Tr = m_TempPitzerRef;
2158  double tinv = 0.0, tln = 0.0, tlin = 0.0, tquad = 0.0;
2159  if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
2160  tlin = T - Tr;
2161  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
2162  tlin = T - Tr;
2163  tquad = T * T - Tr * Tr;
2164  tln = log(T/ Tr);
2165  tinv = 1.0/T - 1.0/Tr;
2166  }
2167 
2168  for (i = 1; i < (m_kk - 1); i++) {
2169  for (j = (i+1); j < m_kk; j++) {
2170 
2171  /*
2172  * Find the counterIJ for the symmetric binary interaction
2173  */
2174  n = m_kk*i + j;
2175  counterIJ = m_CounterIJ[n];
2176 
2177  beta0MX_coeff = m_Beta0MX_ij_coeff.ptrColumn(counterIJ);
2178  beta1MX_coeff = m_Beta1MX_ij_coeff.ptrColumn(counterIJ);
2179  beta2MX_coeff = m_Beta2MX_ij_coeff.ptrColumn(counterIJ);
2180  CphiMX_coeff = m_CphiMX_ij_coeff.ptrColumn(counterIJ);
2181  Theta_coeff = m_Theta_ij_coeff.ptrColumn(counterIJ);
2182 
2183  switch (m_formPitzerTemp) {
2184  case PITZER_TEMP_CONSTANT:
2185  break;
2186  case PITZER_TEMP_LINEAR:
2187 
2188  m_Beta0MX_ij[counterIJ] = beta0MX_coeff[0]
2189  + beta0MX_coeff[1]*tlin;
2190  m_Beta0MX_ij_L[counterIJ] = beta0MX_coeff[1];
2191  m_Beta0MX_ij_LL[counterIJ] = 0.0;
2192 
2193  m_Beta1MX_ij[counterIJ] = beta1MX_coeff[0]
2194  + beta1MX_coeff[1]*tlin;
2195  m_Beta1MX_ij_L[counterIJ] = beta1MX_coeff[1];
2196  m_Beta1MX_ij_LL[counterIJ] = 0.0;
2197 
2198  m_Beta2MX_ij[counterIJ] = beta2MX_coeff[0]
2199  + beta2MX_coeff[1]*tlin;
2200  m_Beta2MX_ij_L[counterIJ] = beta2MX_coeff[1];
2201  m_Beta2MX_ij_LL[counterIJ] = 0.0;
2202 
2203  m_CphiMX_ij[counterIJ] = CphiMX_coeff[0]
2204  + CphiMX_coeff[1]*tlin;
2205  m_CphiMX_ij_L[counterIJ] = CphiMX_coeff[1];
2206  m_CphiMX_ij_LL[counterIJ] = 0.0;
2207 
2208  m_Theta_ij[counterIJ] = Theta_coeff[0] + Theta_coeff[1]*tlin;
2209  m_Theta_ij_L[counterIJ] = Theta_coeff[1];
2210  m_Theta_ij_LL[counterIJ] = 0.0;
2211 
2212  break;
2213 
2214  case PITZER_TEMP_COMPLEX1:
2215  m_Beta0MX_ij[counterIJ] = beta0MX_coeff[0]
2216  + beta0MX_coeff[1]*tlin
2217  + beta0MX_coeff[2]*tquad
2218  + beta0MX_coeff[3]*tinv
2219  + beta0MX_coeff[4]*tln;
2220 
2221  m_Beta1MX_ij[counterIJ] = beta1MX_coeff[0]
2222  + beta1MX_coeff[1]*tlin
2223  + beta1MX_coeff[2]*tquad
2224  + beta1MX_coeff[3]*tinv
2225  + beta1MX_coeff[4]*tln;
2226 
2227  m_Beta2MX_ij[counterIJ] = beta2MX_coeff[0]
2228  + beta2MX_coeff[1]*tlin
2229  + beta2MX_coeff[2]*tquad
2230  + beta2MX_coeff[3]*tinv
2231  + beta2MX_coeff[4]*tln;
2232 
2233  m_CphiMX_ij[counterIJ] = CphiMX_coeff[0]
2234  + CphiMX_coeff[1]*tlin
2235  + CphiMX_coeff[2]*tquad
2236  + CphiMX_coeff[3]*tinv
2237  + CphiMX_coeff[4]*tln;
2238 
2239  m_Theta_ij[counterIJ] = Theta_coeff[0]
2240  + Theta_coeff[1]*tlin
2241  + Theta_coeff[2]*tquad
2242  + Theta_coeff[3]*tinv
2243  + Theta_coeff[4]*tln;
2244 
2245  m_Beta0MX_ij_L[counterIJ] = beta0MX_coeff[1]
2246  + beta0MX_coeff[2]*2.0*T
2247  - beta0MX_coeff[3]/(T*T)
2248  + beta0MX_coeff[4]/T;
2249 
2250  m_Beta1MX_ij_L[counterIJ] = beta1MX_coeff[1]
2251  + beta1MX_coeff[2]*2.0*T
2252  - beta1MX_coeff[3]/(T*T)
2253  + beta1MX_coeff[4]/T;
2254 
2255  m_Beta2MX_ij_L[counterIJ] = beta2MX_coeff[1]
2256  + beta2MX_coeff[2]*2.0*T
2257  - beta2MX_coeff[3]/(T*T)
2258  + beta2MX_coeff[4]/T;
2259 
2260  m_CphiMX_ij_L[counterIJ] = CphiMX_coeff[1]
2261  + CphiMX_coeff[2]*2.0*T
2262  - CphiMX_coeff[3]/(T*T)
2263  + CphiMX_coeff[4]/T;
2264 
2265  m_Theta_ij_L[counterIJ] = Theta_coeff[1]
2266  + Theta_coeff[2]*2.0*T
2267  - Theta_coeff[3]/(T*T)
2268  + Theta_coeff[4]/T;
2269 
2270  doDerivs = 2;
2271  if (doDerivs > 1) {
2272  m_Beta0MX_ij_LL[counterIJ] =
2273  + beta0MX_coeff[2]*2.0
2274  + 2.0*beta0MX_coeff[3]/(T*T*T)
2275  - beta0MX_coeff[4]/(T*T);
2276 
2277  m_Beta1MX_ij_LL[counterIJ] =
2278  + beta1MX_coeff[2]*2.0
2279  + 2.0*beta1MX_coeff[3]/(T*T*T)
2280  - beta1MX_coeff[4]/(T*T);
2281 
2282  m_Beta2MX_ij_LL[counterIJ] =
2283  + beta2MX_coeff[2]*2.0
2284  + 2.0*beta2MX_coeff[3]/(T*T*T)
2285  - beta2MX_coeff[4]/(T*T);
2286 
2287  m_CphiMX_ij_LL[counterIJ] =
2288  + CphiMX_coeff[2]*2.0
2289  + 2.0*CphiMX_coeff[3]/(T*T*T)
2290  - CphiMX_coeff[4]/(T*T);
2291 
2292  m_Theta_ij_LL[counterIJ] =
2293  + Theta_coeff[2]*2.0
2294  + 2.0*Theta_coeff[3]/(T*T*T)
2295  - Theta_coeff[4]/(T*T);
2296  }
2297 
2298 #ifdef DEBUG_HKM
2299  /*
2300  * Turn terms off for debugging
2301  */
2302  //m_Beta0MX_ij_L[counterIJ] = 0;
2303  //m_Beta0MX_ij_LL[counterIJ] = 0;
2304  //m_Beta1MX_ij_L[counterIJ] = 0;
2305  //m_Beta1MX_ij_LL[counterIJ] = 0;
2306  //m_CphiMX_ij_L[counterIJ] = 0;
2307  //m_CphiMX_ij_LL[counterIJ] = 0;
2308 #endif
2309  break;
2310  }
2311  }
2312  }
2313 
2314  // Lambda interactions and Mu_nnn
2315  // i must be neutral for this term to be nonzero. We take advantage of this
2316  // here to lower the operation count.
2317  for (i = 1; i < m_kk; i++) {
2318  if (m_speciesCharge[i] == 0.0) {
2319  for (j = 1; j < m_kk; j++) {
2320  n = i * m_kk + j;
2321  const double* Lambda_coeff = m_Lambda_nj_coeff.ptrColumn(n);
2322  switch (m_formPitzerTemp) {
2323  case PITZER_TEMP_CONSTANT:
2324  m_Lambda_nj(i,j) = Lambda_coeff[0];
2325  break;
2326  case PITZER_TEMP_LINEAR:
2327  m_Lambda_nj(i,j) = Lambda_coeff[0] + Lambda_coeff[1]*tlin;
2328  m_Lambda_nj_L(i,j) = Lambda_coeff[1];
2329  m_Lambda_nj_LL(i,j) = 0.0;
2330  break;
2331  case PITZER_TEMP_COMPLEX1:
2332  m_Lambda_nj(i,j) = Lambda_coeff[0]
2333  + Lambda_coeff[1]*tlin
2334  + Lambda_coeff[2]*tquad
2335  + Lambda_coeff[3]*tinv
2336  + Lambda_coeff[4]*tln;
2337 
2338  m_Lambda_nj_L(i,j) = Lambda_coeff[1]
2339  + Lambda_coeff[2]*2.0*T
2340  - Lambda_coeff[3]/(T*T)
2341  + Lambda_coeff[4]/T;
2342 
2343  m_Lambda_nj_LL(i,j) =
2344  Lambda_coeff[2]*2.0
2345  + 2.0*Lambda_coeff[3]/(T*T*T)
2346  - Lambda_coeff[4]/(T*T);
2347  }
2348 
2349  if (j == i) {
2350  const double* Mu_coeff = m_Mu_nnn_coeff.ptrColumn(i);
2351  switch (m_formPitzerTemp) {
2352  case PITZER_TEMP_CONSTANT:
2353  m_Mu_nnn[i] = Mu_coeff[0];
2354  break;
2355  case PITZER_TEMP_LINEAR:
2356  m_Mu_nnn[i] = Mu_coeff[0] + Mu_coeff[1]*tlin;
2357  m_Mu_nnn_L[i] = Mu_coeff[1];
2358  m_Mu_nnn_LL[i] = 0.0;
2359  break;
2360  case PITZER_TEMP_COMPLEX1:
2361  m_Mu_nnn[i] = Mu_coeff[0]
2362  + Mu_coeff[1]*tlin
2363  + Mu_coeff[2]*tquad
2364  + Mu_coeff[3]*tinv
2365  + Mu_coeff[4]*tln;
2366 
2367  m_Mu_nnn_L[i] = Mu_coeff[1]
2368  + Mu_coeff[2]*2.0*T
2369  - Mu_coeff[3]/(T*T)
2370  + Mu_coeff[4]/T;
2371 
2372  m_Mu_nnn_LL[i] =
2373  Mu_coeff[2]*2.0
2374  + 2.0*Mu_coeff[3]/(T*T*T)
2375  - Mu_coeff[4]/(T*T);
2376  }
2377  }
2378  }
2379  }
2380  }
2381 
2382 
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 ;
2387  const double* Psi_coeff = m_Psi_ijk_coeff.ptrColumn(n);
2388  switch (m_formPitzerTemp) {
2389  case PITZER_TEMP_CONSTANT:
2390  m_Psi_ijk[n] = Psi_coeff[0];
2391  break;
2392  case PITZER_TEMP_LINEAR:
2393  m_Psi_ijk[n] = Psi_coeff[0] + Psi_coeff[1]*tlin;
2394  m_Psi_ijk_L[n] = Psi_coeff[1];
2395  m_Psi_ijk_LL[n] = 0.0;
2396  break;
2397  case PITZER_TEMP_COMPLEX1:
2398  m_Psi_ijk[n] = Psi_coeff[0]
2399  + Psi_coeff[1]*tlin
2400  + Psi_coeff[2]*tquad
2401  + Psi_coeff[3]*tinv
2402  + Psi_coeff[4]*tln;
2403 
2404  m_Psi_ijk_L[n] = Psi_coeff[1]
2405  + Psi_coeff[2]*2.0*T
2406  - Psi_coeff[3]/(T*T)
2407  + Psi_coeff[4]/T;
2408 
2409  m_Psi_ijk_LL[n] =
2410  Psi_coeff[2]*2.0
2411  + 2.0*Psi_coeff[3]/(T*T*T)
2412  - Psi_coeff[4]/(T*T);
2413  }
2414  }
2415  }
2416  }
2417 
2418 }
2419 
2420 /*
2421  * Calculate the Pitzer portion of the activity coefficients.
2422  *
2423  * This is the main routine in the whole module. It calculates the
2424  * molality based activity coefficients for the solutes, and
2425  * the activity of water.
2426  */
2427 void HMWSoln::
2428 s_updatePitzer_lnMolalityActCoeff() const
2429 {
2430 
2431  /*
2432  * HKM -> Assumption is made that the solvent is
2433  * species 0.
2434  */
2435  if (m_indexSolvent != 0) {
2436  printf("Wrong index solvent value!\n");
2437  exit(EXIT_FAILURE);
2438  }
2439 
2440 #ifdef DEBUG_MODE
2441  int printE = 0;
2442  if (temperature() == 323.15) {
2443  printE = 0;
2444  }
2445 #endif
2446  std::string sni, snj, snk;
2447 
2448  /*
2449  * Use the CROPPED molality of the species in solution.
2450  */
2451  const double* molality = DATA_PTR(m_molalitiesCropped);
2452  /*
2453  * These are the charges of the species accessed from class Phase
2454  */
2455  const double* charge = DATA_PTR(m_speciesCharge);
2456 
2457  /*
2458  * These are data inputs about the Pitzer correlation. They come
2459  * from the input file for the Pitzer model.
2460  */
2461  const double* beta0MX = DATA_PTR(m_Beta0MX_ij);
2462  const double* beta1MX = DATA_PTR(m_Beta1MX_ij);
2463  const double* beta2MX = DATA_PTR(m_Beta2MX_ij);
2464  const double* CphiMX = DATA_PTR(m_CphiMX_ij);
2465  const double* thetaij = DATA_PTR(m_Theta_ij);
2466  const double* alpha1MX = DATA_PTR(m_Alpha1MX_ij);
2467  const double* alpha2MX = DATA_PTR(m_Alpha2MX_ij);
2468 
2469  const double* psi_ijk = DATA_PTR(m_Psi_ijk);
2470  //n = k + j * m_kk + i * m_kk * m_kk;
2471 
2472 
2473  double* gamma_Unscaled = DATA_PTR(m_gamma_tmp);
2474  /*
2475  * Local variables defined by Coltrin
2476  */
2477  double etheta[5][5], etheta_prime[5][5], sqrtIs;
2478  /*
2479  * Molality based ionic strength of the solution
2480  */
2481  double Is = 0.0;
2482  /*
2483  * Molarcharge of the solution: In Pitzer's notation,
2484  * this is his variable called "Z".
2485  */
2486  double molarcharge = 0.0;
2487  /*
2488  * molalitysum is the sum of the molalities over all solutes,
2489  * even those with zero charge.
2490  */
2491  double molalitysumUncropped = 0.0;
2492 
2493  double* gfunc = DATA_PTR(m_gfunc_IJ);
2494  double* g2func = DATA_PTR(m_g2func_IJ);
2495  double* hfunc = DATA_PTR(m_hfunc_IJ);
2496  double* h2func = DATA_PTR(m_h2func_IJ);
2497  double* BMX = DATA_PTR(m_BMX_IJ);
2498  double* BprimeMX = DATA_PTR(m_BprimeMX_IJ);
2499  double* BphiMX = DATA_PTR(m_BphiMX_IJ);
2500  double* Phi = DATA_PTR(m_Phi_IJ);
2501  double* Phiprime = DATA_PTR(m_Phiprime_IJ);
2502  double* Phiphi = DATA_PTR(m_PhiPhi_IJ);
2503  double* CMX = DATA_PTR(m_CMX_IJ);
2504 
2505 
2506  double x1, x2;
2507  double Aphi, F, zsqF;
2508  double sum1, sum2, sum3, sum4, sum5, term1;
2509  double sum_m_phi_minus_1, osmotic_coef, lnwateract;
2510 
2511  int z1, z2;
2512  size_t n, i, j, k, m, counterIJ, counterIJ2;
2513 
2514 #ifdef DEBUG_MODE
2515  if (m_debugCalc) {
2516  printf("\n Debugging information from hmw_act \n");
2517  }
2518 #endif
2519  /*
2520  * Make sure the counter variables are setup
2521  */
2522  counterIJ_setup();
2523 
2524  /*
2525  * ---------- Calculate common sums over solutes ---------------------
2526  */
2527  for (n = 1; n < m_kk; n++) {
2528  // ionic strength
2529  Is += charge[n] * charge[n] * molality[n];
2530  // total molar charge
2531  molarcharge += fabs(charge[n]) * molality[n];
2532  molalitysumUncropped += m_molalities[n];
2533  }
2534  Is *= 0.5;
2535 
2536  /*
2537  * Store the ionic molality in the object for reference.
2538  */
2539  m_IionicMolality = Is;
2540  sqrtIs = sqrt(Is);
2541 #ifdef DEBUG_MODE
2542  if (m_debugCalc) {
2543  printf(" Step 1: \n");
2544  printf(" ionic strenth = %14.7le \n total molar "
2545  "charge = %14.7le \n", Is, molarcharge);
2546  }
2547 #endif
2548 
2549  /*
2550  * The following call to calc_lambdas() calculates all 16 elements
2551  * of the elambda and elambda1 arrays, given the value of the
2552  * ionic strength (Is)
2553  */
2554  calc_lambdas(Is);
2555 
2556  /*
2557  * ----- Step 2: Find the coefficients E-theta and -------------------
2558  * E-thetaprime for all combinations of positive
2559  * unlike charges up to 4
2560  */
2561 #ifdef DEBUG_MODE
2562  if (m_debugCalc) {
2563  printf(" Step 2: \n");
2564  }
2565 #endif
2566  for (z1 = 1; z1 <=4; z1++) {
2567  for (z2 =1; z2 <=4; z2++) {
2568  calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
2569 #ifdef DEBUG_MODE
2570  if (m_debugCalc) {
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]);
2573  }
2574 #endif
2575  }
2576  }
2577 
2578 #ifdef DEBUG_MODE
2579  if (m_debugCalc) {
2580  printf(" Step 3: \n");
2581  printf(" Species Species g(x) "
2582  " hfunc(x) \n");
2583  }
2584 #endif
2585 
2586  /*
2587  *
2588  * calculate g(x) and hfunc(x) for each cation-anion pair MX
2589  * In the original literature, hfunc, was called gprime. However,
2590  * it's not the derivative of g(x), so I renamed it.
2591  */
2592  for (i = 1; i < (m_kk - 1); i++) {
2593  for (j = (i+1); j < m_kk; j++) {
2594  /*
2595  * Find the counterIJ for the symmetric binary interaction
2596  */
2597  n = m_kk*i + j;
2598  counterIJ = m_CounterIJ[n];
2599 
2600  /*
2601  * Only loop over oppositely charge species
2602  */
2603  if (charge[i]*charge[j] < 0) {
2604  /*
2605  * x is a reduced function variable
2606  */
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);
2612  } else {
2613  gfunc[counterIJ] = 0.0;
2614  hfunc[counterIJ] = 0.0;
2615  }
2616 
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);
2623  } else {
2624  g2func[counterIJ] = 0.0;
2625  h2func[counterIJ] = 0.0;
2626  }
2627  }
2628  } else {
2629  gfunc[counterIJ] = 0.0;
2630  hfunc[counterIJ] = 0.0;
2631  }
2632 #ifdef DEBUG_MODE
2633  if (m_debugCalc) {
2634  sni = speciesName(i);
2635  snj = speciesName(j);
2636  printf(" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(),
2637  gfunc[counterIJ], hfunc[counterIJ]);
2638  }
2639 #endif
2640  }
2641  }
2642 
2643  /*
2644  * --------- SUBSECTION TO CALCULATE BMX, BprimeMX, BphiMX ----------
2645  * --------- Agrees with Pitzer, Eq. (49), (51), (55)
2646  */
2647 #ifdef DEBUG_MODE
2648  if (m_debugCalc) {
2649  printf(" Step 4: \n");
2650  printf(" Species Species BMX "
2651  "BprimeMX BphiMX \n");
2652  }
2653 #endif
2654 
2655  for (i = 1; i < m_kk - 1; i++) {
2656  for (j = i+1; j < m_kk; j++) {
2657  /*
2658  * Find the counterIJ for the symmetric binary interaction
2659  */
2660  n = m_kk*i + j;
2661  counterIJ = m_CounterIJ[n];
2662 
2663 #ifdef DEBUG_MODE
2664  if (printE) {
2665  if (counterIJ == 2) {
2666  printf("%s %s\n", speciesName(i).c_str(),
2667  speciesName(j).c_str());
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]);
2671  }
2672  }
2673 #endif
2674  /*
2675  * both species have a non-zero charge, and one is positive
2676  * and the other is negative
2677  */
2678  if (charge[i]*charge[j] < 0.0) {
2679  BMX[counterIJ] = beta0MX[counterIJ]
2680  + beta1MX[counterIJ] * gfunc[counterIJ]
2681  + beta2MX[counterIJ] * g2func[counterIJ];
2682 #ifdef DEBUG_MODE
2683  if (m_debugCalc) {
2684  printf("%d %g: %g %g %g %g\n",
2685  counterIJ, BMX[counterIJ], beta0MX[counterIJ],
2686  beta1MX[counterIJ], beta2MX[counterIJ], gfunc[counterIJ]);
2687  }
2688 #endif
2689  if (Is > 1.0E-150) {
2690  BprimeMX[counterIJ] = (beta1MX[counterIJ] * hfunc[counterIJ]/Is +
2691  beta2MX[counterIJ] * h2func[counterIJ]/Is);
2692  } else {
2693  BprimeMX[counterIJ] = 0.0;
2694  }
2695  BphiMX[counterIJ] = BMX[counterIJ] + Is*BprimeMX[counterIJ];
2696  } else {
2697  BMX[counterIJ] = 0.0;
2698  BprimeMX[counterIJ] = 0.0;
2699  BphiMX[counterIJ] = 0.0;
2700  }
2701 #ifdef DEBUG_MODE
2702  if (m_debugCalc) {
2703  sni = speciesName(i);
2704  snj = speciesName(j);
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]);
2708  }
2709 #endif
2710  }
2711  }
2712 
2713  /*
2714  * --------- SUBSECTION TO CALCULATE CMX ----------
2715  * --------- Agrees with Pitzer, Eq. (53).
2716  */
2717 #ifdef DEBUG_MODE
2718  if (m_debugCalc) {
2719  printf(" Step 5: \n");
2720  printf(" Species Species CMX \n");
2721  }
2722 #endif
2723  for (i = 1; i < m_kk-1; i++) {
2724  for (j = i+1; j < m_kk; j++) {
2725  /*
2726  * Find the counterIJ for the symmetric binary interaction
2727  */
2728  n = m_kk*i + j;
2729  counterIJ = m_CounterIJ[n];
2730  /*
2731  * both species have a non-zero charge, and one is positive
2732  * and the other is negative
2733  */
2734  if (charge[i]*charge[j] < 0.0) {
2735  CMX[counterIJ] = CphiMX[counterIJ]/
2736  (2.0* sqrt(fabs(charge[i]*charge[j])));
2737  } else {
2738  CMX[counterIJ] = 0.0;
2739  }
2740 #ifdef DEBUG_MODE
2741  if (printE) {
2742  if (counterIJ == 2) {
2743  printf("%s %s\n", speciesName(i).c_str(),
2744  speciesName(j).c_str());
2745  printf("CphiMX[%d] = %g\n", counterIJ, CphiMX[counterIJ]);
2746  }
2747  }
2748 #endif
2749 #ifdef DEBUG_MODE
2750  if (m_debugCalc) {
2751  sni = speciesName(i);
2752  snj = speciesName(j);
2753  printf(" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(),
2754  CMX[counterIJ]);
2755  }
2756 #endif
2757  }
2758  }
2759 
2760  /*
2761  * ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ----------
2762  * --------- Agrees with Pitzer, Eq. 72, 73, 74
2763  */
2764 #ifdef DEBUG_MODE
2765  if (m_debugCalc) {
2766  printf(" Step 6: \n");
2767  printf(" Species Species Phi_ij "
2768  " Phiprime_ij Phi^phi_ij \n");
2769  }
2770 #endif
2771  for (i = 1; i < m_kk-1; i++) {
2772  for (j = i+1; j < m_kk; j++) {
2773  /*
2774  * Find the counterIJ for the symmetric binary interaction
2775  */
2776  n = m_kk*i + j;
2777  counterIJ = m_CounterIJ[n];
2778  /*
2779  * both species have a non-zero charge, and one is positive
2780  * and the other is negative
2781  */
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];
2788  } else {
2789  Phi[counterIJ] = 0.0;
2790  Phiprime[counterIJ] = 0.0;
2791  Phiphi[counterIJ] = 0.0;
2792  }
2793 #ifdef DEBUG_MODE
2794  if (m_debugCalc) {
2795  sni = speciesName(i);
2796  snj = speciesName(j);
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]);
2800  }
2801 #endif
2802  }
2803  }
2804 
2805  /*
2806  * ------------- SUBSECTION FOR CALCULATION OF F ----------------------
2807  * ------------ Agrees with Pitzer Eqn. (65) --------------------------
2808  */
2809 #ifdef DEBUG_MODE
2810  if (m_debugCalc) {
2811  printf(" Step 7: \n");
2812  }
2813 #endif
2814  // A_Debye_Huckel = 0.5092; (units = sqrt(kg/gmol))
2815  // A_Debye_Huckel = 0.5107; <- This value is used to match GWB data
2816  // ( A * ln(10) = 1.17593)
2817  // Aphi = A_Debye_Huckel * 2.30258509 / 3.0;
2818  Aphi = m_A_Debye / 3.0;
2819  F = -Aphi * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
2820  + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
2821 #ifdef DEBUG_MODE
2822  if (printE) {
2823  printf("Aphi = %20.13g\n", Aphi);
2824  }
2825 #endif
2826 #ifdef DEBUG_MODE
2827  if (m_debugCalc) {
2828  printf(" initial value of F = %10.6f \n", F);
2829  }
2830 #endif
2831  for (i = 1; i < m_kk-1; i++) {
2832  for (j = i+1; j < m_kk; j++) {
2833  /*
2834  * Find the counterIJ for the symmetric binary interaction
2835  */
2836  n = m_kk*i + j;
2837  counterIJ = m_CounterIJ[n];
2838  /*
2839  * both species have a non-zero charge, and one is positive
2840  * and the other is negative
2841  */
2842  if (charge[i]*charge[j] < 0) {
2843  F = F + molality[i]*molality[j] * BprimeMX[counterIJ];
2844  }
2845  /*
2846  * Both species have a non-zero charge, and they
2847  * have the same sign
2848  */
2849  if (charge[i]*charge[j] > 0) {
2850  F = F + molality[i]*molality[j] * Phiprime[counterIJ];
2851  }
2852 #ifdef DEBUG_MODE
2853  if (m_debugCalc) {
2854  printf(" F = %10.6f \n", F);
2855  }
2856 #endif
2857  }
2858  }
2859 #ifdef DEBUG_MODE
2860  if (m_debugCalc) {
2861  printf(" Step 8: Summing in All Contributions to Activity Coefficients \n");
2862  }
2863 #endif
2864 
2865  for (i = 1; i < m_kk; i++) {
2866 
2867  /*
2868  * -------- SUBSECTION FOR CALCULATING THE ACTCOEFF FOR CATIONS -----
2869  * -------- -> equations agree with my notes, Eqn. (118).
2870  * -> Equations agree with Pitzer, eqn.(63)
2871  */
2872  if (charge[i] > 0.0) {
2873 
2874 #ifdef DEBUG_MODE
2875  if (m_debugCalc) {
2876  sni = speciesName(i);
2877  printf(" Contributions to ln(ActCoeff_%s):\n", sni.c_str());
2878  }
2879 #endif
2880  // species i is the cation (positive) to calc the actcoeff
2881  zsqF = charge[i]*charge[i]*F;
2882 #ifdef DEBUG_MODE
2883  if (m_debugCalc) {
2884  printf(" Unary term: z*z*F = %10.5f\n", zsqF);
2885  }
2886 #endif
2887  sum1 = 0.0;
2888  sum2 = 0.0;
2889  sum3 = 0.0;
2890  sum4 = 0.0;
2891  sum5 = 0.0;
2892  for (j = 1; j < m_kk; j++) {
2893  /*
2894  * Find the counterIJ for the symmetric binary interaction
2895  */
2896  n = m_kk*i + j;
2897  counterIJ = m_CounterIJ[n];
2898 
2899  if (charge[j] < 0.0) {
2900  // sum over all anions
2901  sum1 = sum1 + molality[j]*
2902  (2.0*BMX[counterIJ] + molarcharge*CMX[counterIJ]);
2903 #ifdef DEBUG_MODE
2904  if (m_debugCalc) {
2905  snj = speciesName(j) + ":";
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]);
2910  }
2911 #endif
2912  if (j < m_kk-1) {
2913  /*
2914  * This term is the ternary interaction involving the
2915  * non-duplicate sum over double anions, j, k, with
2916  * respect to the cation, i.
2917  */
2918  for (k = j+1; k < m_kk; k++) {
2919  // an inner sum over all anions
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];
2923 #ifdef DEBUG_MODE
2924  if (m_debugCalc) {
2925  if (psi_ijk[n] != 0.0) {
2926  snj = speciesName(j) + "," + speciesName(k) + ":";
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]);
2929  }
2930  }
2931 #endif
2932  }
2933  }
2934  }
2935  }
2936 
2937 
2938  if (charge[j] > 0.0) {
2939  // sum over all cations
2940  if (j != i) {
2941  sum2 = sum2 + molality[j]*(2.0*Phi[counterIJ]);
2942 #ifdef DEBUG_MODE
2943  if (m_debugCalc) {
2944  if ((molality[j] * Phi[counterIJ])!= 0.0) {
2945  snj = speciesName(j) + ":";
2946  printf(" Phi term with %-12s 2 m_j Phi_cc = %10.5f\n", snj.c_str(),
2947  molality[j]*(2.0*Phi[counterIJ]));
2948  }
2949  }
2950 #endif
2951  }
2952  for (k = 1; k < m_kk; k++) {
2953  if (charge[k] < 0.0) {
2954  // two inner sums over anions
2955 
2956  n = k + j * m_kk + i * m_kk * m_kk;
2957  sum2 = sum2 + molality[j]*molality[k]*psi_ijk[n];
2958 #ifdef DEBUG_MODE
2959  if (m_debugCalc) {
2960  if (psi_ijk[n] != 0.0) {
2961  snj = speciesName(j) + "," + speciesName(k) + ":";
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]);
2964  }
2965  }
2966 #endif
2967  /*
2968  * Find the counterIJ for the j,k interaction
2969  */
2970  n = m_kk*j + k;
2971  counterIJ2 = m_CounterIJ[n];
2972  sum4 = sum4 + (fabs(charge[i])*
2973  molality[j]*molality[k]*CMX[counterIJ2]);
2974 #ifdef DEBUG_MODE
2975  if (m_debugCalc) {
2976  if ((molality[j]*molality[k]*CMX[counterIJ2]) != 0.0) {
2977  snj = speciesName(j) + "," + speciesName(k) + ":";
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]);
2980  }
2981  }
2982 #endif
2983  }
2984  }
2985  }
2986 
2987  /*
2988  * Handle neutral j species
2989  */
2990  if (charge[j] == 0) {
2991  sum5 = sum5 + molality[j]*2.0*m_Lambda_nj(j,i);
2992 #ifdef DEBUG_MODE
2993  if (m_debugCalc) {
2994  if ((molality[j]*2.0*m_Lambda_nj(j,i)) != 0.0) {
2995  snj = speciesName(j) + ":";
2996  printf(" Lambda term with %-12s 2 m_j lam_ji = %10.5f\n", snj.c_str(),
2997  molality[j]*2.0*m_Lambda_nj(j,i));
2998  }
2999  }
3000 #endif
3001  /*
3002  * Zeta interaction term
3003  */
3004  for (size_t k = 1; k < m_kk; k++) {
3005  if (charge[k] < 0.0) {
3006  size_t izeta = j;
3007  size_t jzeta = i;
3008  n = izeta * m_kk * m_kk + jzeta * m_kk + k;
3009  double zeta = psi_ijk[n];
3010  if (zeta != 0.0) {
3011  sum5 = sum5 + molality[j]*molality[k]*zeta;
3012 #ifdef DEBUG_MODE
3013  if (m_debugCalc) {
3014  snj = speciesName(j) + "," + speciesName(k) + ":";
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]);
3017  }
3018 #endif
3019  }
3020  }
3021  }
3022  }
3023  }
3024  /*
3025  * Add all of the contributions up to yield the log of the
3026  * solute activity coefficients (molality scale)
3027  */
3028  m_lnActCoeffMolal_Unscaled[i] = zsqF + sum1 + sum2 + sum3 + sum4 + sum5;
3029  gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]);
3030 #ifdef DEBUG_MODE
3031  if (m_debugCalc) {
3032  sni = speciesName(i);
3033  printf(" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f \n",
3034  sni.c_str(), m_lnActCoeffMolal_Unscaled[i], gamma_Unscaled[i]);
3035  }
3036 #endif
3037  }
3038 
3039  /*
3040  * -------- SUBSECTION FOR CALCULATING THE ACTCOEFF FOR ANIONS ------
3041  * -------- -> equations agree with my notes, Eqn. (119).
3042  * -> Equations agree with Pitzer, eqn.(64)
3043  */
3044  if (charge[i] < 0) {
3045 
3046 #ifdef DEBUG_MODE
3047  if (m_debugCalc) {
3048  sni = speciesName(i);
3049  printf(" Contributions to ln(ActCoeff_%s):\n", sni.c_str());
3050  }
3051 #endif
3052 
3053  // species i is an anion (negative)
3054  zsqF = charge[i]*charge[i]*F;
3055 #ifdef DEBUG_MODE
3056  if (m_debugCalc) {
3057  printf(" Unary term: z*z*F = %10.5f\n", zsqF);
3058  }
3059 #endif
3060  sum1 = 0.0;
3061  sum2 = 0.0;
3062  sum3 = 0.0;
3063  sum4 = 0.0;
3064  sum5 = 0.0;
3065  for (j = 1; j < m_kk; j++) {
3066  /*
3067  * Find the counterIJ for the symmetric binary interaction
3068  */
3069  n = m_kk*i + j;
3070  counterIJ = m_CounterIJ[n];
3071 
3072  /*
3073  * For Anions, do the cation interactions.
3074  */
3075  if (charge[j] > 0) {
3076  sum1 = sum1 + molality[j]*
3077  (2.0*BMX[counterIJ]+molarcharge*CMX[counterIJ]);
3078 #ifdef DEBUG_MODE
3079  if (m_debugCalc) {
3080  snj = speciesName(j) + ":";
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]);
3085  }
3086 #endif
3087  if (j < m_kk-1) {
3088  for (k = j+1; k < m_kk; k++) {
3089  // an inner sum over all cations
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];
3093 #ifdef DEBUG_MODE
3094  if (m_debugCalc) {
3095  if (psi_ijk[n] != 0.0) {
3096  snj = speciesName(j) + "," + speciesName(k) + ":";
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]);
3099  }
3100  }
3101 #endif
3102  }
3103  }
3104  }
3105  }
3106 
3107  /*
3108  * For Anions, do the other anion interactions.
3109  */
3110  if (charge[j] < 0.0) {
3111  // sum over all anions
3112  if (j != i) {
3113  sum2 = sum2 + molality[j]*(2.0*Phi[counterIJ]);
3114 #ifdef DEBUG_MODE
3115  if (m_debugCalc) {
3116  if ((molality[j] * Phi[counterIJ])!= 0.0) {
3117  snj = speciesName(j) + ":";
3118  printf(" Phi term with %-12s 2 m_j Phi_aa = %10.5f\n", snj.c_str(),
3119  molality[j]*(2.0*Phi[counterIJ]));
3120  }
3121  }
3122 #endif
3123  }
3124  for (k = 1; k < m_kk; k++) {
3125  if (charge[k] > 0.0) {
3126  // two inner sums over cations
3127  n = k + j * m_kk + i * m_kk * m_kk;
3128  sum2 = sum2 + molality[j]*molality[k]*psi_ijk[n];
3129 #ifdef DEBUG_MODE
3130  if (m_debugCalc) {
3131  if (psi_ijk[n] != 0.0) {
3132  snj = speciesName(j) + "," + speciesName(k) + ":";
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]);
3135  }
3136  }
3137 #endif
3138  /*
3139  * Find the counterIJ for the symmetric binary interaction
3140  */
3141  n = m_kk*j + k;
3142  counterIJ2 = m_CounterIJ[n];
3143  sum4 = sum4 +
3144  (fabs(charge[i])*
3145  molality[j]*molality[k]*CMX[counterIJ2]);
3146 #ifdef DEBUG_MODE
3147  if (m_debugCalc) {
3148  if ((molality[j]*molality[k]*CMX[counterIJ2]) != 0.0) {
3149  snj = speciesName(j) + "," + speciesName(k) + ":";
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]);
3152  }
3153  }
3154 #endif
3155  }
3156  }
3157  }
3158 
3159  /*
3160  * for Anions, do the neutral species interaction
3161  */
3162  if (charge[j] == 0.0) {
3163  sum5 = sum5 + molality[j]*2.0*m_Lambda_nj(j,i);
3164 #ifdef DEBUG_MODE
3165  if (m_debugCalc) {
3166  if ((molality[j]*2.0*m_Lambda_nj(j,i)) != 0.0) {
3167  snj = speciesName(j) + ":";
3168  printf(" Lambda term with %-12s 2 m_j lam_ji = %10.5f\n", snj.c_str(),
3169  molality[j]*2.0*m_Lambda_nj(j,i));
3170  }
3171  }
3172 #endif
3173  /*
3174  * Zeta interaction term
3175  */
3176  for (size_t k = 1; k < m_kk; k++) {
3177  if (charge[k] > 0.0) {
3178  size_t izeta = j;
3179  size_t jzeta = k;
3180  size_t kzeta = i;
3181  n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
3182  double zeta = psi_ijk[n];
3183  if (zeta != 0.0) {
3184  sum5 = sum5 + molality[j]*molality[k]*zeta;
3185 #ifdef DEBUG_MODE
3186  if (m_debugCalc) {
3187  snj = speciesName(j) + "," + speciesName(k) + ":";
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]);
3190  }
3191 #endif
3192  }
3193  }
3194  }
3195  }
3196  }
3197  m_lnActCoeffMolal_Unscaled[i] = zsqF + sum1 + sum2 + sum3 + sum4 + sum5;
3198  gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]);
3199 #ifdef DEBUG_MODE
3200  if (m_debugCalc) {
3201  sni = speciesName(i);
3202  printf(" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f\n",
3203  sni.c_str(), m_lnActCoeffMolal_Unscaled[i], gamma_Unscaled[i]);
3204  }
3205 #endif
3206  }
3207  /*
3208  * ------ SUBSECTION FOR CALCULATING NEUTRAL SOLUTE ACT COEFF -------
3209  * ------ -> equations agree with my notes,
3210  * -> Equations agree with Pitzer,
3211  */
3212  if (charge[i] == 0.0) {
3213 #ifdef DEBUG_MODE
3214  if (m_debugCalc) {
3215  sni = speciesName(i);
3216  printf(" Contributions to ln(ActCoeff_%s):\n", sni.c_str());
3217  }
3218 #endif
3219  sum1 = 0.0;
3220  sum3 = 0.0;
3221  for (j = 1; j < m_kk; j++) {
3222  sum1 = sum1 + molality[j]*2.0*m_Lambda_nj(i,j);
3223 #ifdef DEBUG_MODE
3224  if (m_debugCalc) {
3225  if (m_Lambda_nj(i,j) != 0.0) {
3226  snj = speciesName(j) + ":";
3227  printf(" Lambda_n term on %-16s 2 m_j lambda_n_j = %10.5f\n", snj.c_str(),
3228  molality[j]*2.0*m_Lambda_nj(i,j));
3229  }
3230  }
3231 #endif
3232 
3233  /*
3234  * Zeta term -> we piggyback on the psi term
3235  */
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];
3241 #ifdef DEBUG_MODE
3242  if (m_debugCalc) {
3243  if (psi_ijk[n] != 0.0) {
3244  snj = speciesName(j) + "," + speciesName(k) + ":";
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]);
3247  }
3248  }
3249 #endif
3250  }
3251  }
3252  }
3253  }
3254  sum2 = 3.0 * molality[i]* molality[i] * m_Mu_nnn[i];
3255 #ifdef DEBUG_MODE
3256  if (m_debugCalc) {
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]);
3260  }
3261  }
3262 #endif
3263 
3264  m_lnActCoeffMolal_Unscaled[i] = sum1 + sum2 + sum3;
3265  gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]);
3266 #ifdef DEBUG_MODE
3267  if (m_debugCalc) {
3268  sni = speciesName(i);
3269  printf(" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f\n",
3270  sni.c_str(), m_lnActCoeffMolal_Unscaled[i], gamma_Unscaled[i]);
3271  }
3272 #endif
3273  }
3274 
3275  }
3276 #ifdef DEBUG_MODE
3277  if (m_debugCalc) {
3278  printf(" Step 9: \n");
3279  }
3280 #endif
3281  /*
3282  * -------- SUBSECTION FOR CALCULATING THE OSMOTIC COEFF ---------
3283  * -------- -> equations agree with my notes, Eqn. (117).
3284  * -> Equations agree with Pitzer, eqn.(62)
3285  */
3286  sum1 = 0.0;
3287  sum2 = 0.0;
3288  sum3 = 0.0;
3289  sum4 = 0.0;
3290  sum5 = 0.0;
3291  double sum6 = 0.0;
3292  double sum7 = 0.0;
3293  /*
3294  * term1 is the DH term in the osmotic coefficient expression
3295  * b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer
3296  * implementations.
3297  * Is = Ionic strength on the molality scale (units of (gmol/kg))
3298  * Aphi = A_Debye / 3 (units of sqrt(kg/gmol))
3299  */
3300  term1 = -Aphi * pow(Is,1.5) / (1.0 + 1.2 * sqrt(Is));
3301 
3302  for (j = 1; j < m_kk; j++) {
3303  /*
3304  * Loop Over Cations
3305  */
3306  if (charge[j] > 0.0) {
3307  for (k = 1; k < m_kk; k++) {
3308  if (charge[k] < 0.0) {
3309  /*
3310  * Find the counterIJ for the symmetric j,k binary interaction
3311  */
3312  n = m_kk*j + k;
3313  counterIJ = m_CounterIJ[n];
3314 
3315  sum1 = sum1 + molality[j]*molality[k]*
3316  (BphiMX[counterIJ] + molarcharge*CMX[counterIJ]);
3317  }
3318  }
3319 
3320  for (k = j+1; k < m_kk; k++) {
3321  if (j == (m_kk-1)) {
3322  // we should never reach this step
3323  printf("logic error 1 in Step 9 of hmw_act");
3324  exit(EXIT_FAILURE);
3325  }
3326  if (charge[k] > 0.0) {
3327  /*
3328  * Find the counterIJ for the symmetric j,k binary interaction
3329  * between 2 cations.
3330  */
3331  n = m_kk*j + k;
3332  counterIJ = m_CounterIJ[n];
3333  sum2 = sum2 + molality[j]*molality[k]*Phiphi[counterIJ];
3334  for (m = 1; m < m_kk; m++) {
3335  if (charge[m] < 0.0) {
3336  // species m is an anion
3337  n = m + k * m_kk + j * m_kk * m_kk;
3338  sum2 = sum2 +
3339  molality[j]*molality[k]*molality[m]*psi_ijk[n];
3340  }
3341  }
3342  }
3343  }
3344  }
3345 
3346  /*
3347  * Loop Over Anions
3348  */
3349  if (charge[j] < 0) {
3350  for (k = j+1; k < m_kk; k++) {
3351  if (j == m_kk-1) {
3352  // we should never reach this step
3353  printf("logic error 2 in Step 9 of hmw_act");
3354  exit(EXIT_FAILURE);
3355  }
3356  if (charge[k] < 0) {
3357  /*
3358  * Find the counterIJ for the symmetric j,k binary interaction
3359  * between two anions
3360  */
3361  n = m_kk*j + k;
3362  counterIJ = m_CounterIJ[n];
3363 
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;
3368  sum3 = sum3 +
3369  molality[j]*molality[k]*molality[m]*psi_ijk[n];
3370  }
3371  }
3372  }
3373  }
3374  }
3375 
3376  /*
3377  * Loop Over Neutral Species
3378  */
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);
3383  }
3384  if (charge[k] > 0.0) {
3385  sum5 = sum5 + molality[j]*molality[k]*m_Lambda_nj(j,k);
3386  }
3387  if (charge[k] == 0.0) {
3388  if (k > j) {
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);
3392  }
3393  }
3394  if (charge[k] < 0.0) {
3395  size_t izeta = j;
3396  for (m = 1; m < m_kk; m++) {
3397  if (charge[m] > 0.0) {
3398  size_t jzeta = m;
3399  n = k + jzeta * m_kk + izeta * m_kk * m_kk;
3400  double zeta = psi_ijk[n];
3401  if (zeta != 0.0) {
3402  sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta;
3403  }
3404  }
3405  }
3406  }
3407  }
3408  sum7 += molality[j]*molality[j]*molality[j]*m_Mu_nnn[j];
3409  }
3410  }
3411  sum_m_phi_minus_1 = 2.0 *
3412  (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
3413  /*
3414  * Calculate the osmotic coefficient from
3415  * osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i)
3416  */
3417  if (molalitysumUncropped > 1.0E-150) {
3418  osmotic_coef = 1.0 + (sum_m_phi_minus_1 / molalitysumUncropped);
3419  } else {
3420  osmotic_coef = 1.0;
3421  }
3422 #ifdef DEBUG_MODE
3423  if (printE) {
3424  printf("OsmCoef - 1 = %20.13g\n", osmotic_coef - 1.0);
3425  }
3426 #endif
3427 #ifdef DEBUG_MODE
3428  if (m_debugCalc) {
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);
3434  }
3435 
3436  if (m_debugCalc) {
3437  printf(" Step 10: \n");
3438  }
3439 #endif
3440  lnwateract = -(m_weightSolvent/1000.0) * molalitysumUncropped * osmotic_coef;
3441 
3442  /*
3443  * In Cantera, we define the activity coefficient of the solvent as
3444  *
3445  * act_0 = actcoeff_0 * Xmol_0
3446  *
3447  * We have just computed act_0. However, this routine returns
3448  * ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
3449  */
3450  double xmolSolvent = moleFraction(m_indexSolvent);
3451  double xx = std::max(m_xmolSolventMIN, xmolSolvent);
3452  m_lnActCoeffMolal_Unscaled[0] = lnwateract - log(xx);
3453 #ifdef DEBUG_MODE
3454  if (m_debugCalc) {
3455  double wateract = exp(lnwateract);
3456  printf(" Weight of Solvent = %16.7g\n", m_weightSolvent);
3457  printf(" molalitySumUncropped = %16.7g\n", molalitysumUncropped);
3458  printf(" ln_a_water=%10.6f a_water=%10.6f\n\n",
3459  lnwateract, wateract);
3460  }
3461 #endif
3462 }
3463 
3464 /**
3465  * s_update_dlnMolalityActCoeff_dT() (private, const )
3466  *
3467  * Using internally stored values, this function calculates
3468  * the temperature derivative of the logarithm of the
3469  * activity coefficient for all species in the mechanism.
3470  *
3471  * We assume that the activity coefficients are current.
3472  *
3473  * solvent activity coefficient is on the molality
3474  * scale. It's derivative is too.
3475  */
3476 void HMWSoln::s_update_dlnMolalityActCoeff_dT() const
3477 {
3478  /*
3479  * Zero the unscaled 2nd derivatives
3480  */
3481  m_dlnActCoeffMolaldT_Unscaled.assign(m_kk, 0.0);
3482  /*
3483  * Do the actual calculation of the unscaled temperature derivatives
3484  */
3485  s_updatePitzer_dlnMolalityActCoeff_dT();
3486 
3487  //double xmolSolvent = moleFraction(m_indexSolvent);
3488  //double xx = MAX(m_xmolSolventMIN, xmolSolvent);
3489  // double lnxs = log(xx);
3490 
3491  for (size_t k = 1; k < m_kk; k++) {
3492  if (CROP_speciesCropped_[k] == 2) {
3493  m_dlnActCoeffMolaldT_Unscaled[k] = 0.0;
3494  }
3495  }
3496 
3497  if (CROP_speciesCropped_[0]) {
3498  m_dlnActCoeffMolaldT_Unscaled[0] = 0.0;
3499  }
3500 
3501  /*
3502  * Do the pH scaling to the derivatives
3503  */
3504  s_updateScaling_pHScaling_dT();
3505 
3506 
3507 }
3508 
3509 /*************************************************************************************/
3510 
3511 /*
3512  * Calculate the Pitzer portion of the temperature
3513  * derivative of the log activity coefficients.
3514  * This is an internal routine.
3515  *
3516  * It may be assumed that the
3517  * Pitzer activity coefficient routine is called immediately
3518  * preceding the calling of this routine. Therefore, some
3519  * quantities do not need to be recalculated in this routine.
3520  *
3521  */
3522 void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
3523 {
3524 
3525  /*
3526  * HKM -> Assumption is made that the solvent is
3527  * species 0.
3528  */
3529 #ifdef DEBUG_MODE
3530  m_debugCalc = 0;
3531 #endif
3532  if (m_indexSolvent != 0) {
3533  printf("Wrong index solvent value!\n");
3534  exit(EXIT_FAILURE);
3535  }
3536 
3537  std::string sni, snj, snk;
3538 
3539  const double* molality = DATA_PTR(m_molalitiesCropped);
3540  const double* charge = DATA_PTR(m_speciesCharge);
3541  const double* beta0MX_L = DATA_PTR(m_Beta0MX_ij_L);
3542  const double* beta1MX_L = DATA_PTR(m_Beta1MX_ij_L);
3543  const double* beta2MX_L = DATA_PTR(m_Beta2MX_ij_L);
3544  const double* CphiMX_L = DATA_PTR(m_CphiMX_ij_L);
3545  const double* thetaij_L = DATA_PTR(m_Theta_ij_L);
3546  const double* alpha1MX = DATA_PTR(m_Alpha1MX_ij);
3547  const double* alpha2MX = DATA_PTR(m_Alpha2MX_ij);
3548  const double* psi_ijk_L = DATA_PTR(m_Psi_ijk_L);
3549  double* d_gamma_dT_Unscaled = DATA_PTR(m_gamma_tmp);
3550  /*
3551  * Local variables defined by Coltrin
3552  */
3553  double etheta[5][5], etheta_prime[5][5], sqrtIs;
3554  /*
3555  * Molality based ionic strength of the solution
3556  */
3557  double Is = 0.0;
3558  /*
3559  * Molarcharge of the solution: In Pitzer's notation,
3560  * this is his variable called "Z".
3561  */
3562  double molarcharge = 0.0;
3563  /*
3564  * molalitysum is the sum of the molalities over all solutes,
3565  * even those with zero charge.
3566  */
3567  double molalitysum = 0.0;
3568 
3569  double* gfunc = DATA_PTR(m_gfunc_IJ);
3570  double* g2func = DATA_PTR(m_g2func_IJ);
3571  double* hfunc = DATA_PTR(m_hfunc_IJ);
3572  double* h2func = DATA_PTR(m_h2func_IJ);
3573  double* BMX_L = DATA_PTR(m_BMX_IJ_L);
3574  double* BprimeMX_L= DATA_PTR(m_BprimeMX_IJ_L);
3575  double* BphiMX_L = DATA_PTR(m_BphiMX_IJ_L);
3576  double* Phi_L = DATA_PTR(m_Phi_IJ_L);
3577  double* Phiprime = DATA_PTR(m_Phiprime_IJ);
3578  double* Phiphi_L = DATA_PTR(m_PhiPhi_IJ_L);
3579  double* CMX_L = DATA_PTR(m_CMX_IJ_L);
3580 
3581  double x1, x2;
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;
3585 
3586  int z1, z2;
3587  size_t n, i, j, k, m, counterIJ, counterIJ2;
3588 
3589 #ifdef DEBUG_MODE
3590  if (m_debugCalc) {
3591  printf("\n Debugging information from "
3592  "s_Pitzer_dlnMolalityActCoeff_dT()\n");
3593  }
3594 #endif
3595  /*
3596  * Make sure the counter variables are setup
3597  */
3598  counterIJ_setup();
3599 
3600  /*
3601  * ---------- Calculate common sums over solutes ---------------------
3602  */
3603  for (n = 1; n < m_kk; n++) {
3604  // ionic strength
3605  Is += charge[n] * charge[n] * molality[n];
3606  // total molar charge
3607  molarcharge += fabs(charge[n]) * molality[n];
3608  molalitysum += molality[n];
3609  }
3610  Is *= 0.5;
3611 
3612  /*
3613  * Store the ionic molality in the object for reference.
3614  */
3615  m_IionicMolality = Is;
3616  sqrtIs = sqrt(Is);
3617 #ifdef DEBUG_MODE
3618  if (m_debugCalc) {
3619  printf(" Step 1: \n");
3620  printf(" ionic strenth = %14.7le \n total molar "
3621  "charge = %14.7le \n", Is, molarcharge);
3622  }
3623 #endif
3624 
3625  /*
3626  * The following call to calc_lambdas() calculates all 16 elements
3627  * of the elambda and elambda1 arrays, given the value of the
3628  * ionic strength (Is)
3629  */
3630  calc_lambdas(Is);
3631 
3632  /*
3633  * ----- Step 2: Find the coefficients E-theta and -------------------
3634  * E-thetaprime for all combinations of positive
3635  * unlike charges up to 4
3636  */
3637 #ifdef DEBUG_MODE
3638  if (m_debugCalc) {
3639  printf(" Step 2: \n");
3640  }
3641 #endif
3642  for (z1 = 1; z1 <=4; z1++) {
3643  for (z2 =1; z2 <=4; z2++) {
3644  calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
3645 #ifdef DEBUG_MODE
3646  if (m_debugCalc) {
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]);
3649  }
3650 #endif
3651  }
3652  }
3653 
3654 #ifdef DEBUG_MODE
3655  if (m_debugCalc) {
3656  printf(" Step 3: \n");
3657  printf(" Species Species g(x) "
3658  " hfunc(x) \n");
3659  }
3660 #endif
3661 
3662  /*
3663  *
3664  * calculate g(x) and hfunc(x) for each cation-anion pair MX
3665  * In the original literature, hfunc, was called gprime. However,
3666  * it's not the derivative of g(x), so I renamed it.
3667  */
3668  for (i = 1; i < (m_kk - 1); i++) {
3669  for (j = (i+1); j < m_kk; j++) {
3670  /*
3671  * Find the counterIJ for the symmetric binary interaction
3672  */
3673  n = m_kk*i + j;
3674  counterIJ = m_CounterIJ[n];
3675  /*
3676  * Only loop over oppositely charge species
3677  */
3678  if (charge[i]*charge[j] < 0) {
3679  /*
3680  * x is a reduced function variable
3681  */
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);
3687  } else {
3688  gfunc[counterIJ] = 0.0;
3689  hfunc[counterIJ] = 0.0;
3690  }
3691 
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);
3698  } else {
3699  g2func[counterIJ] = 0.0;
3700  h2func[counterIJ] = 0.0;
3701  }
3702  }
3703  } else {
3704  gfunc[counterIJ] = 0.0;
3705  hfunc[counterIJ] = 0.0;
3706  }
3707 #ifdef DEBUG_MODE
3708  if (m_debugCalc) {
3709  sni = speciesName(i);
3710  snj = speciesName(j);
3711  printf(" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(),
3712  gfunc[counterIJ], hfunc[counterIJ]);
3713  }
3714 #endif
3715  }
3716  }
3717 
3718  /*
3719  * ------- SUBSECTION TO CALCULATE BMX_L, BprimeMX_L, BphiMX_L ----------
3720  * ------- These are now temperature derivatives of the
3721  * previously calculated quantities.
3722  */
3723 #ifdef DEBUG_MODE
3724  if (m_debugCalc) {
3725  printf(" Step 4: \n");
3726  printf(" Species Species BMX "
3727  "BprimeMX BphiMX \n");
3728  }
3729 #endif
3730 
3731  for (i = 1; i < m_kk - 1; i++) {
3732  for (j = i+1; j < m_kk; j++) {
3733  /*
3734  * Find the counterIJ for the symmetric binary interaction
3735  */
3736  n = m_kk*i + j;
3737  counterIJ = m_CounterIJ[n];
3738  /*
3739  * both species have a non-zero charge, and one is positive
3740  * and the other is negative
3741  */
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];
3746 #ifdef DEBUG_MODE
3747  if (m_debugCalc) {
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]);
3751  }
3752 #endif
3753  if (Is > 1.0E-150) {
3754  BprimeMX_L[counterIJ] = (beta1MX_L[counterIJ] * hfunc[counterIJ]/Is +
3755  beta2MX_L[counterIJ] * h2func[counterIJ]/Is);
3756  } else {
3757  BprimeMX_L[counterIJ] = 0.0;
3758  }
3759  BphiMX_L[counterIJ] = BMX_L[counterIJ] + Is*BprimeMX_L[counterIJ];
3760  } else {
3761  BMX_L[counterIJ] = 0.0;
3762  BprimeMX_L[counterIJ] = 0.0;
3763  BphiMX_L[counterIJ] = 0.0;
3764  }
3765 #ifdef DEBUG_MODE
3766  if (m_debugCalc) {
3767  sni = speciesName(i);
3768  snj = speciesName(j);
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]);
3772  }
3773 #endif
3774  }
3775  }
3776 
3777  /*
3778  * --------- SUBSECTION TO CALCULATE CMX_L ----------
3779  * ---------
3780  */
3781 #ifdef DEBUG_MODE
3782  if (m_debugCalc) {
3783  printf(" Step 5: \n");
3784  printf(" Species Species CMX \n");
3785  }
3786 #endif
3787  for (i = 1; i < m_kk-1; i++) {
3788  for (j = i+1; j < m_kk; j++) {
3789  /*
3790  * Find the counterIJ for the symmetric binary interaction
3791  */
3792  n = m_kk*i + j;
3793  counterIJ = m_CounterIJ[n];
3794  /*
3795  * both species have a non-zero charge, and one is positive
3796  * and the other is negative
3797  */
3798  if (charge[i]*charge[j] < 0.0) {
3799  CMX_L[counterIJ] = CphiMX_L[counterIJ]/
3800  (2.0* sqrt(fabs(charge[i]*charge[j])));
3801  } else {
3802  CMX_L[counterIJ] = 0.0;
3803  }
3804 #ifdef DEBUG_MODE
3805  if (m_debugCalc) {
3806  sni = speciesName(i);
3807  snj = speciesName(j);
3808  printf(" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(),
3809  CMX_L[counterIJ]);
3810  }
3811 #endif
3812  }
3813  }
3814 
3815  /*
3816  * ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ----------
3817  * --------
3818  */
3819 #ifdef DEBUG_MODE
3820  if (m_debugCalc) {
3821  printf(" Step 6: \n");
3822  printf(" Species Species Phi_ij "
3823  " Phiprime_ij Phi^phi_ij \n");
3824  }
3825 #endif
3826  for (i = 1; i < m_kk-1; i++) {
3827  for (j = i+1; j < m_kk; j++) {
3828  /*
3829  * Find the counterIJ for the symmetric binary interaction
3830  */
3831  n = m_kk*i + j;
3832  counterIJ = m_CounterIJ[n];
3833  /*
3834  * both species have a non-zero charge, and one is positive
3835  * and the other is negative
3836  */
3837  if (charge[i]*charge[j] > 0) {
3838  z1 = (int) fabs(charge[i]);
3839  z2 = (int) fabs(charge[j]);
3840  //Phi[counterIJ] = thetaij_L[counterIJ] + etheta[z1][z2];
3841  Phi_L[counterIJ] = thetaij_L[counterIJ];
3842  //Phiprime[counterIJ] = etheta_prime[z1][z2];
3843  Phiprime[counterIJ] = 0.0;
3844  Phiphi_L[counterIJ] = Phi_L[counterIJ] + Is * Phiprime[counterIJ];
3845  } else {
3846  Phi_L[counterIJ] = 0.0;
3847  Phiprime[counterIJ] = 0.0;
3848  Phiphi_L[counterIJ] = 0.0;
3849  }
3850 #ifdef DEBUG_MODE
3851  if (m_debugCalc) {
3852  sni = speciesName(i);
3853  snj = speciesName(j);
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]);
3857  }
3858 #endif
3859  }
3860  }
3861 
3862  /*
3863  * ----------- SUBSECTION FOR CALCULATION OF dFdT ---------------------
3864  */
3865 #ifdef DEBUG_MODE
3866  if (m_debugCalc) {
3867  printf(" Step 7: \n");
3868  }
3869 #endif
3870  // A_Debye_Huckel = 0.5092; (units = sqrt(kg/gmol))
3871  // A_Debye_Huckel = 0.5107; <- This value is used to match GWB data
3872  // ( A * ln(10) = 1.17593)
3873  // Aphi = A_Debye_Huckel * 2.30258509 / 3.0;
3874 
3875  double dA_DebyedT = dA_DebyedT_TP();
3876  double dAphidT = dA_DebyedT /3.0;
3877 #ifdef DEBUG_HKM
3878  //dAphidT = 0.0;
3879 #endif
3880  //F = -Aphi * ( sqrt(Is) / (1.0 + 1.2*sqrt(Is))
3881  // + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
3882  //dAphidT = Al / (4.0 * GasConstant * T * T);
3883  dFdT = -dAphidT * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
3884  + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
3885 #ifdef DEBUG_MODE
3886  if (m_debugCalc) {
3887  printf(" initial value of dFdT = %10.6f \n", dFdT);
3888  }
3889 #endif
3890  for (i = 1; i < m_kk-1; i++) {
3891  for (j = i+1; j < m_kk; j++) {
3892  /*
3893  * Find the counterIJ for the symmetric binary interaction
3894  */
3895  n = m_kk*i + j;
3896  counterIJ = m_CounterIJ[n];
3897  /*
3898  * both species have a non-zero charge, and one is positive
3899  * and the other is negative
3900  */
3901  if (charge[i]*charge[j] < 0) {
3902  dFdT = dFdT + molality[i]*molality[j] * BprimeMX_L[counterIJ];
3903  }
3904  /*
3905  * Both species have a non-zero charge, and they
3906  * have the same sign, e.g., both positive or both negative.
3907  */
3908  if (charge[i]*charge[j] > 0) {
3909  dFdT = dFdT + molality[i]*molality[j] * Phiprime[counterIJ];
3910  }
3911 #ifdef DEBUG_MODE
3912  if (m_debugCalc) {
3913  printf(" dFdT = %10.6f \n", dFdT);
3914  }
3915 #endif
3916  }
3917  }
3918 #ifdef DEBUG_MODE
3919  if (m_debugCalc) {
3920  printf(" Step 8: \n");
3921  }
3922 #endif
3923 
3924  for (i = 1; i < m_kk; i++) {
3925 
3926  /*
3927  * -------- SUBSECTION FOR CALCULATING THE dACTCOEFFdT FOR CATIONS -----
3928  * --
3929  */
3930  if (charge[i] > 0) {
3931  // species i is the cation (positive) to calc the actcoeff
3932  zsqdFdT = charge[i]*charge[i]*dFdT;
3933  sum1 = 0.0;
3934  sum2 = 0.0;
3935  sum3 = 0.0;
3936  sum4 = 0.0;
3937  sum5 = 0.0;
3938  for (j = 1; j < m_kk; j++) {
3939  /*
3940  * Find the counterIJ for the symmetric binary interaction
3941  */
3942  n = m_kk*i + j;
3943  counterIJ = m_CounterIJ[n];
3944 
3945  if (charge[j] < 0.0) {
3946  // sum over all anions
3947  sum1 = sum1 + molality[j]*
3948  (2.0*BMX_L[counterIJ] + molarcharge*CMX_L[counterIJ]);
3949  if (j < m_kk-1) {
3950  /*
3951  * This term is the ternary interaction involving the
3952  * non-duplicate sum over double anions, j, k, with
3953  * respect to the cation, i.
3954  */
3955  for (k = j+1; k < m_kk; k++) {
3956  // an inner sum over all anions
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];
3960  }
3961  }
3962  }
3963  }
3964 
3965 
3966  if (charge[j] > 0.0) {
3967  // sum over all cations
3968  if (j != i) {
3969  sum2 = sum2 + molality[j]*(2.0*Phi_L[counterIJ]);
3970  }
3971  for (k = 1; k < m_kk; k++) {
3972  if (charge[k] < 0.0) {
3973  // two inner sums over anions
3974 
3975  n = k + j * m_kk + i * m_kk * m_kk;
3976  sum2 = sum2 + molality[j]*molality[k]*psi_ijk_L[n];
3977  /*
3978  * Find the counterIJ for the j,k interaction
3979  */
3980  n = m_kk*j + k;
3981  counterIJ2 = m_CounterIJ[n];
3982  sum4 = sum4 + (fabs(charge[i])*
3983  molality[j]*molality[k]*CMX_L[counterIJ2]);
3984  }
3985  }
3986  }
3987 
3988  /*
3989  * Handle neutral j species
3990  */
3991  if (charge[j] == 0) {
3992  sum5 = sum5 + molality[j]*2.0*m_Lambda_nj_L(j,i);
3993  }
3994  /*
3995  * Zeta interaction term
3996  */
3997  for (size_t k = 1; k < m_kk; k++) {
3998  if (charge[k] < 0.0) {
3999  size_t izeta = j;
4000  size_t jzeta = i;
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;
4005  }
4006  }
4007  }
4008  }
4009  /*
4010  * Add all of the contributions up to yield the log of the
4011  * solute activity coefficients (molality scale)
4012  */
4013  m_dlnActCoeffMolaldT_Unscaled[i] =
4014  zsqdFdT + sum1 + sum2 + sum3 + sum4 + sum5;
4015  d_gamma_dT_Unscaled[i] = exp(m_dlnActCoeffMolaldT_Unscaled[i]);
4016 #ifdef DEBUG_MODE
4017  if (m_debugCalc) {
4018  sni = speciesName(i);
4019  printf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f \n",
4020  sni.c_str(), m_dlnActCoeffMolaldT_Unscaled[i], d_gamma_dT_Unscaled[i]);
4021  printf(" %12g %12g %12g %12g %12g %12g\n",
4022  zsqdFdT, sum1, sum2, sum3, sum4, sum5);
4023  }
4024 #endif
4025  }
4026 
4027  /*
4028  * ------ SUBSECTION FOR CALCULATING THE dACTCOEFFdT FOR ANIONS ------
4029  *
4030  */
4031  if (charge[i] < 0) {
4032  // species i is an anion (negative)
4033  zsqdFdT = charge[i]*charge[i]*dFdT;
4034  sum1 = 0.0;
4035  sum2 = 0.0;
4036  sum3 = 0.0;
4037  sum4 = 0.0;
4038  sum5 = 0.0;
4039  for (j = 1; j < m_kk; j++) {
4040  /*
4041  * Find the counterIJ for the symmetric binary interaction
4042  */
4043  n = m_kk*i + j;
4044  counterIJ = m_CounterIJ[n];
4045 
4046  /*
4047  * For Anions, do the cation interactions.
4048  */
4049  if (charge[j] > 0) {
4050  sum1 = sum1 + molality[j]*
4051  (2.0*BMX_L[counterIJ] + molarcharge*CMX_L[counterIJ]);
4052  if (j < m_kk-1) {
4053  for (k = j+1; k < m_kk; k++) {
4054  // an inner sum over all cations
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];
4058  }
4059  }
4060  }
4061  }
4062 
4063  /*
4064  * For Anions, do the other anion interactions.
4065  */
4066  if (charge[j] < 0.0) {
4067  // sum over all anions
4068  if (j != i) {
4069  sum2 = sum2 + molality[j]*(2.0*Phi_L[counterIJ]);
4070  }
4071  for (k = 1; k < m_kk; k++) {
4072  if (charge[k] > 0.0) {
4073  // two inner sums over cations
4074  n = k + j * m_kk + i * m_kk * m_kk;
4075  sum2 = sum2 + molality[j]*molality[k]*psi_ijk_L[n];
4076  /*
4077  * Find the counterIJ for the symmetric binary interaction
4078  */
4079  n = m_kk*j + k;
4080  counterIJ2 = m_CounterIJ[n];
4081  sum4 = sum4 +
4082  (fabs(charge[i])*
4083  molality[j]*molality[k]*CMX_L[counterIJ2]);
4084  }
4085  }
4086  }
4087 
4088  /*
4089  * for Anions, do the neutral species interaction
4090  */
4091  if (charge[j] == 0.0) {
4092  sum5 = sum5 + molality[j]*2.0*m_Lambda_nj_L(j,i);
4093  for (size_t k = 1; k < m_kk; k++) {
4094  if (charge[k] > 0.0) {
4095  size_t izeta = j;
4096  size_t jzeta = k;
4097  size_t kzeta = i;
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;
4102  }
4103  }
4104  }
4105  }
4106  }
4107  m_dlnActCoeffMolaldT_Unscaled[i] =
4108  zsqdFdT + sum1 + sum2 + sum3 + sum4 + sum5;
4109  d_gamma_dT_Unscaled[i] = exp(m_dlnActCoeffMolaldT_Unscaled[i]);
4110 #ifdef DEBUG_MODE
4111  if (m_debugCalc) {
4112  sni = speciesName(i);
4113  printf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f\n",
4114  sni.c_str(), m_dlnActCoeffMolaldT_Unscaled[i], d_gamma_dT_Unscaled[i]);
4115  printf(" %12g %12g %12g %12g %12g %12g\n",
4116  zsqdFdT, sum1, sum2, sum3, sum4, sum5);
4117  }
4118 #endif
4119  }
4120  /*
4121  * ------ SUBSECTION FOR CALCULATING NEUTRAL SOLUTE ACT COEFF -------
4122  * ------ -> equations agree with my notes,
4123  * -> Equations agree with Pitzer,
4124  */
4125  if (charge[i] == 0.0) {
4126  sum1 = 0.0;
4127  sum3 = 0.0;
4128  for (j = 1; j < m_kk; j++) {
4129  sum1 = sum1 + molality[j]*2.0*m_Lambda_nj_L(i,j);
4130  /*
4131  * Zeta term -> we piggyback on the psi term
4132  */
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];
4138  }
4139  }
4140  }
4141  }
4142  sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_L[i];
4143  m_dlnActCoeffMolaldT_Unscaled[i] = sum1 + sum2 + sum3;
4144  d_gamma_dT_Unscaled[i] = exp(m_dlnActCoeffMolaldT_Unscaled[i]);
4145 #ifdef DEBUG_MODE
4146  if (m_debugCalc) {
4147  sni = speciesName(i);
4148  printf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f \n",
4149  sni.c_str(), m_dlnActCoeffMolaldT_Unscaled[i], d_gamma_dT_Unscaled[i]);
4150  }
4151 #endif
4152  }
4153 
4154  }
4155 #ifdef DEBUG_MODE
4156  if (m_debugCalc) {
4157  printf(" Step 9: \n");
4158  }
4159 #endif
4160  /*
4161  * ------ SUBSECTION FOR CALCULATING THE d OSMOTIC COEFF dT ---------
4162  *
4163  */
4164  sum1 = 0.0;
4165  sum2 = 0.0;
4166  sum3 = 0.0;
4167  sum4 = 0.0;
4168  sum5 = 0.0;
4169  double sum6 = 0.0;
4170  double sum7 = 0.0;
4171  /*
4172  * term1 is the temperature derivative of the
4173  * DH term in the osmotic coefficient expression
4174  * b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer
4175  * implementations.
4176  * Is = Ionic strength on the molality scale (units of (gmol/kg))
4177  * Aphi = A_Debye / 3 (units of sqrt(kg/gmol))
4178  */
4179  term1 = -dAphidT * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
4180 
4181  for (j = 1; j < m_kk; j++) {
4182  /*
4183  * Loop Over Cations
4184  */
4185  if (charge[j] > 0.0) {
4186  for (k = 1; k < m_kk; k++) {
4187  if (charge[k] < 0.0) {
4188  /*
4189  * Find the counterIJ for the symmetric j,k binary interaction
4190  */
4191  n = m_kk*j + k;
4192  counterIJ = m_CounterIJ[n];
4193 
4194  sum1 = sum1 + molality[j]*molality[k]*
4195  (BphiMX_L[counterIJ] + molarcharge*CMX_L[counterIJ]);
4196  }
4197  }
4198 
4199  for (k = j+1; k < m_kk; k++) {
4200  if (j == (m_kk-1)) {
4201  // we should never reach this step
4202  printf("logic error 1 in Step 9 of hmw_act");
4203  exit(EXIT_FAILURE);
4204  }
4205  if (charge[k] > 0.0) {
4206  /*
4207  * Find the counterIJ for the symmetric j,k binary interaction
4208  * between 2 cations.
4209  */
4210  n = m_kk*j + k;
4211  counterIJ = m_CounterIJ[n];
4212  sum2 = sum2 + molality[j]*molality[k]*Phiphi_L[counterIJ];
4213  for (m = 1; m < m_kk; m++) {
4214  if (charge[m] < 0.0) {
4215  // species m is an anion
4216  n = m + k * m_kk + j * m_kk * m_kk;
4217  sum2 = sum2 +
4218  molality[j]*molality[k]*molality[m]*psi_ijk_L[n];
4219  }
4220  }
4221  }
4222  }
4223  }
4224 
4225  /*
4226  * Loop Over Anions
4227  */
4228  if (charge[j] < 0) {
4229  for (k = j+1; k < m_kk; k++) {
4230  if (j == m_kk-1) {
4231  // we should never reach this step
4232  printf("logic error 2 in Step 9 of hmw_act");
4233  exit(EXIT_FAILURE);
4234  }
4235  if (charge[k] < 0) {
4236  /*
4237  * Find the counterIJ for the symmetric j,k binary interaction
4238  * between two anions
4239  */
4240  n = m_kk*j + k;
4241  counterIJ = m_CounterIJ[n];
4242 
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;
4247  sum3 = sum3 +
4248  molality[j]*molality[k]*molality[m]*psi_ijk_L[n];
4249  }
4250  }
4251  }
4252  }
4253  }
4254 
4255  /*
4256  * Loop Over Neutral Species
4257  */
4258  if (charge[j] == 0) {
4259  for (k = 1; k < m_kk; k++) {
4260  if (charge[k] < 0.0) {
4261  sum4 = sum4 + molality[j]*molality[k]*m_Lambda_nj_L(j,k);
4262  }
4263  if (charge[k] > 0.0) {
4264  sum5 = sum5 + molality[j]*molality[k]*m_Lambda_nj_L(j,k);
4265  }
4266  if (charge[k] == 0.0) {
4267  if (k > j) {
4268  sum6 = sum6 + molality[j]*molality[k]*m_Lambda_nj_L(j,k);
4269  } else if (k == j) {
4270  sum6 = sum6 + 0.5 * molality[j]*molality[k]*m_Lambda_nj_L(j,k);
4271  }
4272  }
4273  if (charge[k] < 0.0) {
4274  size_t izeta = j;
4275  for (m = 1; m < m_kk; m++) {
4276  if (charge[m] > 0.0) {
4277  size_t jzeta = m;
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;
4282  }
4283  }
4284  }
4285  }
4286  }
4287  sum7 += molality[j]*molality[j]*molality[j]*m_Mu_nnn_L[j];
4288  }
4289  }
4290  sum_m_phi_minus_1 = 2.0 *
4291  (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
4292  /*
4293  * Calculate the osmotic coefficient from
4294  * osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i)
4295  */
4296  if (molalitysum > 1.0E-150) {
4297  d_osmotic_coef_dT = 0.0 + (sum_m_phi_minus_1 / molalitysum);
4298  } else {
4299  d_osmotic_coef_dT = 0.0;
4300  }
4301 #ifdef DEBUG_MODE
4302  if (m_debugCalc) {
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);
4308  }
4309 
4310  if (m_debugCalc) {
4311  printf(" Step 10: \n");
4312  }
4313 #endif
4314  d_lnwateract_dT = -(m_weightSolvent/1000.0) * molalitysum * d_osmotic_coef_dT;
4315 
4316  /*
4317  * In Cantera, we define the activity coefficient of the solvent as
4318  *
4319  * act_0 = actcoeff_0 * Xmol_0
4320  *
4321  * We have just computed act_0. However, this routine returns
4322  * ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
4323  */
4324  //double xmolSolvent = moleFraction(m_indexSolvent);
4325  m_dlnActCoeffMolaldT_Unscaled[0] = d_lnwateract_dT;
4326 #ifdef DEBUG_MODE
4327  if (m_debugCalc) {
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);
4331  }
4332 #endif
4333 }
4334 
4335 /**
4336  * This function calculates the temperature second derivative
4337  * of the natural logarithm of the molality activity
4338  * coefficients.
4339  */
4340 void HMWSoln::s_update_d2lnMolalityActCoeff_dT2() const
4341 {
4342  /*
4343  * Zero the unscaled 2nd derivatives
4344  */
4345  m_d2lnActCoeffMolaldT2_Unscaled.assign(m_kk, 0.0);
4346  /*
4347  * Calculate the unscaled 2nd derivatives
4348  */
4349  s_updatePitzer_d2lnMolalityActCoeff_dT2();
4350 
4351  //double xmolSolvent = moleFraction(m_indexSolvent);
4352  //double xx = MAX(m_xmolSolventMIN, xmolSolvent);
4353  //double lnxs = log(xx);
4354 
4355  for (size_t k = 1; k < m_kk; k++) {
4356  if (CROP_speciesCropped_[k] == 2) {
4357  m_d2lnActCoeffMolaldT2_Unscaled[k] = 0.0;
4358  }
4359  }
4360 
4361  if (CROP_speciesCropped_[0]) {
4362  m_d2lnActCoeffMolaldT2_Unscaled[0] = 0.0;
4363  }
4364 
4365  /*
4366  * Scale the 2nd derivatives
4367  */
4368  s_updateScaling_pHScaling_dT2();
4369 }
4370 
4371 /*************************************************************************************/
4372 
4373 /*
4374  * s_updatePitzer_d2lnMolalityActCoeff_dT2() (private, const )
4375  *
4376  * Using internally stored values, this function calculates
4377  * the temperature 2nd derivative of the logarithm of the
4378  * activity coefficient for all species in the mechanism.
4379  * This is an internal routine
4380  *
4381  * We assume that the activity coefficients and first temperature
4382  * derivatives of the activity coefficients are current.
4383  *
4384  * It may be assumed that the
4385  * Pitzer activity coefficient and first deriv routine are called immediately
4386  * preceding the calling of this routine. Therefore, some
4387  * quantities do not need to be recalculated in this routine.
4388  *
4389  * solvent activity coefficient is on the molality
4390  * scale. It's derivatives are too.
4391  */
4392 void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
4393 {
4394 
4395  /*
4396  * HKM -> Assumption is made that the solvent is
4397  * species 0.
4398  */
4399 #ifdef DEBUG_MODE
4400  m_debugCalc = 0;
4401 #endif
4402  if (m_indexSolvent != 0) {
4403  printf("Wrong index solvent value!\n");
4404  exit(EXIT_FAILURE);
4405  }
4406 
4407  std::string sni, snj, snk;
4408 
4409  const double* molality = DATA_PTR(m_molalitiesCropped);
4410  const double* charge = DATA_PTR(m_speciesCharge);
4411  const double* beta0MX_LL= DATA_PTR(m_Beta0MX_ij_LL);
4412  const double* beta1MX_LL= DATA_PTR(m_Beta1MX_ij_LL);
4413  const double* beta2MX_LL= DATA_PTR(m_Beta2MX_ij_LL);
4414  const double* CphiMX_LL = DATA_PTR(m_CphiMX_ij_LL);
4415  const double* thetaij_LL= DATA_PTR(m_Theta_ij_LL);
4416  const double* alpha1MX = DATA_PTR(m_Alpha1MX_ij);
4417  const double* alpha2MX = DATA_PTR(m_Alpha2MX_ij);
4418  const double* psi_ijk_LL= DATA_PTR(m_Psi_ijk_LL);
4419 
4420  /*
4421  * Local variables defined by Coltrin
4422  */
4423  double etheta[5][5], etheta_prime[5][5], sqrtIs;
4424  /*
4425  * Molality based ionic strength of the solution
4426  */
4427  double Is = 0.0;
4428  /*
4429  * Molarcharge of the solution: In Pitzer's notation,
4430  * this is his variable called "Z".
4431  */
4432  double molarcharge = 0.0;
4433  /*
4434  * molalitysum is the sum of the molalities over all solutes,
4435  * even those with zero charge.
4436  */
4437  double molalitysum = 0.0;
4438 
4439  double* gfunc = DATA_PTR(m_gfunc_IJ);
4440  double* g2func = DATA_PTR(m_g2func_IJ);
4441  double* hfunc = DATA_PTR(m_hfunc_IJ);
4442  double* h2func = DATA_PTR(m_h2func_IJ);
4443  double* BMX_LL = DATA_PTR(m_BMX_IJ_LL);
4444  double* BprimeMX_LL=DATA_PTR(m_BprimeMX_IJ_LL);
4445  double* BphiMX_LL= DATA_PTR(m_BphiMX_IJ_LL);
4446  double* Phi_LL = DATA_PTR(m_Phi_IJ_LL);
4447  double* Phiprime = DATA_PTR(m_Phiprime_IJ);
4448  double* Phiphi_LL= DATA_PTR(m_PhiPhi_IJ_LL);
4449  double* CMX_LL = DATA_PTR(m_CMX_IJ_LL);
4450 
4451 
4452  double x1, x2;
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;
4456 
4457  int z1, z2;
4458  size_t n, i, j, k, m, counterIJ, counterIJ2;
4459 
4460 #ifdef DEBUG_MODE
4461  if (m_debugCalc) {
4462  printf("\n Debugging information from "
4463  "s_Pitzer_d2lnMolalityActCoeff_dT2()\n");
4464  }
4465 #endif
4466  /*
4467  * Make sure the counter variables are setup
4468  */
4469  counterIJ_setup();
4470 
4471 
4472  /*
4473  * ---------- Calculate common sums over solutes ---------------------
4474  */
4475  for (n = 1; n < m_kk; n++) {
4476  // ionic strength
4477  Is += charge[n] * charge[n] * molality[n];
4478  // total molar charge
4479  molarcharge += fabs(charge[n]) * molality[n];
4480  molalitysum += molality[n];
4481  }
4482  Is *= 0.5;
4483 
4484  /*
4485  * Store the ionic molality in the object for reference.
4486  */
4487  m_IionicMolality = Is;
4488  sqrtIs = sqrt(Is);
4489 #ifdef DEBUG_MODE
4490  if (m_debugCalc) {
4491  printf(" Step 1: \n");
4492  printf(" ionic strenth = %14.7le \n total molar "
4493  "charge = %14.7le \n", Is, molarcharge);
4494  }
4495 #endif
4496 
4497  /*
4498  * The following call to calc_lambdas() calculates all 16 elements
4499  * of the elambda and elambda1 arrays, given the value of the
4500  * ionic strength (Is)
4501  */
4502  calc_lambdas(Is);
4503 
4504  /*
4505  * ----- Step 2: Find the coefficients E-theta and -------------------
4506  * E-thetaprime for all combinations of positive
4507  * unlike charges up to 4
4508  */
4509 #ifdef DEBUG_MODE
4510  if (m_debugCalc) {
4511  printf(" Step 2: \n");
4512  }
4513 #endif
4514  for (z1 = 1; z1 <=4; z1++) {
4515  for (z2 =1; z2 <=4; z2++) {
4516  calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
4517 #ifdef DEBUG_MODE
4518  if (m_debugCalc) {
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]);
4521  }
4522 #endif
4523  }
4524  }
4525 
4526 #ifdef DEBUG_MODE
4527  if (m_debugCalc) {
4528  printf(" Step 3: \n");
4529  printf(" Species Species g(x) "
4530  " hfunc(x) \n");
4531  }
4532 #endif
4533 
4534  /*
4535  *
4536  * calculate gfunc(x) and hfunc(x) for each cation-anion pair MX
4537  * In the original literature, hfunc, was called gprime. However,
4538  * it's not the derivative of gfunc(x), so I renamed it.
4539  */
4540  for (i = 1; i < (m_kk - 1); i++) {
4541  for (j = (i+1); j < m_kk; j++) {
4542  /*
4543  * Find the counterIJ for the symmetric binary interaction
4544  */
4545  n = m_kk*i + j;
4546  counterIJ = m_CounterIJ[n];
4547  /*
4548  * Only loop over oppositely charge species
4549  */
4550  if (charge[i]*charge[j] < 0) {
4551  /*
4552  * x is a reduced function variable
4553  */
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);
4559  } else {
4560  gfunc[counterIJ] = 0.0;
4561  hfunc[counterIJ] = 0.0;
4562  }
4563 
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);
4570  } else {
4571  g2func[counterIJ] = 0.0;
4572  h2func[counterIJ] = 0.0;
4573  }
4574  }
4575  } else {
4576  gfunc[counterIJ] = 0.0;
4577  hfunc[counterIJ] = 0.0;
4578  }
4579 #ifdef DEBUG_MODE
4580  if (m_debugCalc) {
4581  sni = speciesName(i);
4582  snj = speciesName(j);
4583  printf(" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(),
4584  gfunc[counterIJ], hfunc[counterIJ]);
4585  }
4586 #endif
4587  }
4588  }
4589  /*
4590  * ------- SUBSECTION TO CALCULATE BMX_L, BprimeMX_LL, BphiMX_L ----------
4591  * ------- These are now temperature derivatives of the
4592  * previously calculated quantities.
4593  */
4594 #ifdef DEBUG_MODE
4595  if (m_debugCalc) {
4596  printf(" Step 4: \n");
4597  printf(" Species Species BMX "
4598  "BprimeMX BphiMX \n");
4599  }
4600 #endif
4601 
4602  for (i = 1; i < m_kk - 1; i++) {
4603  for (j = i+1; j < m_kk; j++) {
4604  /*
4605  * Find the counterIJ for the symmetric binary interaction
4606  */
4607  n = m_kk*i + j;
4608  counterIJ = m_CounterIJ[n];
4609  /*
4610  * both species have a non-zero charge, and one is positive
4611  * and the other is negative
4612  */
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];
4617 #ifdef DEBUG_MODE
4618  if (m_debugCalc) {
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]);
4622  }
4623 #endif
4624  if (Is > 1.0E-150) {
4625  BprimeMX_LL[counterIJ] = (beta1MX_LL[counterIJ] * hfunc[counterIJ]/Is +
4626  beta2MX_LL[counterIJ] * h2func[counterIJ]/Is);
4627  } else {
4628  BprimeMX_LL[counterIJ] = 0.0;
4629  }
4630  BphiMX_LL[counterIJ] = BMX_LL[counterIJ] + Is*BprimeMX_LL[counterIJ];
4631  } else {
4632  BMX_LL[counterIJ] = 0.0;
4633  BprimeMX_LL[counterIJ] = 0.0;
4634  BphiMX_LL[counterIJ] = 0.0;
4635  }
4636 #ifdef DEBUG_MODE
4637  if (m_debugCalc) {
4638  sni = speciesName(i);
4639  snj = speciesName(j);
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]);
4643  }
4644 #endif
4645  }
4646  }
4647 
4648  /*
4649  * --------- SUBSECTION TO CALCULATE CMX_LL ----------
4650  * ---------
4651  */
4652 #ifdef DEBUG_MODE
4653  if (m_debugCalc) {
4654  printf(" Step 5: \n");
4655  printf(" Species Species CMX \n");
4656  }
4657 #endif
4658  for (i = 1; i < m_kk-1; i++) {
4659  for (j = i+1; j < m_kk; j++) {
4660  /*
4661  * Find the counterIJ for the symmetric binary interaction
4662  */
4663  n = m_kk*i + j;
4664  counterIJ = m_CounterIJ[n];
4665  /*
4666  * both species have a non-zero charge, and one is positive
4667  * and the other is negative
4668  */
4669  if (charge[i]*charge[j] < 0.0) {
4670  CMX_LL[counterIJ] = CphiMX_LL[counterIJ]/
4671  (2.0* sqrt(fabs(charge[i]*charge[j])));
4672  } else {
4673  CMX_LL[counterIJ] = 0.0;
4674  }
4675 #ifdef DEBUG_MODE
4676  if (m_debugCalc) {
4677  sni = speciesName(i);
4678  snj = speciesName(j);
4679  printf(" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(),
4680  CMX_LL[counterIJ]);
4681  }
4682 #endif
4683  }
4684  }
4685 
4686  /*
4687  * ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ----------
4688  * --------
4689  */
4690 #ifdef DEBUG_MODE
4691  if (m_debugCalc) {
4692  printf(" Step 6: \n");
4693  printf(" Species Species Phi_ij "
4694  " Phiprime_ij Phi^phi_ij \n");
4695  }
4696 #endif
4697  for (i = 1; i < m_kk-1; i++) {
4698  for (j = i+1; j < m_kk; j++) {
4699  /*
4700  * Find the counterIJ for the symmetric binary interaction
4701  */
4702  n = m_kk*i + j;
4703  counterIJ = m_CounterIJ[n];
4704  /*
4705  * both species have a non-zero charge, and one is positive
4706  * and the other is negative
4707  */
4708  if (charge[i]*charge[j] > 0) {
4709  z1 = (int) fabs(charge[i]);
4710  z2 = (int) fabs(charge[j]);
4711  //Phi[counterIJ] = thetaij[counterIJ] + etheta[z1][z2];
4712  //Phi_L[counterIJ] = thetaij_L[counterIJ];
4713  Phi_LL[counterIJ] = thetaij_LL[counterIJ];
4714  //Phiprime[counterIJ] = etheta_prime[z1][z2];
4715  Phiprime[counterIJ] = 0.0;
4716  //Phiphi[counterIJ] = Phi[counterIJ] + Is * Phiprime[counterIJ];
4717  //Phiphi_L[counterIJ] = Phi_L[counterIJ] + Is * Phiprime[counterIJ];
4718  Phiphi_LL[counterIJ] = Phi_LL[counterIJ];
4719  } else {
4720  Phi_LL[counterIJ] = 0.0;
4721  Phiprime[counterIJ] = 0.0;
4722  Phiphi_LL[counterIJ] = 0.0;
4723  }
4724 #ifdef DEBUG_MODE
4725  if (m_debugCalc) {
4726  sni = speciesName(i);
4727  snj = speciesName(j);
4728  //printf(" %-16s %-16s %10.6f %10.6f %10.6f \n",
4729  // sni.c_str(), snj.c_str(),
4730  // Phi_L[counterIJ], Phiprime[counterIJ], Phiphi_L[counterIJ] );
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]);
4734  }
4735 #endif
4736  }
4737  }
4738 
4739  /*
4740  * ----------- SUBSECTION FOR CALCULATION OF d2FdT2 ---------------------
4741  */
4742 #ifdef DEBUG_MODE
4743  if (m_debugCalc) {
4744  printf(" Step 7: \n");
4745  }
4746 #endif
4747  // A_Debye_Huckel = 0.5092; (units = sqrt(kg/gmol))
4748  // A_Debye_Huckel = 0.5107; <- This value is used to match GWB data
4749  // ( A * ln(10) = 1.17593)
4750  // Aphi = A_Debye_Huckel * 2.30258509 / 3.0;
4751  // Aphi = m_A_Debye / 3.0;
4752 
4753  //double dA_DebyedT = dA_DebyedT_TP();
4754  //double dAphidT = dA_DebyedT /3.0;
4755  double d2AphidT2 = d2A_DebyedT2_TP() / 3.0;
4756 #ifdef DEBUG_HKM
4757  //d2AphidT2 = 0.0;
4758 #endif
4759  //F = -Aphi * ( sqrt(Is) / (1.0 + 1.2*sqrt(Is))
4760  // + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
4761  //dAphidT = Al / (4.0 * GasConstant * T * T);
4762  //dFdT = -dAphidT * ( sqrt(Is) / (1.0 + 1.2*sqrt(Is))
4763  // + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
4764  d2FdT2 = -d2AphidT2 * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
4765  + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
4766 #ifdef DEBUG_MODE
4767  if (m_debugCalc) {
4768  printf(" initial value of d2FdT2 = %10.6f \n", d2FdT2);
4769  }
4770 #endif
4771  for (i = 1; i < m_kk-1; i++) {
4772  for (j = i+1; j < m_kk; j++) {
4773  /*
4774  * Find the counterIJ for the symmetric binary interaction
4775  */
4776  n = m_kk*i + j;
4777  counterIJ = m_CounterIJ[n];
4778  /*
4779  * both species have a non-zero charge, and one is positive
4780  * and the other is negative
4781  */
4782  if (charge[i]*charge[j] < 0) {
4783  d2FdT2 = d2FdT2 + molality[i]*molality[j] * BprimeMX_LL[counterIJ];
4784  }
4785  /*
4786  * Both species have a non-zero charge, and they
4787  * have the same sign, e.g., both positive or both negative.
4788  */
4789  if (charge[i]*charge[j] > 0) {
4790  d2FdT2 = d2FdT2 + molality[i]*molality[j] * Phiprime[counterIJ];
4791  }
4792 #ifdef DEBUG_MODE
4793  if (m_debugCalc) {
4794  printf(" d2FdT2 = %10.6f \n", d2FdT2);
4795  }
4796 #endif
4797  }
4798  }
4799 #ifdef DEBUG_MODE
4800  if (m_debugCalc) {
4801  printf(" Step 8: \n");
4802  }
4803 #endif
4804 
4805  for (i = 1; i < m_kk; i++) {
4806 
4807  /*
4808  * -------- SUBSECTION FOR CALCULATING THE dACTCOEFFdT FOR CATIONS -----
4809  * --
4810  */
4811  if (charge[i] > 0) {
4812  // species i is the cation (positive) to calc the actcoeff
4813  zsqd2FdT2 = charge[i]*charge[i]*d2FdT2;
4814  sum1 = 0.0;
4815  sum2 = 0.0;
4816  sum3 = 0.0;
4817  sum4 = 0.0;
4818  sum5 = 0.0;
4819  for (j = 1; j < m_kk; j++) {
4820  /*
4821  * Find the counterIJ for the symmetric binary interaction
4822  */
4823  n = m_kk*i + j;
4824  counterIJ = m_CounterIJ[n];
4825 
4826  if (charge[j] < 0.0) {
4827  // sum over all anions
4828  sum1 = sum1 + molality[j]*
4829  (2.0*BMX_LL[counterIJ] + molarcharge*CMX_LL[counterIJ]);
4830  if (j < m_kk-1) {
4831  /*
4832  * This term is the ternary interaction involving the
4833  * non-duplicate sum over double anions, j, k, with
4834  * respect to the cation, i.
4835  */
4836  for (k = j+1; k < m_kk; k++) {
4837  // an inner sum over all anions
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];
4841  }
4842  }
4843  }
4844  }
4845 
4846 
4847  if (charge[j] > 0.0) {
4848  // sum over all cations
4849  if (j != i) {
4850  sum2 = sum2 + molality[j]*(2.0*Phi_LL[counterIJ]);
4851  }
4852  for (k = 1; k < m_kk; k++) {
4853  if (charge[k] < 0.0) {
4854  // two inner sums over anions
4855 
4856  n = k + j * m_kk + i * m_kk * m_kk;
4857  sum2 = sum2 + molality[j]*molality[k]*psi_ijk_LL[n];
4858  /*
4859  * Find the counterIJ for the j,k interaction
4860  */
4861  n = m_kk*j + k;
4862  counterIJ2 = m_CounterIJ[n];
4863  sum4 = sum4 + (fabs(charge[i])*
4864  molality[j]*molality[k]*CMX_LL[counterIJ2]);
4865  }
4866  }
4867  }
4868 
4869  /*
4870  * Handle neutral j species
4871  */
4872  if (charge[j] == 0) {
4873  sum5 = sum5 + molality[j]*2.0*m_Lambda_nj_LL(j,i);
4874  /*
4875  * Zeta interaction term
4876  */
4877  for (size_t k = 1; k < m_kk; k++) {
4878  if (charge[k] < 0.0) {
4879  size_t izeta = j;
4880  size_t jzeta = i;
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;
4885  }
4886  }
4887  }
4888  }
4889  }
4890  /*
4891  * Add all of the contributions up to yield the log of the
4892  * solute activity coefficients (molality scale)
4893  */
4894  m_d2lnActCoeffMolaldT2_Unscaled[i] =
4895  zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5;
4896 #ifdef DEBUG_MODE
4897  if (m_debugCalc) {
4898  sni = speciesName(i);
4899  printf(" %-16s d2lngammadT2[i]=%10.6f \n",
4900  sni.c_str(), m_d2lnActCoeffMolaldT2_Unscaled[i]);
4901  printf(" %12g %12g %12g %12g %12g %12g\n",
4902  zsqd2FdT2, sum1, sum2, sum3, sum4, sum5);
4903  }
4904 #endif
4905  }
4906 
4907 
4908  /*
4909  * ------ SUBSECTION FOR CALCULATING THE d2ACTCOEFFdT2 FOR ANIONS ------
4910  *
4911  */
4912  if (charge[i] < 0) {
4913  // species i is an anion (negative)
4914  zsqd2FdT2 = charge[i]*charge[i]*d2FdT2;
4915  sum1 = 0.0;
4916  sum2 = 0.0;
4917  sum3 = 0.0;
4918  sum4 = 0.0;
4919  sum5 = 0.0;
4920  for (j = 1; j < m_kk; j++) {
4921  /*
4922  * Find the counterIJ for the symmetric binary interaction
4923  */
4924  n = m_kk*i + j;
4925  counterIJ = m_CounterIJ[n];
4926 
4927  /*
4928  * For Anions, do the cation interactions.
4929  */
4930  if (charge[j] > 0) {
4931  sum1 = sum1 + molality[j]*
4932  (2.0*BMX_LL[counterIJ] + molarcharge*CMX_LL[counterIJ]);
4933  if (j < m_kk-1) {
4934  for (k = j+1; k < m_kk; k++) {
4935  // an inner sum over all cations
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];
4939  }
4940  }
4941  }
4942  }
4943 
4944  /*
4945  * For Anions, do the other anion interactions.
4946  */
4947  if (charge[j] < 0.0) {
4948  // sum over all anions
4949  if (j != i) {
4950  sum2 = sum2 + molality[j]*(2.0*Phi_LL[counterIJ]);
4951  }
4952  for (k = 1; k < m_kk; k++) {
4953  if (charge[k] > 0.0) {
4954  // two inner sums over cations
4955  n = k + j * m_kk + i * m_kk * m_kk;
4956  sum2 = sum2 + molality[j]*molality[k]*psi_ijk_LL[n];
4957  /*
4958  * Find the counterIJ for the symmetric binary interaction
4959  */
4960  n = m_kk*j + k;
4961  counterIJ2 = m_CounterIJ[n];
4962  sum4 = sum4 +
4963  (fabs(charge[i])*
4964  molality[j]*molality[k]*CMX_LL[counterIJ2]);
4965  }
4966  }
4967  }
4968 
4969  /*
4970  * for Anions, do the neutral species interaction
4971  */
4972  if (charge[j] == 0.0) {
4973  sum5 = sum5 + molality[j]*2.0*m_Lambda_nj_LL(j,i);
4974  /*
4975  * Zeta interaction term
4976  */
4977  for (size_t k = 1; k < m_kk; k++) {
4978  if (charge[k] > 0.0) {
4979  size_t izeta = j;
4980  size_t jzeta = k;
4981  size_t kzeta = i;
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;
4986  }
4987  }
4988  }
4989  }
4990  }
4991  m_d2lnActCoeffMolaldT2_Unscaled[i] =
4992  zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5;
4993 #ifdef DEBUG_MODE
4994  if (m_debugCalc) {
4995  sni = speciesName(i);
4996  printf(" %-16s d2lngammadT2[i]=%10.6f\n",
4997  sni.c_str(), m_d2lnActCoeffMolaldT2_Unscaled[i]);
4998  printf(" %12g %12g %12g %12g %12g %12g\n",
4999  zsqd2FdT2, sum1, sum2, sum3, sum4, sum5);
5000  }
5001 #endif
5002  }
5003  /*
5004  * ------ SUBSECTION FOR CALCULATING NEUTRAL SOLUTE ACT COEFF -------
5005  * ------ -> equations agree with my notes,
5006  * -> Equations agree with Pitzer,
5007  */
5008  if (charge[i] == 0.0) {
5009  sum1 = 0.0;
5010  sum3 = 0.0;
5011  for (j = 1; j < m_kk; j++) {
5012  sum1 = sum1 + molality[j]*2.0*m_Lambda_nj_LL(i,j);
5013  /*
5014  * Zeta term -> we piggyback on the psi term
5015  */
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];
5021  }
5022  }
5023  }
5024  }
5025  sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_LL[i];
5026  m_d2lnActCoeffMolaldT2_Unscaled[i] = sum1 + sum2 + sum3;
5027 #ifdef DEBUG_MODE
5028  if (m_debugCalc) {
5029  sni = speciesName(i);
5030  printf(" %-16s d2lngammadT2[i]=%10.6f \n",
5031  sni.c_str(), m_d2lnActCoeffMolaldT2_Unscaled[i]);
5032  }
5033 #endif
5034  }
5035 
5036  }
5037 #ifdef DEBUG_MODE
5038  if (m_debugCalc) {
5039  printf(" Step 9: \n");
5040  }
5041 #endif
5042 
5043  /*
5044  * ------ SUBSECTION FOR CALCULATING THE d2 OSMOTIC COEFF dT2 ---------
5045  *
5046  */
5047  sum1 = 0.0;
5048  sum2 = 0.0;
5049  sum3 = 0.0;
5050  sum4 = 0.0;
5051  sum5 = 0.0;
5052  double sum6 = 0.0;
5053  double sum7 = 0.0;
5054  /*
5055  * term1 is the temperature derivative of the
5056  * DH term in the osmotic coefficient expression
5057  * b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer
5058  * implementations.
5059  * Is = Ionic strength on the molality scale (units of (gmol/kg))
5060  * Aphi = A_Debye / 3 (units of sqrt(kg/gmol))
5061  */
5062  term1 = -d2AphidT2 * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
5063 
5064  for (j = 1; j < m_kk; j++) {
5065  /*
5066  * Loop Over Cations
5067  */
5068  if (charge[j] > 0.0) {
5069  for (k = 1; k < m_kk; k++) {
5070  if (charge[k] < 0.0) {
5071  /*
5072  * Find the counterIJ for the symmetric j,k binary interaction
5073  */
5074  n = m_kk*j + k;
5075  counterIJ = m_CounterIJ[n];
5076 
5077  sum1 = sum1 + molality[j]*molality[k]*
5078  (BphiMX_LL[counterIJ] + molarcharge*CMX_LL[counterIJ]);
5079  }
5080  }
5081 
5082  for (k = j+1; k < m_kk; k++) {
5083  if (j == (m_kk-1)) {
5084  // we should never reach this step
5085  printf("logic error 1 in Step 9 of hmw_act");
5086  exit(EXIT_FAILURE);
5087  }
5088  if (charge[k] > 0.0) {
5089  /*
5090  * Find the counterIJ for the symmetric j,k binary interaction
5091  * between 2 cations.
5092  */
5093  n = m_kk*j + k;
5094  counterIJ = m_CounterIJ[n];
5095  sum2 = sum2 + molality[j]*molality[k]*Phiphi_LL[counterIJ];
5096  for (m = 1; m < m_kk; m++) {
5097  if (charge[m] < 0.0) {
5098  // species m is an anion
5099  n = m + k * m_kk + j * m_kk * m_kk;
5100  sum2 = sum2 +
5101  molality[j]*molality[k]*molality[m]*psi_ijk_LL[n];
5102  }
5103  }
5104  }
5105  }
5106  }
5107 
5108  /*
5109  * Loop Over Anions
5110  */
5111  if (charge[j] < 0) {
5112  for (k = j+1; k < m_kk; k++) {
5113  if (j == m_kk-1) {
5114  // we should never reach this step
5115  printf("logic error 2 in Step 9 of hmw_act");
5116  exit(EXIT_FAILURE);
5117  }
5118  if (charge[k] < 0) {
5119  /*
5120  * Find the counterIJ for the symmetric j,k binary interaction
5121  * between two anions
5122  */
5123  n = m_kk*j + k;
5124  counterIJ = m_CounterIJ[n];
5125 
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;
5130  sum3 = sum3 +
5131  molality[j]*molality[k]*molality[m]*psi_ijk_LL[n];
5132  }
5133  }
5134  }
5135  }
5136  }
5137 
5138  /*
5139  * Loop Over Neutral Species
5140  */
5141  if (charge[j] == 0) {
5142  for (k = 1; k < m_kk; k++) {
5143  if (charge[k] < 0.0) {
5144  sum4 = sum4 + molality[j]*molality[k]*m_Lambda_nj_LL(j,k);
5145  }
5146  if (charge[k] > 0.0) {
5147  sum5 = sum5 + molality[j]*molality[k]*m_Lambda_nj_LL(j,k);
5148  }
5149  if (charge[k] == 0.0) {
5150  if (k > j) {
5151  sum6 = sum6 + molality[j]*molality[k]*m_Lambda_nj_LL(j,k);
5152  } else if (k == j) {
5153  sum6 = sum6 + 0.5 * molality[j]*molality[k]*m_Lambda_nj_LL(j,k);
5154  }
5155  }
5156  if (charge[k] < 0.0) {
5157  size_t izeta = j;
5158  for (m = 1; m < m_kk; m++) {
5159  if (charge[m] > 0.0) {
5160  size_t jzeta = m;
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;
5165  }
5166  }
5167  }
5168  }
5169  }
5170 
5171  sum7 += molality[j] * molality[j] * molality[j] * m_Mu_nnn_LL[j];
5172  }
5173  }
5174  sum_m_phi_minus_1 = 2.0 *
5175  (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
5176  /*
5177  * Calculate the osmotic coefficient from
5178  * osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i)
5179  */
5180  if (molalitysum > 1.0E-150) {
5181  d2_osmotic_coef_dT2 = 0.0 + (sum_m_phi_minus_1 / molalitysum);
5182  } else {
5183  d2_osmotic_coef_dT2 = 0.0;
5184  }
5185 #ifdef DEBUG_MODE
5186  if (m_debugCalc) {
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);
5192  }
5193 
5194  if (m_debugCalc) {
5195  printf(" Step 10: \n");
5196  }
5197 #endif
5198  d2_lnwateract_dT2 = -(m_weightSolvent/1000.0) * molalitysum * d2_osmotic_coef_dT2;
5199 
5200  /*
5201  * In Cantera, we define the activity coefficient of the solvent as
5202  *
5203  * act_0 = actcoeff_0 * Xmol_0
5204  *
5205  * We have just computed act_0. However, this routine returns
5206  * ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
5207  */
5208  m_d2lnActCoeffMolaldT2_Unscaled[0] = d2_lnwateract_dT2;
5209 
5210 #ifdef DEBUG_MODE
5211  if (m_debugCalc) {
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);
5215  }
5216 #endif
5217 }
5218 
5219 /********************************************************************************************/
5220 
5221 /*
5222  * s_update_dlnMolalityActCoeff_dP() (private, const )
5223  *
5224  * Using internally stored values, this function calculates
5225  * the pressure derivative of the logarithm of the
5226  * activity coefficient for all species in the mechanism.
5227  *
5228  * We assume that the activity coefficients are current.
5229  *
5230  * solvent activity coefficient is on the molality
5231  * scale. Its derivative is too.
5232  */
5233 void HMWSoln::s_update_dlnMolalityActCoeff_dP() const
5234 {
5235  m_dlnActCoeffMolaldP_Unscaled.assign(m_kk, 0.0);
5236  s_updatePitzer_dlnMolalityActCoeff_dP();
5237 
5238 
5239  for (size_t k = 1; k < m_kk; k++) {
5240  if (CROP_speciesCropped_[k] == 2) {
5241  m_dlnActCoeffMolaldP_Unscaled[k] = 0.0;
5242  }
5243  }
5244 
5245  if (CROP_speciesCropped_[0]) {
5246  m_dlnActCoeffMolaldP_Unscaled[0] = 0.0;
5247  }
5248 
5249 
5250  s_updateScaling_pHScaling_dP();
5251 }
5252 
5253 /*
5254  * s_updatePitzer_dlnMolalityActCoeff_dP() (private, const )
5255  *
5256  * Using internally stored values, this function calculates
5257  * the pressure derivative of the logarithm of the
5258  * activity coefficient for all species in the mechanism.
5259  * This is an internal routine
5260  *
5261  * We assume that the activity coefficients are current.
5262  *
5263  * It may be assumed that the
5264  * Pitzer activity coefficient and first deriv routine are called immediately
5265  * preceding the calling of this routine. Therefore, some
5266  * quantities do not need to be recalculated in this routine.
5267  *
5268  * solvent activity coefficient is on the molality
5269  * scale. Its derivatives are too.
5270  */
5271 void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
5272 {
5273 
5274  /*
5275  * HKM -> Assumption is made that the solvent is
5276  * species 0.
5277  */
5278 #ifdef DEBUG_MODE
5279  m_debugCalc = 0;
5280 #endif
5281  if (m_indexSolvent != 0) {
5282  printf("Wrong index solvent value!\n");
5283  exit(EXIT_FAILURE);
5284  }
5285 
5286  std::string sni, snj, snk;
5287 
5288  const double* molality = DATA_PTR(m_molalitiesCropped);
5289  const double* charge = DATA_PTR(m_speciesCharge);
5290  const double* beta0MX_P = DATA_PTR(m_Beta0MX_ij_P);
5291  const double* beta1MX_P = DATA_PTR(m_Beta1MX_ij_P);
5292  const double* beta2MX_P = DATA_PTR(m_Beta2MX_ij_P);
5293  const double* CphiMX_P = DATA_PTR(m_CphiMX_ij_P);
5294  const double* thetaij_P = DATA_PTR(m_Theta_ij_P);
5295  const double* alpha1MX = DATA_PTR(m_Alpha1MX_ij);
5296  const double* alpha2MX = DATA_PTR(m_Alpha2MX_ij);
5297  const double* psi_ijk_P = DATA_PTR(m_Psi_ijk_P);
5298 
5299  /*
5300  * Local variables defined by Coltrin
5301  */
5302  double etheta[5][5], etheta_prime[5][5], sqrtIs;
5303  /*
5304  * Molality based ionic strength of the solution
5305  */
5306  double Is = 0.0;
5307  /*
5308  * Molarcharge of the solution: In Pitzer's notation,
5309  * this is his variable called "Z".
5310  */
5311  double molarcharge = 0.0;
5312  /*
5313  * molalitysum is the sum of the molalities over all solutes,
5314  * even those with zero charge.
5315  */
5316  double molalitysum = 0.0;
5317 
5318  double* gfunc = DATA_PTR(m_gfunc_IJ);
5319  double* g2func = DATA_PTR(m_g2func_IJ);
5320  double* hfunc = DATA_PTR(m_hfunc_IJ);
5321  double* h2func = DATA_PTR(m_h2func_IJ);
5322  double* BMX_P = DATA_PTR(m_BMX_IJ_P);
5323  double* BprimeMX_P= DATA_PTR(m_BprimeMX_IJ_P);
5324  double* BphiMX_P = DATA_PTR(m_BphiMX_IJ_P);
5325  double* Phi_P = DATA_PTR(m_Phi_IJ_P);
5326  double* Phiprime = DATA_PTR(m_Phiprime_IJ);
5327  double* Phiphi_P = DATA_PTR(m_PhiPhi_IJ_P);
5328  double* CMX_P = DATA_PTR(m_CMX_IJ_P);
5329 
5330  double x1, x2;
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;
5334 
5335  int z1, z2;
5336  size_t n, i, j, k, m, counterIJ, counterIJ2;
5337 
5338  double currTemp = temperature();
5339  double currPres = pressure();
5340 
5341 #ifdef DEBUG_MODE
5342  if (m_debugCalc) {
5343  printf("\n Debugging information from "
5344  "s_Pitzer_dlnMolalityActCoeff_dP()\n");
5345  }
5346 #endif
5347  /*
5348  * Make sure the counter variables are setup
5349  */
5350  counterIJ_setup();
5351 
5352  /*
5353  * ---------- Calculate common sums over solutes ---------------------
5354  */
5355  for (n = 1; n < m_kk; n++) {
5356  // ionic strength
5357  Is += charge[n] * charge[n] * molality[n];
5358  // total molar charge
5359  molarcharge += fabs(charge[n]) * molality[n];
5360  molalitysum += molality[n];
5361  }
5362  Is *= 0.5;
5363 
5364  /*
5365  * Store the ionic molality in the object for reference.
5366  */
5367  m_IionicMolality = Is;
5368  sqrtIs = sqrt(Is);
5369 #ifdef DEBUG_MODE
5370  if (m_debugCalc) {
5371  printf(" Step 1: \n");
5372  printf(" ionic strenth = %14.7le \n total molar "
5373  "charge = %14.7le \n", Is, molarcharge);
5374  }
5375 #endif
5376 
5377  /*
5378  * The following call to calc_lambdas() calculates all 16 elements
5379  * of the elambda and elambda1 arrays, given the value of the
5380  * ionic strength (Is)
5381  */
5382  calc_lambdas(Is);
5383 
5384 
5385  /*
5386  * ----- Step 2: Find the coefficients E-theta and -------------------
5387  * E-thetaprime for all combinations of positive
5388  * unlike charges up to 4
5389  */
5390 #ifdef DEBUG_MODE
5391  if (m_debugCalc) {
5392  printf(" Step 2: \n");
5393  }
5394 #endif
5395  for (z1 = 1; z1 <=4; z1++) {
5396  for (z2 =1; z2 <=4; z2++) {
5397  calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
5398 #ifdef DEBUG_MODE
5399  if (m_debugCalc) {
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]);
5402  }
5403 #endif
5404  }
5405  }
5406 
5407 #ifdef DEBUG_MODE
5408  if (m_debugCalc) {
5409  printf(" Step 3: \n");
5410  printf(" Species Species g(x) "
5411  " hfunc(x) \n");
5412  }
5413 #endif
5414 
5415  /*
5416  *
5417  * calculate g(x) and hfunc(x) for each cation-anion pair MX
5418  * In the original literature, hfunc, was called gprime. However,
5419  * it's not the derivative of g(x), so I renamed it.
5420  */
5421  for (i = 1; i < (m_kk - 1); i++) {
5422  for (j = (i+1); j < m_kk; j++) {
5423  /*
5424  * Find the counterIJ for the symmetric binary interaction
5425  */
5426  n = m_kk*i + j;
5427  counterIJ = m_CounterIJ[n];
5428  /*
5429  * Only loop over oppositely charge species
5430  */
5431  if (charge[i]*charge[j] < 0) {
5432  /*
5433  * x is a reduced function variable
5434  */
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);
5440  } else {
5441  gfunc[counterIJ] = 0.0;
5442  hfunc[counterIJ] = 0.0;
5443  }
5444 
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);
5451  } else {
5452  g2func[counterIJ] = 0.0;
5453  h2func[counterIJ] = 0.0;
5454  }
5455  }
5456  } else {
5457  gfunc[counterIJ] = 0.0;
5458  hfunc[counterIJ] = 0.0;
5459  }
5460 #ifdef DEBUG_MODE
5461  if (m_debugCalc) {
5462  sni = speciesName(i);
5463  snj = speciesName(j);
5464  printf(" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(),
5465  gfunc[counterIJ], hfunc[counterIJ]);
5466  }
5467 #endif
5468  }
5469  }
5470 
5471  /*
5472  * ------- SUBSECTION TO CALCULATE BMX_P, BprimeMX_P, BphiMX_P ----------
5473  * ------- These are now temperature derivatives of the
5474  * previously calculated quantities.
5475  */
5476 #ifdef DEBUG_MODE
5477  if (m_debugCalc) {
5478  printf(" Step 4: \n");
5479  printf(" Species Species BMX "
5480  "BprimeMX BphiMX \n");
5481  }
5482 #endif
5483 
5484  for (i = 1; i < m_kk - 1; i++) {
5485  for (j = i+1; j < m_kk; j++) {
5486  /*
5487  * Find the counterIJ for the symmetric binary interaction
5488  */
5489  n = m_kk*i + j;
5490  counterIJ = m_CounterIJ[n];
5491  /*
5492  * both species have a non-zero charge, and one is positive
5493  * and the other is negative
5494  */
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];
5499 #ifdef DEBUG_MODE
5500  if (m_debugCalc) {
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]);
5504  }
5505 #endif
5506  if (Is > 1.0E-150) {
5507  BprimeMX_P[counterIJ] = (beta1MX_P[counterIJ] * hfunc[counterIJ]/Is +
5508  beta2MX_P[counterIJ] * h2func[counterIJ]/Is);
5509  } else {
5510  BprimeMX_P[counterIJ] = 0.0;
5511  }
5512  BphiMX_P[counterIJ] = BMX_P[counterIJ] + Is*BprimeMX_P[counterIJ];
5513  } else {
5514  BMX_P[counterIJ] = 0.0;
5515  BprimeMX_P[counterIJ] = 0.0;
5516  BphiMX_P[counterIJ] = 0.0;
5517  }
5518 #ifdef DEBUG_MODE
5519  if (m_debugCalc) {
5520  sni = speciesName(i);
5521  snj = speciesName(j);
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]);
5525  }
5526 #endif
5527  }
5528  }
5529 
5530  /*
5531  * --------- SUBSECTION TO CALCULATE CMX_P ----------
5532  * ---------
5533  */
5534 #ifdef DEBUG_MODE
5535  if (m_debugCalc) {
5536  printf(" Step 5: \n");
5537  printf(" Species Species CMX \n");
5538  }
5539 #endif
5540  for (i = 1; i < m_kk-1; i++) {
5541  for (j = i+1; j < m_kk; j++) {
5542  /*
5543  * Find the counterIJ for the symmetric binary interaction
5544  */
5545  n = m_kk*i + j;
5546  counterIJ = m_CounterIJ[n];
5547  /*
5548  * both species have a non-zero charge, and one is positive
5549  * and the other is negative
5550  */
5551  if (charge[i]*charge[j] < 0.0) {
5552  CMX_P[counterIJ] = CphiMX_P[counterIJ]/
5553  (2.0* sqrt(fabs(charge[i]*charge[j])));
5554  } else {
5555  CMX_P[counterIJ] = 0.0;
5556  }
5557 #ifdef DEBUG_MODE
5558  if (m_debugCalc) {
5559  sni = speciesName(i);
5560  snj = speciesName(j);
5561  printf(" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(),
5562  CMX_P[counterIJ]);
5563  }
5564 #endif
5565  }
5566  }
5567 
5568  /*
5569  * ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ----------
5570  * --------
5571  */
5572 #ifdef DEBUG_MODE
5573  if (m_debugCalc) {
5574  printf(" Step 6: \n");
5575  printf(" Species Species Phi_ij "
5576  " Phiprime_ij Phi^phi_ij \n");
5577  }
5578 #endif
5579  for (i = 1; i < m_kk-1; i++) {
5580  for (j = i+1; j < m_kk; j++) {
5581  /*
5582  * Find the counterIJ for the symmetric binary interaction
5583  */
5584  n = m_kk*i + j;
5585  counterIJ = m_CounterIJ[n];
5586  /*
5587  * both species have a non-zero charge, and one is positive
5588  * and the other is negative
5589  */
5590  if (charge[i]*charge[j] > 0) {
5591  z1 = (int) fabs(charge[i]);
5592  z2 = (int) fabs(charge[j]);
5593  //Phi[counterIJ] = thetaij_L[counterIJ] + etheta[z1][z2];
5594  Phi_P[counterIJ] = thetaij_P[counterIJ];
5595  //Phiprime[counterIJ] = etheta_prime[z1][z2];
5596  Phiprime[counterIJ] = 0.0;
5597  Phiphi_P[counterIJ] = Phi_P[counterIJ] + Is * Phiprime[counterIJ];
5598  } else {
5599  Phi_P[counterIJ] = 0.0;
5600  Phiprime[counterIJ] = 0.0;
5601  Phiphi_P[counterIJ] = 0.0;
5602  }
5603 #ifdef DEBUG_MODE
5604  if (m_debugCalc) {
5605  sni = speciesName(i);
5606  snj = speciesName(j);
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]);
5610  }
5611 #endif
5612  }
5613  }
5614 
5615  /*
5616  * ----------- SUBSECTION FOR CALCULATION OF dFdT ---------------------
5617  */
5618 #ifdef DEBUG_MODE
5619  if (m_debugCalc) {
5620  printf(" Step 7: \n");
5621  }
5622 #endif
5623  // A_Debye_Huckel = 0.5092; (units = sqrt(kg/gmol))
5624  // A_Debye_Huckel = 0.5107; <- This value is used to match GWB data
5625  // ( A * ln(10) = 1.17593)
5626  // Aphi = A_Debye_Huckel * 2.30258509 / 3.0;
5627 
5628  double dA_DebyedP = dA_DebyedP_TP(currTemp, currPres);
5629  double dAphidP = dA_DebyedP /3.0;
5630 #ifdef DEBUG_MODE
5631  //dAphidT = 0.0;
5632 #endif
5633  //F = -Aphi * ( sqrt(Is) / (1.0 + 1.2*sqrt(Is))
5634  // + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
5635  //dAphidT = Al / (4.0 * GasConstant * T * T);
5636  dFdP = -dAphidP * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
5637  + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
5638 #ifdef DEBUG_MODE
5639  if (m_debugCalc) {
5640  printf(" initial value of dFdP = %10.6f \n", dFdP);
5641  }
5642 #endif
5643  for (i = 1; i < m_kk-1; i++) {
5644  for (j = i+1; j < m_kk; j++) {
5645  /*
5646  * Find the counterIJ for the symmetric binary interaction
5647  */
5648  n = m_kk*i + j;
5649  counterIJ = m_CounterIJ[n];
5650  /*
5651  * both species have a non-zero charge, and one is positive
5652  * and the other is negative
5653  */
5654  if (charge[i]*charge[j] < 0) {
5655  dFdP = dFdP + molality[i]*molality[j] * BprimeMX_P[counterIJ];
5656  }
5657  /*
5658  * Both species have a non-zero charge, and they
5659  * have the same sign, e.g., both positive or both negative.
5660  */
5661  if (charge[i]*charge[j] > 0) {
5662  dFdP = dFdP + molality[i]*molality[j] * Phiprime[counterIJ];
5663  }
5664 #ifdef DEBUG_MODE
5665  if (m_debugCalc) {
5666  printf(" dFdP = %10.6f \n", dFdP);
5667  }
5668 #endif
5669  }
5670  }
5671 #ifdef DEBUG_MODE
5672  if (m_debugCalc) {
5673  printf(" Step 8: \n");
5674  }
5675 #endif
5676 
5677 
5678  for (i = 1; i < m_kk; i++) {
5679 
5680  /*
5681  * -------- SUBSECTION FOR CALCULATING THE dACTCOEFFdP FOR CATIONS -----
5682  * --
5683  */
5684  if (charge[i] > 0) {
5685  // species i is the cation (positive) to calc the actcoeff
5686  zsqdFdP = charge[i]*charge[i]*dFdP;
5687  sum1 = 0.0;
5688  sum2 = 0.0;
5689  sum3 = 0.0;
5690  sum4 = 0.0;
5691  sum5 = 0.0;
5692  for (j = 1; j < m_kk; j++) {
5693  /*
5694  * Find the counterIJ for the symmetric binary interaction
5695  */
5696  n = m_kk*i + j;
5697  counterIJ = m_CounterIJ[n];
5698 
5699  if (charge[j] < 0.0) {
5700  // sum over all anions
5701  sum1 = sum1 + molality[j]*
5702  (2.0*BMX_P[counterIJ] + molarcharge*CMX_P[counterIJ]);
5703  if (j < m_kk-1) {
5704  /*
5705  * This term is the ternary interaction involving the
5706  * non-duplicate sum over double anions, j, k, with
5707  * respect to the cation, i.
5708  */
5709  for (k = j+1; k < m_kk; k++) {
5710  // an inner sum over all anions
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];
5714  }
5715  }
5716  }
5717  }
5718 
5719 
5720  if (charge[j] > 0.0) {
5721  // sum over all cations
5722  if (j != i) {
5723  sum2 = sum2 + molality[j]*(2.0*Phi_P[counterIJ]);
5724  }
5725  for (k = 1; k < m_kk; k++) {
5726  if (charge[k] < 0.0) {
5727  // two inner sums over anions
5728 
5729  n = k + j * m_kk + i * m_kk * m_kk;
5730  sum2 = sum2 + molality[j]*molality[k]*psi_ijk_P[n];
5731  /*
5732  * Find the counterIJ for the j,k interaction
5733  */
5734  n = m_kk*j + k;
5735  counterIJ2 = m_CounterIJ[n];
5736  sum4 = sum4 + (fabs(charge[i])*
5737  molality[j]*molality[k]*CMX_P[counterIJ2]);
5738  }
5739  }
5740  }
5741 
5742  /*
5743  * for Anions, do the neutral species interaction
5744  */
5745  if (charge[j] == 0) {
5746  sum5 = sum5 + molality[j]*2.0*m_Lambda_nj_P(j,i);
5747  /*
5748  * Zeta interaction term
5749  */
5750  for (size_t k = 1; k < m_kk; k++) {
5751  if (charge[k] < 0.0) {
5752  size_t izeta = j;
5753  size_t jzeta = i;
5754  n = izeta * m_kk * m_kk + jzeta * m_kk + k;
5755  double zeta_P = psi_ijk_P[n];
5756  if (zeta_P != 0.0) {
5757  sum5 = sum5 + molality[j]*molality[k]*zeta_P;
5758  }
5759  }
5760  }
5761  }
5762  }
5763 
5764  /*
5765  * Add all of the contributions up to yield the log of the
5766  * solute activity coefficients (molality scale)
5767  */
5768  m_dlnActCoeffMolaldP_Unscaled[i] =
5769  zsqdFdP + sum1 + sum2 + sum3 + sum4 + sum5;
5770 
5771 #ifdef DEBUG_MODE
5772  if (m_debugCalc) {
5773  sni = speciesName(i);
5774  printf(" %-16s lngamma[i]=%10.6f \n",
5775  sni.c_str(), m_dlnActCoeffMolaldP_Unscaled[i]);
5776  printf(" %12g %12g %12g %12g %12g %12g\n",
5777  zsqdFdP, sum1, sum2, sum3, sum4, sum5);
5778  }
5779 #endif
5780  }
5781 
5782 
5783  /*
5784  * ------ SUBSECTION FOR CALCULATING THE dACTCOEFFdP FOR ANIONS ------
5785  *
5786  */
5787  if (charge[i] < 0) {
5788  // species i is an anion (negative)
5789  zsqdFdP = charge[i]*charge[i]*dFdP;
5790  sum1 = 0.0;
5791  sum2 = 0.0;
5792  sum3 = 0.0;
5793  sum4 = 0.0;
5794  sum5 = 0.0;
5795  for (j = 1; j < m_kk; j++) {
5796  /*
5797  * Find the counterIJ for the symmetric binary interaction
5798  */
5799  n = m_kk*i + j;
5800  counterIJ = m_CounterIJ[n];
5801 
5802  /*
5803  * For Anions, do the cation interactions.
5804  */
5805  if (charge[j] > 0) {
5806  sum1 = sum1 + molality[j]*
5807  (2.0*BMX_P[counterIJ] + molarcharge*CMX_P[counterIJ]);
5808  if (j < m_kk-1) {
5809  for (k = j+1; k < m_kk; k++) {
5810  // an inner sum over all cations
5811  if (charge[k] > 0) {
5812  n = k + j * m_kk + i * m_kk * m_kk;
5813  sum3 = sum3 + molality[j]*molality[k]*psi_ijk_P[n];
5814  }
5815  }
5816  }
5817  }
5818 
5819  /*
5820  * For Anions, do the other anion interactions.
5821  */
5822  if (charge[j] < 0.0) {
5823  // sum over all anions
5824  if (j != i) {
5825  sum2 = sum2 + molality[j]*(2.0*Phi_P[counterIJ]);
5826  }
5827  for (k = 1; k < m_kk; k++) {
5828  if (charge[k] > 0.0) {
5829  // two inner sums over cations
5830  n = k + j * m_kk + i * m_kk * m_kk;
5831  sum2 = sum2 + molality[j]*molality[k]*psi_ijk_P[n];
5832  /*
5833  * Find the counterIJ for the symmetric binary interaction
5834  */
5835  n = m_kk*j + k;
5836  counterIJ2 = m_CounterIJ[n];
5837  sum4 = sum4 +
5838  (fabs(charge[i])*
5839  molality[j]*molality[k]*CMX_P[counterIJ2]);
5840  }
5841  }
5842  }
5843 
5844  /*
5845  * for Anions, do the neutral species interaction
5846  */
5847  if (charge[j] == 0.0) {
5848  sum5 = sum5 + molality[j]*2.0*m_Lambda_nj_P(j,i);
5849  /*
5850  * Zeta interaction term
5851  */
5852  for (size_t k = 1; k < m_kk; k++) {
5853  if (charge[k] > 0.0) {
5854  size_t izeta = j;
5855  size_t jzeta = k;
5856  size_t kzeta = i;
5857  n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
5858  double zeta_P = psi_ijk_P[n];
5859  if (zeta_P != 0.0) {
5860  sum5 = sum5 + molality[j]*molality[k]*zeta_P;
5861  }
5862  }
5863  }
5864  }
5865  }
5866  m_dlnActCoeffMolaldP_Unscaled[i] =
5867  zsqdFdP + sum1 + sum2 + sum3 + sum4 + sum5;
5868 #ifdef DEBUG_MODE
5869  if (m_debugCalc) {
5870  sni = speciesName(i);
5871  printf(" %-16s lndactcoeffmolaldP[i]=%10.6f \n",
5872  sni.c_str(), m_dlnActCoeffMolaldP_Unscaled[i]);
5873  printf(" %12g %12g %12g %12g %12g %12g\n",
5874  zsqdFdP, sum1, sum2, sum3, sum4, sum5);
5875  }
5876 #endif
5877  }
5878 
5879  /*
5880  * ------ SUBSECTION FOR CALCULATING d NEUTRAL SOLUTE ACT COEFF dP -------
5881  */
5882  if (charge[i] == 0.0) {
5883  sum1 = 0.0;
5884  sum3 = 0.0;
5885  for (j = 1; j < m_kk; j++) {
5886  sum1 += molality[j]*2.0*m_Lambda_nj_P(i,j);
5887  /*
5888  * Zeta term -> we piggyback on the psi term
5889  */
5890  if (charge[j] > 0.0) {
5891  for (k = 1; k < m_kk; k++) {
5892  if (charge[k] < 0.0) {
5893  n = k + j * m_kk + i * m_kk * m_kk;
5894  sum3 = sum3 + molality[j]*molality[k]*psi_ijk_P[n];
5895  }
5896  }
5897  }
5898  }
5899  sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_P[i];
5900  m_dlnActCoeffMolaldP_Unscaled[i] = sum1 + sum2 + sum3;
5901 #ifdef DEBUG_MODE
5902  if (m_debugCalc) {
5903  sni = speciesName(i);
5904  printf(" %-16s dlnActCoeffMolaldP[i]=%10.6f \n",
5905  sni.c_str(), m_dlnActCoeffMolaldP_Unscaled[i]);
5906  }
5907 #endif
5908  }
5909 
5910  }
5911 #ifdef DEBUG_MODE
5912  if (m_debugCalc) {
5913  printf(" Step 9: \n");
5914  }
5915 #endif
5916 
5917  /*
5918  * ------ SUBSECTION FOR CALCULATING THE d OSMOTIC COEFF dP ---------
5919  *
5920  */
5921  sum1 = 0.0;
5922  sum2 = 0.0;
5923  sum3 = 0.0;
5924  sum4 = 0.0;
5925  sum5 = 0.0;
5926  double sum6 = 0.0;
5927  double sum7 = 0.0;
5928  /*
5929  * term1 is the temperature derivative of the
5930  * DH term in the osmotic coefficient expression
5931  * b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer
5932  * implementations.
5933  * Is = Ionic strength on the molality scale (units of (gmol/kg))
5934  * Aphi = A_Debye / 3 (units of sqrt(kg/gmol))
5935  */
5936  term1 = -dAphidP * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
5937 
5938  for (j = 1; j < m_kk; j++) {
5939  /*
5940  * Loop Over Cations
5941  */
5942  if (charge[j] > 0.0) {
5943  for (k = 1; k < m_kk; k++) {
5944  if (charge[k] < 0.0) {
5945  /*
5946  * Find the counterIJ for the symmetric j,k binary interaction
5947  */
5948  n = m_kk*j + k;
5949  counterIJ = m_CounterIJ[n];
5950 
5951  sum1 = sum1 + molality[j]*molality[k]*
5952  (BphiMX_P[counterIJ] + molarcharge*CMX_P[counterIJ]);
5953  }
5954  }
5955 
5956  for (k = j+1; k < m_kk; k++) {
5957  if (j == (m_kk-1)) {
5958  // we should never reach this step
5959  printf("logic error 1 in Step 9 of hmw_act");
5960  exit(EXIT_FAILURE);
5961  }
5962  if (charge[k] > 0.0) {
5963  /*
5964  * Find the counterIJ for the symmetric j,k binary interaction
5965  * between 2 cations.
5966  */
5967  n = m_kk*j + k;
5968  counterIJ = m_CounterIJ[n];
5969  sum2 = sum2 + molality[j]*molality[k]*Phiphi_P[counterIJ];
5970  for (m = 1; m < m_kk; m++) {
5971  if (charge[m] < 0.0) {
5972  // species m is an anion
5973  n = m + k * m_kk + j * m_kk * m_kk;
5974  sum2 = sum2 +
5975  molality[j]*molality[k]*molality[m]*psi_ijk_P[n];
5976  }
5977  }
5978  }
5979  }
5980  }
5981 
5982 
5983  /*
5984  * Loop Over Anions
5985  */
5986  if (charge[j] < 0) {
5987  for (k = j+1; k < m_kk; k++) {
5988  if (j == m_kk-1) {
5989  // we should never reach this step
5990  printf("logic error 2 in Step 9 of hmw_act");
5991  exit(EXIT_FAILURE);
5992  }
5993  if (charge[k] < 0) {
5994  /*
5995  * Find the counterIJ for the symmetric j,k binary interaction
5996  * between two anions
5997  */
5998  n = m_kk*j + k;
5999  counterIJ = m_CounterIJ[n];
6000 
6001  sum3 = sum3 + molality[j]*molality[k]*Phiphi_P[counterIJ];
6002  for (m = 1; m < m_kk; m++) {
6003  if (charge[m] > 0.0) {
6004  n = m + k * m_kk + j * m_kk * m_kk;
6005  sum3 = sum3 +
6006  molality[j]*molality[k]*molality[m]*psi_ijk_P[n];
6007  }
6008  }
6009  }
6010  }
6011  }
6012 
6013  /*
6014  * Loop Over Neutral Species
6015  */
6016  if (charge[j] == 0) {
6017  for (k = 1; k < m_kk; k++) {
6018  if (charge[k] < 0.0) {
6019  sum4 = sum4 + molality[j]*molality[k]*m_Lambda_nj_P(j,k);
6020  }
6021  if (charge[k] > 0.0) {
6022  sum5 = sum5 + molality[j]*molality[k]*m_Lambda_nj_P(j,k);
6023  }
6024  if (charge[k] == 0.0) {
6025  if (k > j) {
6026  sum6 = sum6 + molality[j]*molality[k]*m_Lambda_nj_P(j,k);
6027  } else if (k == j) {
6028  sum6 = sum6 + 0.5 * molality[j]*molality[k]*m_Lambda_nj_P(j,k);
6029  }
6030  }
6031  if (charge[k] < 0.0) {
6032  size_t izeta = j;
6033  for (m = 1; m < m_kk; m++) {
6034  if (charge[m] > 0.0) {
6035  size_t jzeta = m;
6036  n = k + jzeta * m_kk + izeta * m_kk * m_kk;
6037  double zeta_P = psi_ijk_P[n];
6038  if (zeta_P != 0.0) {
6039  sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_P;
6040  }
6041  }
6042  }
6043  }
6044  }
6045 
6046  sum7 += molality[j] * molality[j] * molality[j] * m_Mu_nnn_P[j];
6047  }
6048  }
6049  sum_m_phi_minus_1 = 2.0 *
6050  (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
6051 
6052  /*
6053  * Calculate the osmotic coefficient from
6054  * osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i)
6055  */
6056  if (molalitysum > 1.0E-150) {
6057  d_osmotic_coef_dP = 0.0 + (sum_m_phi_minus_1 / molalitysum);
6058  } else {
6059  d_osmotic_coef_dP = 0.0;
6060  }
6061 #ifdef DEBUG_MODE
6062  if (m_debugCalc) {
6063  printf(" term1=%10.6f sum1=%10.6f sum2=%10.6f "
6064  "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n",
6065  term1, sum1, sum2, sum3, sum4, sum5);
6066  printf(" sum_m_phi_minus_1=%10.6f d_osmotic_coef_dP =%10.6f\n",
6067  sum_m_phi_minus_1, d_osmotic_coef_dP);
6068  }
6069 
6070  if (m_debugCalc) {
6071  printf(" Step 10: \n");
6072  }
6073 #endif
6074  d_lnwateract_dP = -(m_weightSolvent/1000.0) * molalitysum * d_osmotic_coef_dP;
6075 
6076 
6077  /*
6078  * In Cantera, we define the activity coefficient of the solvent as
6079  *
6080  * act_0 = actcoeff_0 * Xmol_0
6081  *
6082  * We have just computed act_0. However, this routine returns
6083  * ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
6084  */
6085  //double xmolSolvent = moleFraction(m_indexSolvent);
6086  m_dlnActCoeffMolaldP_Unscaled[0] = d_lnwateract_dP;
6087 #ifdef DEBUG_MODE
6088  if (m_debugCalc) {
6089  double d_wateract_dP = exp(d_lnwateract_dP);
6090  printf(" d_ln_a_water_dP = %10.6f d_a_water_dP=%10.6f\n\n",
6091  d_lnwateract_dP, d_wateract_dP);
6092  }
6093 #endif
6094 
6095 }
6096 
6097 /**********************************************************************************************/
6098 
6099 /*
6100  * Calculate the lambda interactions.
6101  *
6102  * Calculate E-lambda terms for charge combinations of like sign,
6103  * using method of Pitzer (1975).
6104  *
6105  * This code snippet is included from Bethke, Appendix 2.
6106  */
6107 void HMWSoln::calc_lambdas(double is) const
6108 {
6109  double aphi, dj, jfunc, jprime, t, x, zprod;
6110  int i, ij, j;
6111  /*
6112  * Coefficients c1-c4 are used to approximate
6113  * the integral function "J";
6114  * aphi is the Debye-Huckel constant at 25 C
6115  */
6116 
6117  double c1 = 4.581, c2 = 0.7237, c3 = 0.0120, c4 = 0.528;
6118 
6119  aphi = 0.392; /* Value at 25 C */
6120 #ifdef DEBUG_MODE
6121  if (m_debugCalc) {
6122  printf(" Is = %g\n", is);
6123  }
6124 #endif
6125  if (is < 1.0E-150) {
6126  for (i = 0; i < 17; i++) {
6127  elambda[i] = 0.0;
6128  elambda1[i] = 0.0;
6129  }
6130  return;
6131  }
6132  /*
6133  * Calculate E-lambda terms for charge combinations of like sign,
6134  * using method of Pitzer (1975). Charges up to 4 are calculated.
6135  */
6136 
6137  for (i=1; i<=4; i++) {
6138  for (j=i; j<=4; j++) {
6139  ij = i*j;
6140  /*
6141  * calculate the product of the charges
6142  */
6143  zprod = (double)ij;
6144  /*
6145  * calculate Xmn (A1) from Harvie, Weare (1980).
6146  */
6147  x = 6.0* zprod * aphi * sqrt(is); /* eqn 23 */
6148 
6149  jfunc = x / (4.0 + c1*pow(x,-c2)*exp(-c3*pow(x,c4))); /* eqn 47 */
6150 
6151  t = c3 * c4 * pow(x,c4);
6152  dj = c1* pow(x,(-c2-1.0)) * (c2+t) * exp(-c3*pow(x,c4));
6153  jprime = (jfunc/x)*(1.0 + jfunc*dj);
6154 
6155  elambda[ij] = zprod*jfunc / (4.0*is); /* eqn 14 */
6156  elambda1[ij] = (3.0*zprod*zprod*aphi*jprime/(4.0*sqrt(is))
6157  - elambda[ij])/is;
6158 #ifdef DEBUG_MODE
6159  if (m_debugCalc) {
6160  printf(" ij = %d, elambda = %g, elambda1 = %g\n",
6161  ij, elambda[ij], elambda1[ij]);
6162  }
6163 #endif
6164  }
6165  }
6166 }
6167 
6168 /*
6169  * Calculate the etheta interaction.
6170  * This interaction accounts for the mixing effects of like-signed
6171  * ions with different charges. There is fairly extensive literature
6172  * on this effect. See the notes.
6173  * This interaction will be nonzero for species with the same charge.
6174  *
6175  * This code snippet is included from Bethke, Appendix 2.
6176  */
6177 void HMWSoln::calc_thetas(int z1, int z2,
6178  double* etheta, double* etheta_prime) const
6179 {
6180  int i, j;
6181  double f1, f2;
6182 
6183  /*
6184  * Calculate E-theta(i) and E-theta'(I) using method of
6185  * Pitzer (1987)
6186  */
6187  i = abs(z1);
6188  j = abs(z2);
6189 
6190 #ifdef DEBUG_MODE
6191  if (i > 4 || j > 4) {
6192  printf("we shouldn't be here\n");
6193  exit(EXIT_FAILURE);
6194  }
6195 #endif
6196 
6197  if ((i == 0) || (j == 0)) {
6198  printf("ERROR calc_thetas called with one species being neutral\n");
6199  exit(EXIT_FAILURE);
6200  }
6201 
6202  /*
6203  * Check to see if the charges are of opposite sign. If they are of
6204  * opposite sign then their etheta interaction is zero.
6205  */
6206  if (z1*z2 < 0) {
6207  *etheta = 0.0;
6208  *etheta_prime = 0.0;
6209  }
6210  /*
6211  * Actually calculate the interaction.
6212  */
6213  else {
6214  f1 = (double)i / (2.0 * j);
6215  f2 = (double)j / (2.0 * i);
6216  *etheta = elambda[i*j] - f1*elambda[j*j] - f2*elambda[i*i];
6217  *etheta_prime = elambda1[i*j] - f1*elambda1[j*j] - f2*elambda1[i*i];
6218  }
6219 }
6220 
6221 // This function will be called to update the internally stored
6222 // natural logarithm of the molality activity coefficients
6223 /*
6224  * Normally they are all one. However, sometimes they are not,
6225  * due to stability schemes
6226  *
6227  * gamma_k_molar = gamma_k_molal / Xmol_solvent
6228  *
6229  * gamma_o_molar = gamma_o_molal
6230  */
6231 void HMWSoln::s_updateIMS_lnMolalityActCoeff() const
6232 {
6233  double tmp;
6234  /*
6235  * Calculate the molalities. Currently, the molalities
6236  * may not be current with respect to the contents of the
6237  * State objects' data.
6238  */
6239  calcMolalities();
6240 
6241  double xmolSolvent = moleFraction(m_indexSolvent);
6242  double xx = std::max(m_xmolSolventMIN, xmolSolvent);
6243  if (IMS_typeCutoff_ == 0) {
6244  for (size_t k = 1; k < m_kk; k++) {
6245  IMS_lnActCoeffMolal_[k]= 0.0;
6246  }
6247  IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
6248  return;
6249  } else if (IMS_typeCutoff_ == 1) {
6250  if (xmolSolvent > 3.0 * IMS_X_o_cutoff_/2.0) {
6251  for (size_t k = 1; k < m_kk; k++) {
6252  IMS_lnActCoeffMolal_[k]= 0.0;
6253  }
6254  IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
6255  return;
6256  } else if (xmolSolvent < IMS_X_o_cutoff_/2.0) {
6257  tmp = log(xx * IMS_gamma_k_min_);
6258  for (size_t k = 1; k < m_kk; k++) {
6259  IMS_lnActCoeffMolal_[k]= tmp;
6260  }
6261  IMS_lnActCoeffMolal_[m_indexSolvent] = log(IMS_gamma_o_min_);
6262  return;
6263  } else {
6264  /*
6265  * If we are in the middle region, calculate the connecting polynomials
6266  */
6267  double xminus = xmolSolvent - IMS_X_o_cutoff_/2.0;
6268  double xminus2 = xminus * xminus;
6269  double xminus3 = xminus2 * xminus;
6270  double x_o_cut2 = IMS_X_o_cutoff_ * IMS_X_o_cutoff_;
6271  double x_o_cut3 = x_o_cut2 * IMS_X_o_cutoff_;
6272 
6273  double h2 = 3.5 * xminus2 / IMS_X_o_cutoff_ - 2.0 * xminus3 / x_o_cut2;
6274  double h2_prime = 7.0 * xminus / IMS_X_o_cutoff_ - 6.0 * xminus2 / x_o_cut2;
6275 
6276  double h1 = (1.0 - 3.0 * xminus2 / x_o_cut2 + 2.0 * xminus3/ x_o_cut3);
6277  double h1_prime = (- 6.0 * xminus / x_o_cut2 + 6.0 * xminus2/ x_o_cut3);
6278 
6279  double h1_g = h1 / IMS_gamma_o_min_;
6280  double h1_g_prime = h1_prime / IMS_gamma_o_min_;
6281 
6282  double alpha = 1.0 / (exp(1.0) * IMS_gamma_k_min_);
6283  double h1_f = h1 * alpha;
6284  double h1_f_prime = h1_prime * alpha;
6285 
6286  double f = h2 + h1_f;
6287  double f_prime = h2_prime + h1_f_prime;
6288 
6289  double g = h2 + h1_g;
6290  double g_prime = h2_prime + h1_g_prime;
6291 
6292  tmp = (xmolSolvent/ g * g_prime + (1.0-xmolSolvent) / f * f_prime);
6293  double lngammak = -1.0 - log(f) + tmp * xmolSolvent;
6294  double lngammao =-log(g) - tmp * (1.0-xmolSolvent);
6295 
6296  tmp = log(xmolSolvent) + lngammak;
6297  for (size_t k = 1; k < m_kk; k++) {
6298  IMS_lnActCoeffMolal_[k]= tmp;
6299  }
6300  IMS_lnActCoeffMolal_[m_indexSolvent] = lngammao;
6301  }
6302  }
6303  // Exponentials - trial 2
6304  else if (IMS_typeCutoff_ == 2) {
6305  if (xmolSolvent > IMS_X_o_cutoff_) {
6306  for (size_t k = 1; k < m_kk; k++) {
6307  IMS_lnActCoeffMolal_[k]= 0.0;
6308  }
6309  IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
6310  return;
6311  } else {
6312 
6313  double xoverc = xmolSolvent/IMS_cCut_;
6314  double eterm = std::exp(-xoverc);
6315 
6316  double fptmp = IMS_bfCut_ - IMS_afCut_ / IMS_cCut_ - IMS_bfCut_*xoverc
6317  + 2.0*IMS_dfCut_*xmolSolvent - IMS_dfCut_*xmolSolvent*xoverc;
6318  double f_prime = 1.0 + eterm*fptmp;
6319  double f = xmolSolvent + IMS_efCut_
6320  + eterm * (IMS_afCut_ + xmolSolvent * (IMS_bfCut_ + IMS_dfCut_*xmolSolvent));
6321 
6322  double gptmp = IMS_bgCut_ - IMS_agCut_ / IMS_cCut_ - IMS_bgCut_*xoverc
6323  + 2.0*IMS_dgCut_*xmolSolvent - IMS_dgCut_*xmolSolvent*xoverc;
6324  double g_prime = 1.0 + eterm*gptmp;
6325  double g = xmolSolvent + IMS_egCut_
6326  + eterm * (IMS_agCut_ + xmolSolvent * (IMS_bgCut_ + IMS_dgCut_*xmolSolvent));
6327 
6328  tmp = (xmolSolvent / g * g_prime + (1.0 - xmolSolvent) / f * f_prime);
6329  double lngammak = -1.0 - log(f) + tmp * xmolSolvent;
6330  double lngammao =-log(g) - tmp * (1.0-xmolSolvent);
6331 
6332  tmp = log(xx) + lngammak;
6333  for (size_t k = 1; k < m_kk; k++) {
6334  IMS_lnActCoeffMolal_[k]= tmp;
6335  }
6336  IMS_lnActCoeffMolal_[m_indexSolvent] = lngammao;
6337  }
6338  }
6339  return;
6340 }
6341 
6342 
6343 /**
6344  * This routine prints out the input pitzer coefficients for the
6345  * current mechanism
6346  */
6347 void HMWSoln::printCoeffs() const
6348 {
6349  size_t i, j, k;
6350  std::string sni, snj;
6351  calcMolalities();
6352  const double* charge = DATA_PTR(m_speciesCharge);
6353  double* molality = DATA_PTR(m_molalitiesCropped);
6354  double* moleF = DATA_PTR(m_tmpV);
6355  /*
6356  * Update the coefficients wrt Temperature
6357  * Calculate the derivatives as well
6358  */
6359  s_updatePitzer_CoeffWRTemp(2);
6360  getMoleFractions(moleF);
6361 
6362  printf("Index Name MoleF MolalityCropped Charge\n");
6363  for (k = 0; k < m_kk; k++) {
6364  sni = speciesName(k);
6365  printf("%2s %-16s %14.7le %14.7le %5.1f \n",
6366  int2str(k).c_str(), sni.c_str(), moleF[k], molality[k], charge[k]);
6367  }
6368 
6369  printf("\n Species Species beta0MX "
6370  "beta1MX beta2MX CphiMX alphaMX thetaij \n");
6371  for (i = 1; i < m_kk - 1; i++) {
6372  sni = speciesName(i);
6373  for (j = i+1; j < m_kk; j++) {
6374  snj = speciesName(j);
6375  size_t n = i * m_kk + j;
6376  size_t ct = m_CounterIJ[n];
6377  printf(" %-16s %-16s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f \n",
6378  sni.c_str(), snj.c_str(),
6379  m_Beta0MX_ij[ct], m_Beta1MX_ij[ct],
6380  m_Beta2MX_ij[ct], m_CphiMX_ij[ct],
6381  m_Alpha1MX_ij[ct], m_Theta_ij[ct]);
6382 
6383 
6384  }
6385  }
6386 
6387  printf("\n Species Species Species "
6388  "psi \n");
6389  for (i = 1; i < m_kk; i++) {
6390  sni = speciesName(i);
6391  for (j = 1; j < m_kk; j++) {
6392  snj = speciesName(j);
6393  for (k = 1; k < m_kk; k++) {
6394  std::string snk = speciesName(k);
6395  size_t n = k + j * m_kk + i * m_kk * m_kk;
6396  if (m_Psi_ijk[n] != 0.0) {
6397  printf(" %-16s %-16s %-16s %9.5f \n",
6398  sni.c_str(), snj.c_str(),
6399  snk.c_str(), m_Psi_ijk[n]);
6400  }
6401  }
6402  }
6403  }
6404 }
6405 
6406 //! Apply the current phScale to a set of activity Coefficients or activities
6407 /*!
6408  * See the Eq3/6 Manual for a thorough discussion.
6409  *
6410  * @param acMolality input/Output vector containing the molality based
6411  * activity coefficients. length: m_kk.
6412  */
6413 void HMWSoln::applyphScale(doublereal* acMolality) const
6414 {
6415  if (m_pHScalingType == PHSCALE_PITZER) {
6416  return;
6417  }
6418  AssertTrace(m_pHScalingType == PHSCALE_NBS);
6419  doublereal lnGammaClMs2 = s_NBS_CLM_lnMolalityActCoeff();
6420  doublereal lnGammaCLMs1 = m_lnActCoeffMolal_Unscaled[m_indexCLM];
6421  doublereal afac = -1.0 *(lnGammaClMs2 - lnGammaCLMs1);
6422  for (size_t k = 0; k < m_kk; k++) {
6423  acMolality[k] *= exp(m_speciesCharge[k] * afac);
6424  }
6425 }
6426 
6427 // Apply the current phScale to a set of activity Coefficients or activities
6428 /*
6429  * See the Eq3/6 Manual for a thorough discussion.
6430  *
6431  * @param acMolality input/Output vector containing the molality based
6432  * activity coefficients. length: m_kk.
6433  */
6434 void HMWSoln::s_updateScaling_pHScaling() const
6435 {
6436  if (m_pHScalingType == PHSCALE_PITZER) {
6437  m_lnActCoeffMolal_Scaled = m_lnActCoeffMolal_Unscaled;
6438  return;
6439  }
6440  AssertTrace(m_pHScalingType == PHSCALE_NBS);
6441  doublereal lnGammaClMs2 = s_NBS_CLM_lnMolalityActCoeff();
6442  doublereal lnGammaCLMs1 = m_lnActCoeffMolal_Unscaled[m_indexCLM];
6443  doublereal afac = -1.0 *(lnGammaClMs2 - lnGammaCLMs1);
6444  for (size_t k = 0; k < m_kk; k++) {
6445  m_lnActCoeffMolal_Scaled[k] = m_lnActCoeffMolal_Unscaled[k] + m_speciesCharge[k] * afac;
6446  }
6447 }
6448 
6449 // Apply the current phScale to a set of derivativies of the activity Coefficients
6450 // wrt temperature
6451 /*
6452  * See the Eq3/6 Manual for a thorough discussion of the need
6453  *
6454  */
6455 void HMWSoln::s_updateScaling_pHScaling_dT() const
6456 {
6457  if (m_pHScalingType == PHSCALE_PITZER) {
6458  m_dlnActCoeffMolaldT_Scaled = m_dlnActCoeffMolaldT_Unscaled;
6459  return;
6460  }
6461  AssertTrace(m_pHScalingType == PHSCALE_NBS);
6462  doublereal dlnGammaClM_dT_s2 = s_NBS_CLM_dlnMolalityActCoeff_dT();
6463  doublereal dlnGammaCLM_dT_s1 = m_dlnActCoeffMolaldT_Unscaled[m_indexCLM];
6464  doublereal afac = -1.0 *(dlnGammaClM_dT_s2 - dlnGammaCLM_dT_s1);
6465  for (size_t k = 0; k < m_kk; k++) {
6466  m_dlnActCoeffMolaldT_Scaled[k] = m_dlnActCoeffMolaldT_Unscaled[k] + m_speciesCharge[k] * afac;
6467  }
6468 }
6469 
6470 // Apply the current phScale to a set of 2nd derivatives of the activity Coefficients
6471 // wrt temperature
6472 /*
6473  * See the Eq3/6 Manual for a thorough discussion of the need
6474  *
6475  */
6476 void HMWSoln::s_updateScaling_pHScaling_dT2() const
6477 {
6478  if (m_pHScalingType == PHSCALE_PITZER) {
6479  m_d2lnActCoeffMolaldT2_Scaled = m_d2lnActCoeffMolaldT2_Unscaled;
6480  return;
6481  }
6482  AssertTrace(m_pHScalingType == PHSCALE_NBS);
6483  doublereal d2lnGammaClM_dT2_s2 = s_NBS_CLM_d2lnMolalityActCoeff_dT2();
6484  doublereal d2lnGammaCLM_dT2_s1 = m_d2lnActCoeffMolaldT2_Unscaled[m_indexCLM];
6485  doublereal afac = -1.0 *(d2lnGammaClM_dT2_s2 - d2lnGammaCLM_dT2_s1);
6486  for (size_t k = 0; k < m_kk; k++) {
6487  m_d2lnActCoeffMolaldT2_Scaled[k] = m_d2lnActCoeffMolaldT2_Unscaled[k] + m_speciesCharge[k] * afac;
6488  }
6489 }
6490 
6491 // Apply the current phScale to a set of derivatives of the activity Coefficients
6492 // wrt pressure
6493 /*
6494  * See the Eq3/6 Manual for a thorough discussion of the need
6495  */
6496 void HMWSoln::s_updateScaling_pHScaling_dP() const
6497 {
6498  if (m_pHScalingType == PHSCALE_PITZER) {
6499  m_dlnActCoeffMolaldP_Scaled = m_dlnActCoeffMolaldP_Unscaled;
6500  return;
6501  }
6502  AssertTrace(m_pHScalingType == PHSCALE_NBS);
6503  doublereal dlnGammaClM_dP_s2 = s_NBS_CLM_dlnMolalityActCoeff_dP();
6504  doublereal dlnGammaCLM_dP_s1 = m_dlnActCoeffMolaldP_Unscaled[m_indexCLM];
6505  doublereal afac = -1.0 *(dlnGammaClM_dP_s2 - dlnGammaCLM_dP_s1);
6506  for (size_t k = 0; k < m_kk; k++) {
6507  m_dlnActCoeffMolaldP_Scaled[k] = m_dlnActCoeffMolaldP_Unscaled[k] + m_speciesCharge[k] * afac;
6508  }
6509 }
6510 
6511 // Calculate the temperature derivative of the Chlorine activity coefficient
6512 /*
6513  * We assume here that the m_IionicMolality variable is up to date.
6514  */
6515 doublereal HMWSoln::s_NBS_CLM_lnMolalityActCoeff() const
6516 {
6517  doublereal sqrtIs = sqrt(m_IionicMolality);
6518  doublereal A = m_A_Debye;
6519  doublereal lnGammaClMs2 = - A * sqrtIs /(1.0 + 1.5 * sqrtIs);
6520  return lnGammaClMs2;
6521 }
6522 
6523 // Calculate the temperature derivative of the Chlorine activity coefficient
6524 /*
6525  * We assume here that the m_IionicMolality variable is up to date.
6526  */
6527 doublereal HMWSoln::s_NBS_CLM_dlnMolalityActCoeff_dT() const
6528 {
6529  doublereal sqrtIs = sqrt(m_IionicMolality);
6530  doublereal dAdT = dA_DebyedT_TP();
6531  doublereal d_lnGammaClM_dT = - dAdT * sqrtIs /(1.0 + 1.5 * sqrtIs);
6532  return d_lnGammaClM_dT;
6533 }
6534 
6535 // Calculate the second temperature derivative of the Chlorine activity coefficient
6536 /*
6537  * We assume here that the m_IionicMolality variable is up to date.
6538  */
6539 doublereal HMWSoln::s_NBS_CLM_d2lnMolalityActCoeff_dT2() const
6540 {
6541  doublereal sqrtIs = sqrt(m_IionicMolality);
6542  doublereal d2AdT2 = d2A_DebyedT2_TP();
6543  doublereal d_lnGammaClM_dT2 = - d2AdT2 * sqrtIs /(1.0 + 1.5 * sqrtIs);
6544  return d_lnGammaClM_dT2;
6545 }
6546 
6547 // Calculate the pressure derivative of the Chlorine activity coefficient
6548 /*
6549  * We assume here that the m_IionicMolality variable is up to date.
6550  */
6551 doublereal HMWSoln::s_NBS_CLM_dlnMolalityActCoeff_dP() const
6552 {
6553  doublereal sqrtIs = sqrt(m_IionicMolality);
6554  doublereal dAdP = dA_DebyedP_TP();
6555  doublereal d_lnGammaClM_dP = - dAdP * sqrtIs /(1.0 + 1.5 * sqrtIs);
6556  return d_lnGammaClM_dP;
6557 }
6558 
6559 int HMWSoln::debugPrinting()
6560 {
6561 #ifdef DEBUG_MODE
6562  return m_debugCalc;
6563 #else
6564  return 0;
6565 #endif
6566 }
6567 
6568 }
6569 /*****************************************************************************/
Generated by   doxygen 1.8.2