Cantera  2.1.2
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"
26 #include <cstdio>
27 
28 namespace Cantera
29 {
30 
33  m_formPitzer(PITZERFORM_BASE),
34  m_formPitzerTemp(PITZER_TEMP_CONSTANT),
35  m_formGC(2),
36  m_IionicMolality(0.0),
37  m_maxIionicStrength(100.0),
38  m_TempPitzerRef(298.15),
39  m_IionicMolalityStoich(0.0),
40  m_form_A_Debye(A_DEBYE_WATER),
41  m_A_Debye(1.172576), // units = sqrt(kg/gmol)
42  m_waterSS(0),
43  m_densWaterSS(1000.),
44  m_waterProps(0),
45  m_molalitiesAreCropped(false),
46  IMS_typeCutoff_(0),
47  IMS_X_o_cutoff_(0.2),
48  IMS_gamma_o_min_(1.0E-5),
49  IMS_gamma_k_min_(10.0),
50  IMS_cCut_(0.05),
51  IMS_slopefCut_(0.6),
52  IMS_slopegCut_(0.0),
53  IMS_dfCut_(0.0),
54  IMS_efCut_(0.0),
55  IMS_afCut_(0.0),
56  IMS_bfCut_(0.0),
57  IMS_dgCut_(0.0),
58  IMS_egCut_(0.0),
59  IMS_agCut_(0.0),
60  IMS_bgCut_(0.0),
61  MC_X_o_cutoff_(0.0),
62  MC_X_o_min_(0.0),
63  MC_slopepCut_(0.0),
64  MC_dpCut_(0.0),
65  MC_epCut_(0.0),
66  MC_apCut_(0.0),
67  MC_bpCut_(0.0),
68  MC_cpCut_(0.0),
69  CROP_ln_gamma_o_min(-6.0),
70  CROP_ln_gamma_o_max(3.0),
71  CROP_ln_gamma_k_min(-5.0),
72  CROP_ln_gamma_k_max(15.0),
73  m_debugCalc(0)
74 {
75  for (size_t i = 0; i < 17; i++) {
76  elambda[i] = 0.0;
77  elambda1[i] = 0.0;
78  }
79 }
80 
81 HMWSoln::HMWSoln(const std::string& inputFile, const std::string& id_) :
83  m_formPitzer(PITZERFORM_BASE),
84  m_formPitzerTemp(PITZER_TEMP_CONSTANT),
85  m_formGC(2),
86  m_IionicMolality(0.0),
87  m_maxIionicStrength(100.0),
88  m_TempPitzerRef(298.15),
89  m_IionicMolalityStoich(0.0),
90  m_form_A_Debye(A_DEBYE_WATER),
91  m_A_Debye(1.172576), // units = sqrt(kg/gmol)
92  m_waterSS(0),
93  m_densWaterSS(1000.),
94  m_waterProps(0),
95  m_molalitiesAreCropped(false),
96  IMS_typeCutoff_(0),
97  IMS_X_o_cutoff_(0.2),
98  IMS_gamma_o_min_(1.0E-5),
99  IMS_gamma_k_min_(10.0),
100  IMS_cCut_(0.05),
101  IMS_slopefCut_(0.6),
102  IMS_slopegCut_(0.0),
103  IMS_dfCut_(0.0),
104  IMS_efCut_(0.0),
105  IMS_afCut_(0.0),
106  IMS_bfCut_(0.0),
107  IMS_dgCut_(0.0),
108  IMS_egCut_(0.0),
109  IMS_agCut_(0.0),
110  IMS_bgCut_(0.0),
111  MC_X_o_cutoff_(0.0),
112  MC_X_o_min_(0.0),
113  MC_slopepCut_(0.0),
114  MC_dpCut_(0.0),
115  MC_epCut_(0.0),
116  MC_apCut_(0.0),
117  MC_bpCut_(0.0),
118  MC_cpCut_(0.0),
119  CROP_ln_gamma_o_min(-6.0),
120  CROP_ln_gamma_o_max(3.0),
121  CROP_ln_gamma_k_min(-5.0),
122  CROP_ln_gamma_k_max(15.0),
123  m_debugCalc(0)
124 {
125  for (int i = 0; i < 17; i++) {
126  elambda[i] = 0.0;
127  elambda1[i] = 0.0;
128  }
129  initThermoFile(inputFile, id_);
130 }
131 
132 HMWSoln::HMWSoln(XML_Node& phaseRoot, const std::string& id_) :
133  MolalityVPSSTP(),
134  m_formPitzer(PITZERFORM_BASE),
135  m_formPitzerTemp(PITZER_TEMP_CONSTANT),
136  m_formGC(2),
137  m_IionicMolality(0.0),
138  m_maxIionicStrength(100.0),
139  m_TempPitzerRef(298.15),
140  m_IionicMolalityStoich(0.0),
141  m_form_A_Debye(A_DEBYE_WATER),
142  m_A_Debye(1.172576), // units = sqrt(kg/gmol)
143  m_waterSS(0),
144  m_densWaterSS(1000.),
145  m_waterProps(0),
146  m_molalitiesAreCropped(false),
147  IMS_typeCutoff_(0),
148  IMS_X_o_cutoff_(0.2),
149  IMS_gamma_o_min_(1.0E-5),
150  IMS_gamma_k_min_(10.0),
151  IMS_cCut_(0.05),
152  IMS_slopefCut_(0.6),
153  IMS_slopegCut_(0.0),
154  IMS_dfCut_(0.0),
155  IMS_efCut_(0.0),
156  IMS_afCut_(0.0),
157  IMS_bfCut_(0.0),
158  IMS_dgCut_(0.0),
159  IMS_egCut_(0.0),
160  IMS_agCut_(0.0),
161  IMS_bgCut_(0.0),
162  MC_X_o_cutoff_(0.0),
163  MC_X_o_min_(0.0),
164  MC_slopepCut_(0.0),
165  MC_dpCut_(0.0),
166  MC_epCut_(0.0),
167  MC_apCut_(0.0),
168  MC_bpCut_(0.0),
169  MC_cpCut_(0.0),
170  CROP_ln_gamma_o_min(-6.0),
171  CROP_ln_gamma_o_max(3.0),
172  CROP_ln_gamma_k_min(-5.0),
173  CROP_ln_gamma_k_max(15.0),
174  m_debugCalc(0)
175 {
176  for (int i = 0; i < 17; i++) {
177  elambda[i] = 0.0;
178  elambda1[i] = 0.0;
179  }
180  importPhase(*findXMLPhase(&phaseRoot, id_), this);
181 }
182 
184  MolalityVPSSTP(),
185  m_formPitzer(PITZERFORM_BASE),
186  m_formPitzerTemp(PITZER_TEMP_CONSTANT),
187  m_formGC(2),
188  m_IionicMolality(0.0),
189  m_maxIionicStrength(100.0),
190  m_TempPitzerRef(298.15),
191  m_IionicMolalityStoich(0.0),
192  m_form_A_Debye(A_DEBYE_WATER),
193  m_A_Debye(1.172576), // units = sqrt(kg/gmol)
194  m_waterSS(0),
195  m_densWaterSS(1000.),
196  m_waterProps(0),
197  m_molalitiesAreCropped(false),
198  IMS_typeCutoff_(0),
199  IMS_X_o_cutoff_(0.2),
200  IMS_gamma_o_min_(1.0E-5),
201  IMS_gamma_k_min_(10.0),
202  IMS_cCut_(0.05),
203  IMS_slopefCut_(0.6),
204  IMS_slopegCut_(0.0),
205  IMS_dfCut_(0.0),
206  IMS_efCut_(0.0),
207  IMS_afCut_(0.0),
208  IMS_bfCut_(0.0),
209  IMS_dgCut_(0.0),
210  IMS_egCut_(0.0),
211  IMS_agCut_(0.0),
212  IMS_bgCut_(0.0),
213  MC_X_o_cutoff_(0.0),
214  MC_X_o_min_(0.0),
215  MC_slopepCut_(0.0),
216  MC_dpCut_(0.0),
217  MC_epCut_(0.0),
218  MC_apCut_(0.0),
219  MC_bpCut_(0.0),
220  MC_cpCut_(0.0),
221  CROP_ln_gamma_o_min(-6.0),
222  CROP_ln_gamma_o_max(3.0),
223  CROP_ln_gamma_k_min(-5.0),
224  CROP_ln_gamma_k_max(15.0),
225  m_debugCalc(0)
226 {
227  /*
228  * Use the assignment operator to do the brunt
229  * of the work for the copy constructor.
230  */
231  *this = b;
232 }
233 
236 {
237  if (&b != this) {
241  m_formGC = b.m_formGC;
242  m_Aionic = b.m_Aionic;
248  m_A_Debye = b.m_A_Debye;
249 
250  // This is an internal shallow copy of the PDSS_Water pointer
251  m_waterSS = providePDSS(0);
252  if (!m_waterSS) {
253  throw CanteraError("HMWSoln::operator=()", "Dynamic cast to PDSS_Water failed");
254  }
255 
257 
258  if (m_waterProps) {
259  delete m_waterProps;
260  m_waterProps = 0;
261  }
262  if (b.m_waterProps) {
263  m_waterProps = new WaterProps(dynamic_cast<PDSS_Water*>(m_waterSS));
264  }
265 
266  m_pp = b.m_pp;
267  m_tmpV = b.m_tmpV;
296  m_Psi_ijk = b.m_Psi_ijk;
306 
307  m_Mu_nnn = b.m_Mu_nnn;
312 
315 
318 
321 
324 
328 
333  m_BMX_IJ = b.m_BMX_IJ;
345  m_Phi_IJ = b.m_Phi_IJ;
354  m_CMX_IJ = b.m_CMX_IJ;
359 
365  IMS_cCut_ = b.IMS_cCut_;
367  IMS_dfCut_ = b.IMS_dfCut_;
368  IMS_efCut_ = b.IMS_efCut_;
369  IMS_afCut_ = b.IMS_afCut_;
370  IMS_bfCut_ = b.IMS_bfCut_;
372  IMS_dgCut_ = b.IMS_dgCut_;
373  IMS_egCut_ = b.IMS_egCut_;
374  IMS_agCut_ = b.IMS_agCut_;
375  IMS_bgCut_ = b.IMS_bgCut_;
379  MC_dpCut_ = b.MC_dpCut_;
380  MC_epCut_ = b.MC_epCut_;
381  MC_apCut_ = b.MC_apCut_;
382  MC_bpCut_ = b.MC_bpCut_;
383  MC_cpCut_ = b.MC_cpCut_;
384  CROP_ln_gamma_o_min = b.CROP_ln_gamma_o_min;
385  CROP_ln_gamma_o_max = b.CROP_ln_gamma_o_max;
386  CROP_ln_gamma_k_min = b.CROP_ln_gamma_k_min;
387  CROP_ln_gamma_k_max = b.CROP_ln_gamma_k_max;
389 
391  }
392  return *this;
393 }
394 
395 HMWSoln::HMWSoln(int testProb) :
396  MolalityVPSSTP(),
397  m_formPitzer(PITZERFORM_BASE),
398  m_formPitzerTemp(PITZER_TEMP_CONSTANT),
399  m_formGC(2),
400  m_IionicMolality(0.0),
401  m_maxIionicStrength(100.0),
402  m_TempPitzerRef(298.15),
403  m_IionicMolalityStoich(0.0),
404  m_form_A_Debye(A_DEBYE_WATER),
405  m_A_Debye(1.172576), // units = sqrt(kg/gmol)
406  m_waterSS(0),
407  m_densWaterSS(1000.),
408  m_waterProps(0),
409  m_molalitiesAreCropped(false),
410  IMS_typeCutoff_(0),
411  IMS_X_o_cutoff_(0.2),
412  IMS_gamma_o_min_(1.0E-5),
413  IMS_gamma_k_min_(10.0),
414  IMS_cCut_(0.05),
415  IMS_slopefCut_(0.6),
416  IMS_slopegCut_(0.0),
417  IMS_dfCut_(0.0),
418  IMS_efCut_(0.0),
419  IMS_afCut_(0.0),
420  IMS_bfCut_(0.0),
421  IMS_dgCut_(0.0),
422  IMS_egCut_(0.0),
423  IMS_agCut_(0.0),
424  IMS_bgCut_(0.0),
425  MC_X_o_cutoff_(0.0),
426  MC_X_o_min_(0.0),
427  MC_slopepCut_(0.0),
428  MC_dpCut_(0.0),
429  MC_epCut_(0.0),
430  MC_apCut_(0.0),
431  MC_bpCut_(0.0),
432  MC_cpCut_(0.0),
433  CROP_ln_gamma_o_min(-6.0),
434  CROP_ln_gamma_o_max(3.0),
435  CROP_ln_gamma_k_min(-5.0),
436  CROP_ln_gamma_k_max(15.0),
437  m_debugCalc(0)
438 {
439  if (testProb != 1) {
440  printf("unknown test problem\n");
441  exit(EXIT_FAILURE);
442  }
443 
444  initThermoFile("HMW_NaCl.xml", "");
445 
446  size_t i = speciesIndex("Cl-");
447  size_t j = speciesIndex("H+");
448  size_t n = i * m_kk + j;
449  size_t ct = m_CounterIJ[n];
450  m_Beta0MX_ij[ct] = 0.1775;
451  m_Beta1MX_ij[ct] = 0.2945;
452  m_CphiMX_ij[ct] = 0.0008;
453  m_Alpha1MX_ij[ct]= 2.000;
454 
455 
456  i = speciesIndex("Cl-");
457  j = speciesIndex("Na+");
458  n = i * m_kk + j;
459  ct = m_CounterIJ[n];
460  m_Beta0MX_ij[ct] = 0.0765;
461  m_Beta1MX_ij[ct] = 0.2664;
462  m_CphiMX_ij[ct] = 0.00127;
463  m_Alpha1MX_ij[ct]= 2.000;
464 
465 
466  i = speciesIndex("Cl-");
467  j = speciesIndex("OH-");
468  n = i * m_kk + j;
469  ct = m_CounterIJ[n];
470  m_Theta_ij[ct] = -0.05;
471 
472  i = speciesIndex("H+");
473  j = speciesIndex("Na+");
474  n = i * m_kk + j;
475  ct = m_CounterIJ[n];
476  m_Theta_ij[ct] = 0.036;
477 
478  i = speciesIndex("Na+");
479  j = speciesIndex("OH-");
480  n = i * m_kk + j;
481  ct = m_CounterIJ[n];
482  m_Beta0MX_ij[ct] = 0.0864;
483  m_Beta1MX_ij[ct] = 0.253;
484  m_CphiMX_ij[ct] = 0.0044;
485  m_Alpha1MX_ij[ct]= 2.000;
486 
487  i = speciesIndex("Cl-");
488  j = speciesIndex("H+");
489  size_t k = speciesIndex("Na+");
490  double param = -0.004;
491  n = i * m_kk *m_kk + j * m_kk + k ;
492  m_Psi_ijk[n] = param;
493  m_Psi_ijk_coeff(0,n) = param;
494  n = i * m_kk *m_kk + k * m_kk + j ;
495  m_Psi_ijk[n] = param;
496  m_Psi_ijk_coeff(0,n) = param;
497  n = j * m_kk *m_kk + i * m_kk + k ;
498  m_Psi_ijk[n] = param;
499  m_Psi_ijk_coeff(0,n) = param;
500  n = j * m_kk *m_kk + k * m_kk + i ;
501  m_Psi_ijk[n] = param;
502  m_Psi_ijk_coeff(0,n) = param;
503  n = k * m_kk *m_kk + j * m_kk + i ;
504  m_Psi_ijk[n] = param;
505  m_Psi_ijk_coeff(0,n) = param;
506  n = k * m_kk *m_kk + i * m_kk + j ;
507  m_Psi_ijk[n] = param;
508  m_Psi_ijk_coeff(0,n) = param;
509 
510  i = speciesIndex("Cl-");
511  j = speciesIndex("Na+");
512  k = speciesIndex("OH-");
513  param = -0.006;
514  n = i * m_kk *m_kk + j * m_kk + k ;
515  m_Psi_ijk[n] = param;
516  m_Psi_ijk_coeff(0,n) = param;
517  n = i * m_kk *m_kk + k * m_kk + j ;
518  m_Psi_ijk[n] = param;
519  m_Psi_ijk_coeff(0,n) = param;
520  n = j * m_kk *m_kk + i * m_kk + k ;
521  m_Psi_ijk[n] = param;
522  m_Psi_ijk_coeff(0,n) = param;
523  n = j * m_kk *m_kk + k * m_kk + i ;
524  m_Psi_ijk[n] = param;
525  m_Psi_ijk_coeff(0,n) = param;
526  n = k * m_kk *m_kk + j * m_kk + i ;
527  m_Psi_ijk[n] = param;
528  m_Psi_ijk_coeff(0,n) = param;
529  n = k * m_kk *m_kk + i * m_kk + j ;
530  m_Psi_ijk[n] = param;
531  m_Psi_ijk_coeff(0,n) = param;
532 
533  printCoeffs();
534 }
535 
537 {
538  if (m_waterProps) {
539  delete m_waterProps;
540  m_waterProps = 0;
541  }
542 }
543 
545 {
546  return new HMWSoln(*this);
547 }
548 
549 int HMWSoln::eosType() const
550 {
551  int res;
552  switch (m_formGC) {
553  case 0:
554  res = cHMWSoln0;
555  break;
556  case 1:
557  res = cHMWSoln1;
558  break;
559  case 2:
560  res = cHMWSoln2;
561  break;
562  default:
563  throw CanteraError("eosType", "Unknown type");
564  }
565  return res;
566 }
567 
568 //
569 // -------- Molar Thermodynamic Properties of the Solution ---------------
570 //
571 doublereal HMWSoln::enthalpy_mole() const
572 {
575  return mean_X(DATA_PTR(m_tmpV));
576 }
577 
578 doublereal HMWSoln::relative_enthalpy() const
579 {
581  double hbar = mean_X(DATA_PTR(m_tmpV));
583  double RT = GasConstant * temperature();
584  for (size_t k = 0; k < m_kk; k++) {
585  m_gamma_tmp[k] *= RT;
586  }
587  double h0bar = mean_X(DATA_PTR(m_gamma_tmp));
588  return hbar - h0bar;
589 }
590 
592 {
593  double L = relative_enthalpy();
595  double xanion = 0.0;
596  size_t kcation = npos;
597  double xcation = 0.0;
598  size_t kanion = npos;
599  for (size_t k = 0; k < m_kk; k++) {
600  if (charge(k) > 0.0) {
601  if (m_tmpV[k] > xanion) {
602  xanion = m_tmpV[k];
603  kanion = k;
604  }
605  } else if (charge(k) < 0.0) {
606  if (m_tmpV[k] > xcation) {
607  xcation = m_tmpV[k];
608  kcation = k;
609  }
610  }
611  }
612  if (kcation == npos || kanion == npos) {
613  return L;
614  }
615  double xuse = xcation;
616  double factor = 1;
617  if (xanion < xcation) {
618  xuse = xanion;
619  if (charge(kcation) != 1.0) {
620  factor = charge(kcation);
621  }
622  } else {
623  if (charge(kanion) != 1.0) {
624  factor = charge(kanion);
625  }
626  }
627  xuse = xuse / factor;
628  return L / xuse;
629 }
630 
631 doublereal HMWSoln::intEnergy_mole() const
632 {
633  double hh = enthalpy_mole();
634  double pres = pressure();
635  double molarV = 1.0/molarDensity();
636  return hh - pres * molarV;
637 }
638 
639 doublereal HMWSoln::entropy_mole() const
640 {
642  return mean_X(DATA_PTR(m_tmpV));
643 }
644 
645 doublereal HMWSoln::gibbs_mole() const
646 {
648  return mean_X(DATA_PTR(m_tmpV));
649 }
650 
651 doublereal HMWSoln::cp_mole() const
652 {
654  return mean_X(DATA_PTR(m_tmpV));
655 }
656 
657 doublereal HMWSoln::cv_mole() const
658 {
659  double kappa_t = isothermalCompressibility();
660  double beta = thermalExpansionCoeff();
661  double cp = cp_mole();
662  double tt = temperature();
663  double molarV = molarVolume();
664  return cp - beta * beta * tt * molarV / kappa_t;
665 }
666 
667 //
668 // ------- Mechanical Equation of State Properties ------------------------
669 //
670 
671 doublereal HMWSoln::pressure() const
672 {
673  return m_Pcurrent;
674 }
675 
676 void HMWSoln::setPressure(doublereal p)
677 {
678  setState_TP(temperature(), p);
679 }
680 
682 {
683  double* vbar = &m_pp[0];
685  double* x = &m_tmpV[0];
686  getMoleFractions(x);
687  doublereal vtotal = 0.0;
688  for (size_t i = 0; i < m_kk; i++) {
689  vtotal += vbar[i] * x[i];
690  }
691  doublereal dd = meanMolecularWeight() / vtotal;
692  Phase::setDensity(dd);
693 }
694 
696 {
697  throw CanteraError("HMWSoln::isothermalCompressibility",
698  "unimplemented");
699  return 0.0;
700 }
701 
703 {
704  throw CanteraError("HMWSoln::thermalExpansionCoeff",
705  "unimplemented");
706  return 0.0;
707 }
708 
709 double HMWSoln::density() const
710 {
711  // calcDensity();
712  return Phase::density();
713 }
714 
715 void HMWSoln::setDensity(const doublereal rho)
716 {
717  double dens_old = density();
718 
719  if (rho != dens_old) {
720  throw CanteraError("HMWSoln::setDensity",
721  "Density is not an independent variable");
722  }
723 }
724 
725 void HMWSoln::setMolarDensity(const doublereal rho)
726 {
727  throw CanteraError("HMWSoln::setMolarDensity",
728  "Density is not an independent variable");
729 }
730 
731 void HMWSoln::setTemperature(const doublereal temp)
732 {
733  setState_TP(temp, m_Pcurrent);
734 }
735 
736 void HMWSoln::setState_TP(doublereal temp, doublereal pres)
737 {
738  Phase::setTemperature(temp);
739  /*
740  * Store the current pressure
741  */
742  m_Pcurrent = pres;
743 
744  /*
745  * update the standard state thermo
746  * -> This involves calling the water function and setting the pressure
747  */
749  /*
750  * Store the internal density of the water SS.
751  * Note, we would have to do this for all other
752  * species if they had pressure dependent properties.
753  */
755  /*
756  * Calculate all of the other standard volumes
757  * -> note these are constant for now
758  */
759  calcDensity();
760 }
761 
762 //
763 // ------- Activities and Activity Concentrations
764 //
765 
766 void HMWSoln::getActivityConcentrations(doublereal* c) const
767 {
768  double cs_solvent = standardConcentration();
769  getActivities(c);
770  c[0] *= cs_solvent;
771  if (m_kk > 1) {
772  double cs_solute = standardConcentration(1);
773  for (size_t k = 1; k < m_kk; k++) {
774  c[k] *= cs_solute;
775  }
776  }
777 }
778 
779 doublereal HMWSoln::standardConcentration(size_t k) const
780 {
782  double mvSolvent = m_tmpV[m_indexSolvent];
783  if (k > 0) {
784  return m_Mnaught / mvSolvent;
785  }
786  return 1.0 / mvSolvent;
787 }
788 
789 doublereal HMWSoln::logStandardConc(size_t k) const
790 {
791  double c_solvent = standardConcentration(k);
792  return log(c_solvent);
793 }
794 
795 void HMWSoln::getUnitsStandardConc(double* uA, int k, int sizeUA) const
796 {
797  for (int i = 0; i < sizeUA; i++) {
798  if (i == 0) {
799  uA[0] = 1.0;
800  }
801  if (i == 1) {
802  uA[1] = -int(nDim());
803  }
804  if (i == 2) {
805  uA[2] = 0.0;
806  }
807  if (i == 3) {
808  uA[3] = 0.0;
809  }
810  if (i == 4) {
811  uA[4] = 0.0;
812  }
813  if (i == 5) {
814  uA[5] = 0.0;
815  }
816  }
817 }
818 
819 void HMWSoln::getActivities(doublereal* ac) const
820 {
822  /*
823  * Update the molality array, m_molalities()
824  * This requires an update due to mole fractions
825  */
826  s_update_lnMolalityActCoeff();
827  /*
828  * Now calculate the array of activities.
829  */
830  for (size_t k = 0; k < m_kk; k++) {
831  if (k != m_indexSolvent) {
832  ac[k] = m_molalities[k] * exp(m_lnActCoeffMolal_Scaled[k]);
833  }
834  }
835  double xmolSolvent = moleFraction(m_indexSolvent);
836  ac[m_indexSolvent] =
837  exp(m_lnActCoeffMolal_Scaled[m_indexSolvent]) * xmolSolvent;
838  /*
839  * Apply the pH scale
840  */
841  //applyphScale(ac);
842 }
843 
844 void HMWSoln::
845 getUnscaledMolalityActivityCoefficients(doublereal* acMolality) const
846 {
848  A_Debye_TP(-1.0, -1.0);
849  s_update_lnMolalityActCoeff();
850  std::copy(m_lnActCoeffMolal_Unscaled.begin(), m_lnActCoeffMolal_Unscaled.end(), acMolality);
851  for (size_t k = 0; k < m_kk; k++) {
852  acMolality[k] = exp(acMolality[k]);
853  }
854 }
855 
856 //
857 // ------ Partial Molar Properties of the Solution -----------------
858 //
859 
860 void HMWSoln::getChemPotentials(doublereal* mu) const
861 {
862  double xx;
863  /*
864  * First get the standard chemical potentials in
865  * molar form.
866  * -> this requires updates of standard state as a function
867  * of T and P
868  */
870  /*
871  * Update the activity coefficients
872  * This also updates the internal molality array.
873  */
874  s_update_lnMolalityActCoeff();
875  doublereal RT = GasConstant * temperature();
876  double xmolSolvent = moleFraction(m_indexSolvent);
877  for (size_t k = 0; k < m_kk; k++) {
878  if (m_indexSolvent != k) {
879  xx = std::max(m_molalities[k], SmallNumber);
880  mu[k] += RT * (log(xx) + m_lnActCoeffMolal_Scaled[k]);
881  }
882  }
883  xx = std::max(xmolSolvent, SmallNumber);
884  mu[m_indexSolvent] +=
885  RT * (log(xx) + m_lnActCoeffMolal_Scaled[m_indexSolvent]);
886 }
887 
888 void HMWSoln::getPartialMolarEnthalpies(doublereal* hbar) const
889 {
890  /*
891  * Get the nondimensional standard state enthalpies
892  */
893  getEnthalpy_RT(hbar);
894  /*
895  * dimensionalize it.
896  */
897  double T = temperature();
898  double RT = GasConstant * T;
899  for (size_t k = 0; k < m_kk; k++) {
900  hbar[k] *= RT;
901  }
902  /*
903  * Update the activity coefficients, This also update the
904  * internally stored molalities.
905  */
906  s_update_lnMolalityActCoeff();
908  double RTT = RT * T;
909  for (size_t k = 0; k < m_kk; k++) {
910  hbar[k] -= RTT * m_dlnActCoeffMolaldT_Scaled[k];
911  }
912 }
913 
914 void HMWSoln::
915 getPartialMolarEntropies(doublereal* sbar) const
916 {
917  /*
918  * Get the standard state entropies at the temperature
919  * and pressure of the solution.
920  */
921  getEntropy_R(sbar);
922  /*
923  * Dimensionalize the entropies
924  */
925  doublereal R = GasConstant;
926  for (size_t k = 0; k < m_kk; k++) {
927  sbar[k] *= R;
928  }
929  /*
930  * Update the activity coefficients, This also update the
931  * internally stored molalities.
932  */
933  s_update_lnMolalityActCoeff();
934  /*
935  * First we will add in the obvious dependence on the T
936  * term out front of the log activity term
937  */
938  doublereal mm;
939  for (size_t k = 0; k < m_kk; k++) {
940  if (k != m_indexSolvent) {
941  mm = std::max(SmallNumber, m_molalities[k]);
942  sbar[k] -= R * (log(mm) + m_lnActCoeffMolal_Scaled[k]);
943  }
944  }
945  double xmolSolvent = moleFraction(m_indexSolvent);
946  mm = std::max(SmallNumber, xmolSolvent);
947  sbar[m_indexSolvent] -= R *(log(mm) + m_lnActCoeffMolal_Scaled[m_indexSolvent]);
948  /*
949  * Check to see whether activity coefficients are temperature
950  * dependent. If they are, then calculate the their temperature
951  * derivatives and add them into the result.
952  */
954  double RT = R * temperature();
955  for (size_t k = 0; k < m_kk; k++) {
956  sbar[k] -= RT * m_dlnActCoeffMolaldT_Scaled[k];
957  }
958 }
959 
960 void HMWSoln::getPartialMolarVolumes(doublereal* vbar) const
961 {
962  /*
963  * Get the standard state values in m^3 kmol-1
964  */
965  getStandardVolumes(vbar);
966  /*
967  * Update the derivatives wrt the activity coefficients.
968  */
969  s_update_lnMolalityActCoeff();
971  double T = temperature();
972  double RT = GasConstant * T;
973  for (size_t k = 0; k < m_kk; k++) {
974  vbar[k] += RT * m_dlnActCoeffMolaldP_Scaled[k];
975  }
976 }
977 
978 void HMWSoln::getPartialMolarCp(doublereal* cpbar) const
979 {
980  /*
981  * Get the nondimensional gibbs standard state of the
982  * species at the T and P of the solution.
983  */
984  getCp_R(cpbar);
985 
986  for (size_t k = 0; k < m_kk; k++) {
987  cpbar[k] *= GasConstant;
988  }
989  /*
990  * Update the activity coefficients, This also update the
991  * internally stored molalities.
992  */
993  s_update_lnMolalityActCoeff();
996  double T = temperature();
997  double RT = GasConstant * T;
998  double RTT = RT * T;
999  for (size_t k = 0; k < m_kk; k++) {
1000  cpbar[k] -= (2.0 * RT * m_dlnActCoeffMolaldT_Scaled[k] +
1002  }
1003 }
1004 
1005 /*
1006  * -------------- Utilities -------------------------------
1007  */
1008 
1009 void HMWSoln::setParameters(int n, doublereal* const c)
1010 {
1011  warn_deprecated("HMWSoln::setParameters");
1012 }
1013 
1014 void HMWSoln::getParameters(int& n, doublereal* const c) const
1015 {
1016  warn_deprecated("HMWSoln::getParameters");
1017 }
1018 
1020 {
1021 }
1022 
1023 doublereal HMWSoln::satPressure(doublereal t) {
1024  double p_old = pressure();
1025  double t_old = temperature();
1026  double pres = m_waterSS->satPressure(t);
1027  /*
1028  * Set the underlying object back to its original state.
1029  */
1030  m_waterSS->setState_TP(t_old, p_old);
1031  return pres;
1032 }
1033 
1034 double HMWSoln::A_Debye_TP(double tempArg, double presArg) const
1035 {
1036  double T = temperature();
1037  double A;
1038  if (tempArg != -1.0) {
1039  T = tempArg;
1040  }
1041  double P = pressure();
1042  if (presArg != -1.0) {
1043  P = presArg;
1044  }
1045 
1046  switch (m_form_A_Debye) {
1047  case A_DEBYE_CONST:
1048  A = m_A_Debye;
1049  break;
1050  case A_DEBYE_WATER:
1051  A = m_waterProps->ADebye(T, P, 0);
1052  m_A_Debye = A;
1053  break;
1054  default:
1055  printf("shouldn't be here\n");
1056  exit(EXIT_FAILURE);
1057  }
1058  return A;
1059 }
1060 
1061 double HMWSoln::dA_DebyedT_TP(double tempArg, double presArg) const
1062 {
1063  doublereal T = temperature();
1064  if (tempArg != -1.0) {
1065  T = tempArg;
1066  }
1067  doublereal P = pressure();
1068  if (presArg != -1.0) {
1069  P = presArg;
1070  }
1071  doublereal dAdT;
1072  switch (m_form_A_Debye) {
1073  case A_DEBYE_CONST:
1074  dAdT = 0.0;
1075  break;
1076  case A_DEBYE_WATER:
1077  dAdT = m_waterProps->ADebye(T, P, 1);
1078  //dAdT = WaterProps::ADebye(T, P, 1);
1079  break;
1080  default:
1081  printf("shouldn't be here\n");
1082  exit(EXIT_FAILURE);
1083  }
1084  return dAdT;
1085 }
1086 
1087 double HMWSoln::dA_DebyedP_TP(double tempArg, double presArg) const
1088 {
1089  double T = temperature();
1090  if (tempArg != -1.0) {
1091  T = tempArg;
1092  }
1093  double P = pressure();
1094  if (presArg != -1.0) {
1095  P = presArg;
1096  }
1097  double dAdP;
1098  switch (m_form_A_Debye) {
1099  case A_DEBYE_CONST:
1100  dAdP = 0.0;
1101  break;
1102  case A_DEBYE_WATER:
1103  dAdP = m_waterProps->ADebye(T, P, 3);
1104  break;
1105  default:
1106  printf("shouldn't be here\n");
1107  exit(EXIT_FAILURE);
1108  }
1109  return dAdP;
1110 }
1111 
1112 double HMWSoln::ADebye_L(double tempArg, double presArg) const
1113 {
1114  double dAdT = dA_DebyedT_TP();
1115  double dAphidT = dAdT /3.0;
1116  double T = temperature();
1117  if (tempArg != -1.0) {
1118  T = tempArg;
1119  }
1120  return dAphidT * (4.0 * GasConstant * T * T);
1121 }
1122 
1123 double HMWSoln::ADebye_V(double tempArg, double presArg) const
1124 {
1125  double dAdP = dA_DebyedP_TP();
1126  double dAphidP = dAdP /3.0;
1127  double T = temperature();
1128  if (tempArg != -1.0) {
1129  T = tempArg;
1130  }
1131  return - dAphidP * (4.0 * GasConstant * T);
1132 }
1133 
1134 double HMWSoln::ADebye_J(double tempArg, double presArg) const
1135 {
1136  double T = temperature();
1137  if (tempArg != -1.0) {
1138  T = tempArg;
1139  }
1140  double A_L = ADebye_L(T, presArg);
1141  double d2 = d2A_DebyedT2_TP(T, presArg);
1142  double d2Aphi = d2 / 3.0;
1143  return 2.0 * A_L / T + 4.0 * GasConstant * T * T *d2Aphi;
1144 }
1145 
1146 double HMWSoln::d2A_DebyedT2_TP(double tempArg, double presArg) const
1147 {
1148  double T = temperature();
1149  if (tempArg != -1.0) {
1150  T = tempArg;
1151  }
1152  double P = pressure();
1153  if (presArg != -1.0) {
1154  P = presArg;
1155  }
1156  double d2AdT2;
1157  switch (m_form_A_Debye) {
1158  case A_DEBYE_CONST:
1159  d2AdT2 = 0.0;
1160  break;
1161  case A_DEBYE_WATER:
1162  d2AdT2 = m_waterProps->ADebye(T, P, 2);
1163  break;
1164  default:
1165  printf("shouldn't be here\n");
1166  exit(EXIT_FAILURE);
1167  }
1168  return d2AdT2;
1169 }
1170 
1171 /*
1172  * ---------- Other Property Functions
1173  */
1174 double HMWSoln::AionicRadius(int k) const
1175 {
1176  return m_Aionic[k];
1177 }
1178 
1179 /*
1180  * ------------ Private and Restricted Functions ------------------
1181  */
1182 
1183 doublereal HMWSoln::err(const std::string& msg) const
1184 {
1185  throw CanteraError("HMWSoln",
1186  "Unfinished func called: " + msg);
1187  return 0.0;
1188 }
1189 
1191 {
1192  m_kk = nSpecies();
1193 
1194  /*
1195  * Resize lengths equal to the number of species in
1196  * the phase.
1197  */
1198  m_electrolyteSpeciesType.resize(m_kk, cEST_polarNeutral);
1199  m_speciesSize.resize(m_kk);
1200  m_speciesCharge_Stoich.resize(m_kk, 0.0);
1201  m_Aionic.resize(m_kk, 0.0);
1202 
1203  m_pp.resize(m_kk, 0.0);
1204  m_tmpV.resize(m_kk, 0.0);
1205  m_molalitiesCropped.resize(m_kk, 0.0);
1206 
1207  size_t maxCounterIJlen = 1 + (m_kk-1) * (m_kk-2) / 2;
1208 
1209  /*
1210  * Figure out the size of the temperature coefficient
1211  * arrays
1212  */
1213  int TCoeffLength = 1;
1214  if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
1215  TCoeffLength = 2;
1216  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
1217  TCoeffLength = 5;
1218  }
1219 
1220  m_Beta0MX_ij.resize(maxCounterIJlen, 0.0);
1221  m_Beta0MX_ij_L.resize(maxCounterIJlen, 0.0);
1222  m_Beta0MX_ij_LL.resize(maxCounterIJlen, 0.0);
1223  m_Beta0MX_ij_P.resize(maxCounterIJlen, 0.0);
1224  m_Beta0MX_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
1225 
1226  m_Beta1MX_ij.resize(maxCounterIJlen, 0.0);
1227  m_Beta1MX_ij_L.resize(maxCounterIJlen, 0.0);
1228  m_Beta1MX_ij_LL.resize(maxCounterIJlen, 0.0);
1229  m_Beta1MX_ij_P.resize(maxCounterIJlen, 0.0);
1230  m_Beta1MX_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
1231 
1232  m_Beta2MX_ij.resize(maxCounterIJlen, 0.0);
1233  m_Beta2MX_ij_L.resize(maxCounterIJlen, 0.0);
1234  m_Beta2MX_ij_LL.resize(maxCounterIJlen, 0.0);
1235  m_Beta2MX_ij_P.resize(maxCounterIJlen, 0.0);
1236  m_Beta2MX_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
1237 
1238  m_CphiMX_ij.resize(maxCounterIJlen, 0.0);
1239  m_CphiMX_ij_L.resize(maxCounterIJlen, 0.0);
1240  m_CphiMX_ij_LL.resize(maxCounterIJlen, 0.0);
1241  m_CphiMX_ij_P.resize(maxCounterIJlen, 0.0);
1242  m_CphiMX_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
1243 
1244  m_Alpha1MX_ij.resize(maxCounterIJlen, 2.0);
1245  m_Alpha2MX_ij.resize(maxCounterIJlen, 12.0);
1246  m_Theta_ij.resize(maxCounterIJlen, 0.0);
1247  m_Theta_ij_L.resize(maxCounterIJlen, 0.0);
1248  m_Theta_ij_LL.resize(maxCounterIJlen, 0.0);
1249  m_Theta_ij_P.resize(maxCounterIJlen, 0.0);
1250  m_Theta_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
1251 
1252  size_t n = m_kk*m_kk*m_kk;
1253  m_Psi_ijk.resize(n, 0.0);
1254  m_Psi_ijk_L.resize(n, 0.0);
1255  m_Psi_ijk_LL.resize(n, 0.0);
1256  m_Psi_ijk_P.resize(n, 0.0);
1257  m_Psi_ijk_coeff.resize(TCoeffLength, n, 0.0);
1258 
1259  m_Lambda_nj.resize(m_kk, m_kk, 0.0);
1260  m_Lambda_nj_L.resize(m_kk, m_kk, 0.0);
1261  m_Lambda_nj_LL.resize(m_kk, m_kk, 0.0);
1262  m_Lambda_nj_P.resize(m_kk, m_kk, 0.0);
1263  m_Lambda_nj_coeff.resize(TCoeffLength, m_kk * m_kk, 0.0);
1264 
1265  m_Mu_nnn.resize(m_kk, 0.0);
1266  m_Mu_nnn_L.resize(m_kk, 0.0);
1267  m_Mu_nnn_LL.resize(m_kk, 0.0);
1268  m_Mu_nnn_P.resize(m_kk, 0.0);
1269  m_Mu_nnn_coeff.resize(TCoeffLength, m_kk, 0.0);
1270 
1271  m_lnActCoeffMolal_Scaled.resize(m_kk, 0.0);
1272  m_dlnActCoeffMolaldT_Scaled.resize(m_kk, 0.0);
1273  m_d2lnActCoeffMolaldT2_Scaled.resize(m_kk, 0.0);
1274  m_dlnActCoeffMolaldP_Scaled.resize(m_kk, 0.0);
1275 
1276  m_lnActCoeffMolal_Unscaled.resize(m_kk, 0.0);
1277  m_dlnActCoeffMolaldT_Unscaled.resize(m_kk, 0.0);
1278  m_d2lnActCoeffMolaldT2_Unscaled.resize(m_kk, 0.0);
1279  m_dlnActCoeffMolaldP_Unscaled.resize(m_kk, 0.0);
1280 
1281  m_CounterIJ.resize(m_kk*m_kk, 0);
1282 
1283  m_gfunc_IJ.resize(maxCounterIJlen, 0.0);
1284  m_g2func_IJ.resize(maxCounterIJlen, 0.0);
1285  m_hfunc_IJ.resize(maxCounterIJlen, 0.0);
1286  m_h2func_IJ.resize(maxCounterIJlen, 0.0);
1287  m_BMX_IJ.resize(maxCounterIJlen, 0.0);
1288  m_BMX_IJ_L.resize(maxCounterIJlen, 0.0);
1289  m_BMX_IJ_LL.resize(maxCounterIJlen, 0.0);
1290  m_BMX_IJ_P.resize(maxCounterIJlen, 0.0);
1291  m_BprimeMX_IJ.resize(maxCounterIJlen, 0.0);
1292  m_BprimeMX_IJ_L.resize(maxCounterIJlen, 0.0);
1293  m_BprimeMX_IJ_LL.resize(maxCounterIJlen, 0.0);
1294  m_BprimeMX_IJ_P.resize(maxCounterIJlen, 0.0);
1295  m_BphiMX_IJ.resize(maxCounterIJlen, 0.0);
1296  m_BphiMX_IJ_L.resize(maxCounterIJlen, 0.0);
1297  m_BphiMX_IJ_LL.resize(maxCounterIJlen, 0.0);
1298  m_BphiMX_IJ_P.resize(maxCounterIJlen, 0.0);
1299  m_Phi_IJ.resize(maxCounterIJlen, 0.0);
1300  m_Phi_IJ_L.resize(maxCounterIJlen, 0.0);
1301  m_Phi_IJ_LL.resize(maxCounterIJlen, 0.0);
1302  m_Phi_IJ_P.resize(maxCounterIJlen, 0.0);
1303  m_Phiprime_IJ.resize(maxCounterIJlen, 0.0);
1304  m_PhiPhi_IJ.resize(maxCounterIJlen, 0.0);
1305  m_PhiPhi_IJ_L.resize(maxCounterIJlen, 0.0);
1306  m_PhiPhi_IJ_LL.resize(maxCounterIJlen, 0.0);
1307  m_PhiPhi_IJ_P.resize(maxCounterIJlen, 0.0);
1308  m_CMX_IJ.resize(maxCounterIJlen, 0.0);
1309  m_CMX_IJ_L.resize(maxCounterIJlen, 0.0);
1310  m_CMX_IJ_LL.resize(maxCounterIJlen, 0.0);
1311  m_CMX_IJ_P.resize(maxCounterIJlen, 0.0);
1312 
1313  m_gamma_tmp.resize(m_kk, 0.0);
1314 
1315  IMS_lnActCoeffMolal_.resize(m_kk, 0.0);
1316  CROP_speciesCropped_.resize(m_kk, 0);
1317 
1318  counterIJ_setup();
1319 }
1320 
1321 void HMWSoln::s_update_lnMolalityActCoeff() const
1322 {
1323  /*
1324  * Calculate the molalities. Currently, the molalities
1325  * may not be current with respect to the contents of the
1326  * State objects' data.
1327  */
1328  calcMolalities();
1329  /*
1330  * Calculate a cropped set of molalities that will be used
1331  * in all activity coefficient calculations.
1332  */
1334  /*
1335  * Calculate the stoichiometric ionic charge. This isn't used in the
1336  * Pitzer formulation.
1337  */
1338  m_IionicMolalityStoich = 0.0;
1339  for (size_t k = 0; k < m_kk; k++) {
1340  double z_k = charge(k);
1341  double zs_k1 = m_speciesCharge_Stoich[k];
1342  if (z_k == zs_k1) {
1343  m_IionicMolalityStoich += m_molalities[k] * z_k * z_k;
1344  } else {
1345  double zs_k2 = z_k - zs_k1;
1347  += m_molalities[k] * (zs_k1 * zs_k1 + zs_k2 * zs_k2);
1348  }
1349  }
1350 
1351  /*
1352  * Update the temperature dependence of the pitzer coefficients
1353  * and their derivatives
1354  */
1356 
1357  /*
1358  * Calculate the IMS cutoff factors
1359  */
1361 
1362  /*
1363  * Now do the main calculation.
1364  */
1366 
1367  double xmolSolvent = moleFraction(m_indexSolvent);
1368  double xx = std::max(m_xmolSolventMIN, xmolSolvent);
1369  double lnActCoeffMolal0 = - log(xx) + (xx - 1.0)/xx;
1370  double lnxs = log(xx);
1371 
1372  for (size_t k = 1; k < m_kk; k++) {
1373  CROP_speciesCropped_[k] = 0;
1375  if (m_lnActCoeffMolal_Unscaled[k] > (CROP_ln_gamma_k_max- 2.5 *lnxs)) {
1376  CROP_speciesCropped_[k] = 2;
1377  m_lnActCoeffMolal_Unscaled[k] = CROP_ln_gamma_k_max - 2.5 * lnxs;
1378  }
1379  if (m_lnActCoeffMolal_Unscaled[k] < (CROP_ln_gamma_k_min - 2.5 *lnxs)) {
1380  // -1.0 and -1.5 caused multiple solutions
1381  CROP_speciesCropped_[k] = 2;
1382  m_lnActCoeffMolal_Unscaled[k] = CROP_ln_gamma_k_min - 2.5 * lnxs;
1383  }
1384  }
1385  CROP_speciesCropped_[0] = 0;
1386  m_lnActCoeffMolal_Unscaled[0] += (IMS_lnActCoeffMolal_[0] - lnActCoeffMolal0);
1387  if (m_lnActCoeffMolal_Unscaled[0] < CROP_ln_gamma_o_min) {
1388  CROP_speciesCropped_[0] = 2;
1389  m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_min;
1390  }
1391  if (m_lnActCoeffMolal_Unscaled[0] > CROP_ln_gamma_o_max) {
1392  CROP_speciesCropped_[0] = 2;
1393  // -0.5 caused multiple solutions
1394  m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_max;
1395  }
1396  if (m_lnActCoeffMolal_Unscaled[0] > CROP_ln_gamma_o_max - 0.5 * lnxs) {
1397  CROP_speciesCropped_[0] = 2;
1398  m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_max - 0.5 * lnxs;
1399  }
1400 
1401  /*
1402  * Now do the pH Scaling
1403  */
1405 }
1406 
1408 {
1409  doublereal Imax = 0.0, Itmp;
1410  doublereal Iac_max;
1411  m_molalitiesAreCropped = false;
1412 
1413  for (size_t k = 0; k < m_kk; k++) {
1415  Itmp = m_molalities[k] * charge(k) * charge(k);
1416  if (Itmp > Imax) {
1417  Imax = Itmp;
1418  }
1419  }
1420 
1421  int cropMethod = 1;
1422 
1423 
1424  if (cropMethod == 0) {
1425 
1426  /*
1427  * Quick return
1428  */
1429  if (Imax < m_maxIionicStrength) {
1430  return;
1431  }
1432 
1433  m_molalitiesAreCropped = true;
1434 
1435  for (size_t i = 1; i < (m_kk - 1); i++) {
1436  double charge_i = charge(i);
1437  double abs_charge_i = fabs(charge_i);
1438  if (charge_i == 0.0) {
1439  continue;
1440  }
1441  for (size_t j = (i+1); j < m_kk; j++) {
1442  double charge_j = charge(j);
1443  double abs_charge_j = fabs(charge_j);
1444  /*
1445  * Find the counterIJ for the symmetric binary interaction
1446  */
1447  // n = m_kk*i + j;
1448  // counterIJ = m_CounterIJ[n];
1449  /*
1450  * Only loop over oppositely charge species
1451  */
1452  if (charge_i * charge_j < 0) {
1453  Iac_max = m_maxIionicStrength;
1454 
1456  Imax = m_molalitiesCropped[i] * abs_charge_i * abs_charge_i;
1457  if (Imax > Iac_max) {
1458  m_molalitiesCropped[i] = Iac_max / (abs_charge_i * abs_charge_i);
1459  }
1460  Imax = m_molalitiesCropped[j] * fabs(abs_charge_j * abs_charge_i);
1461  if (Imax > Iac_max) {
1462  m_molalitiesCropped[j] = Iac_max / (abs_charge_j * abs_charge_i);
1463  }
1464  } else {
1465  Imax = m_molalitiesCropped[j] * abs_charge_j * abs_charge_j;
1466  if (Imax > Iac_max) {
1467  m_molalitiesCropped[j] = Iac_max / (abs_charge_j * abs_charge_j);
1468  }
1469  Imax = m_molalitiesCropped[i] * abs_charge_j * abs_charge_i;
1470  if (Imax > Iac_max) {
1471  m_molalitiesCropped[i] = Iac_max / (abs_charge_j * abs_charge_i);
1472  }
1473  }
1474  }
1475  }
1476  }
1477 
1478  /*
1479  * Do this loop 10 times until we have achieved charge neutrality
1480  * in the cropped molalities
1481  */
1482  for (int times = 0; times< 10; times++) {
1483  double anion_charge = 0.0;
1484  double cation_charge = 0.0;
1485  size_t anion_contrib_max_i = npos;
1486  double anion_contrib_max = -1.0;
1487  size_t cation_contrib_max_i = npos;
1488  double cation_contrib_max = -1.0;
1489  for (size_t i = 0; i < m_kk; i++) {
1490  double charge_i = charge(i);
1491  if (charge_i < 0.0) {
1492  double anion_contrib = - m_molalitiesCropped[i] * charge_i;
1493  anion_charge += anion_contrib ;
1494  if (anion_contrib > anion_contrib_max) {
1495  anion_contrib_max = anion_contrib;
1496  anion_contrib_max_i = i;
1497  }
1498  } else if (charge_i > 0.0) {
1499  double cation_contrib = m_molalitiesCropped[i] * charge_i;
1500  cation_charge += cation_contrib ;
1501  if (cation_contrib > cation_contrib_max) {
1502  cation_contrib_max = cation_contrib;
1503  cation_contrib_max_i = i;
1504  }
1505  }
1506  }
1507  double total_charge = cation_charge - anion_charge;
1508  if (total_charge > 1.0E-8) {
1509  double desiredCrop = total_charge/charge(cation_contrib_max_i);
1510  double maxCrop = 0.66 * m_molalitiesCropped[cation_contrib_max_i];
1511  if (desiredCrop < maxCrop) {
1512  m_molalitiesCropped[cation_contrib_max_i] -= desiredCrop;
1513  break;
1514  } else {
1515  m_molalitiesCropped[cation_contrib_max_i] -= maxCrop;
1516  }
1517  } else if (total_charge < -1.0E-8) {
1518  double desiredCrop = total_charge/charge(anion_contrib_max_i);
1519  double maxCrop = 0.66 * m_molalitiesCropped[anion_contrib_max_i];
1520  if (desiredCrop < maxCrop) {
1521  m_molalitiesCropped[anion_contrib_max_i] -= desiredCrop;
1522  break;
1523  } else {
1524  m_molalitiesCropped[anion_contrib_max_i] -= maxCrop;
1525  }
1526  } else {
1527  break;
1528  }
1529  }
1530  }
1531 
1532  if (cropMethod == 1) {
1533  double* molF = DATA_PTR(m_gamma_tmp);
1534  getMoleFractions(molF);
1535  double xmolSolvent = molF[m_indexSolvent];
1536  if (xmolSolvent >= MC_X_o_cutoff_) {
1537  return;
1538  }
1539 
1540  m_molalitiesAreCropped = true;
1541 
1542  double poly = MC_apCut_ + MC_bpCut_ * xmolSolvent + MC_dpCut_* xmolSolvent * xmolSolvent;
1543  double p = xmolSolvent + MC_epCut_ + exp(- xmolSolvent/ MC_cpCut_) * poly;
1544  double denomInv = 1.0/ (m_Mnaught * p);
1545 
1546  for (size_t k = 0; k < m_kk; k++) {
1547  m_molalitiesCropped[k] = molF[k] * denomInv ;
1548  }
1549 
1550  // Do a further check to see if the Ionic strength is below a max value
1551  // Reduce the molalities to enforce this. Note, this algorithm preserves
1552  // the charge neutrality of the solution after cropping.
1553  Itmp = 0.0;
1554  for (size_t k = 0; k < m_kk; k++) {
1555  Itmp += m_molalitiesCropped[k] * charge(k) * charge(k);
1556  }
1557  if (Itmp > m_maxIionicStrength) {
1558  double ratio = Itmp / m_maxIionicStrength;
1559  for (size_t k = 0; k < m_kk; k++) {
1560  if (charge(k) != 0.0) {
1561  m_molalitiesCropped[k] *= ratio;
1562  }
1563  }
1564  }
1565  }
1566 
1567 }
1568 
1570 {
1571  size_t n, nc, i, j;
1572  m_CounterIJ.resize(m_kk * m_kk);
1573  int counter = 0;
1574  for (i = 0; i < m_kk; i++) {
1575  n = i;
1576  nc = m_kk * i;
1577  m_CounterIJ[n] = 0;
1578  m_CounterIJ[nc] = 0;
1579  }
1580  for (i = 1; i < (m_kk - 1); i++) {
1581  n = m_kk * i + i;
1582  m_CounterIJ[n] = 0;
1583  for (j = (i+1); j < m_kk; j++) {
1584  n = m_kk * j + i;
1585  nc = m_kk * i + j;
1586  counter++;
1587  m_CounterIJ[n] = counter;
1588  m_CounterIJ[nc] = counter;
1589  }
1590  }
1591 }
1592 
1593 void HMWSoln::s_updatePitzer_CoeffWRTemp(int doDerivs) const
1594 {
1595 
1596  size_t i, j, n, counterIJ;
1597  const double* beta0MX_coeff;
1598  const double* beta1MX_coeff;
1599  const double* beta2MX_coeff;
1600  const double* CphiMX_coeff;
1601  const double* Theta_coeff;
1602  double T = temperature();
1603  double Tr = m_TempPitzerRef;
1604  double tinv = 0.0, tln = 0.0, tlin = 0.0, tquad = 0.0;
1605  if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
1606  tlin = T - Tr;
1607  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
1608  tlin = T - Tr;
1609  tquad = T * T - Tr * Tr;
1610  tln = log(T/ Tr);
1611  tinv = 1.0/T - 1.0/Tr;
1612  }
1613 
1614  for (i = 1; i < (m_kk - 1); i++) {
1615  for (j = (i+1); j < m_kk; j++) {
1616 
1617  /*
1618  * Find the counterIJ for the symmetric binary interaction
1619  */
1620  n = m_kk*i + j;
1621  counterIJ = m_CounterIJ[n];
1622 
1623  beta0MX_coeff = m_Beta0MX_ij_coeff.ptrColumn(counterIJ);
1624  beta1MX_coeff = m_Beta1MX_ij_coeff.ptrColumn(counterIJ);
1625  beta2MX_coeff = m_Beta2MX_ij_coeff.ptrColumn(counterIJ);
1626  CphiMX_coeff = m_CphiMX_ij_coeff.ptrColumn(counterIJ);
1627  Theta_coeff = m_Theta_ij_coeff.ptrColumn(counterIJ);
1628 
1629  switch (m_formPitzerTemp) {
1630  case PITZER_TEMP_CONSTANT:
1631  break;
1632  case PITZER_TEMP_LINEAR:
1633 
1634  m_Beta0MX_ij[counterIJ] = beta0MX_coeff[0]
1635  + beta0MX_coeff[1]*tlin;
1636  m_Beta0MX_ij_L[counterIJ] = beta0MX_coeff[1];
1637  m_Beta0MX_ij_LL[counterIJ] = 0.0;
1638 
1639  m_Beta1MX_ij[counterIJ] = beta1MX_coeff[0]
1640  + beta1MX_coeff[1]*tlin;
1641  m_Beta1MX_ij_L[counterIJ] = beta1MX_coeff[1];
1642  m_Beta1MX_ij_LL[counterIJ] = 0.0;
1643 
1644  m_Beta2MX_ij[counterIJ] = beta2MX_coeff[0]
1645  + beta2MX_coeff[1]*tlin;
1646  m_Beta2MX_ij_L[counterIJ] = beta2MX_coeff[1];
1647  m_Beta2MX_ij_LL[counterIJ] = 0.0;
1648 
1649  m_CphiMX_ij[counterIJ] = CphiMX_coeff[0]
1650  + CphiMX_coeff[1]*tlin;
1651  m_CphiMX_ij_L[counterIJ] = CphiMX_coeff[1];
1652  m_CphiMX_ij_LL[counterIJ] = 0.0;
1653 
1654  m_Theta_ij[counterIJ] = Theta_coeff[0] + Theta_coeff[1]*tlin;
1655  m_Theta_ij_L[counterIJ] = Theta_coeff[1];
1656  m_Theta_ij_LL[counterIJ] = 0.0;
1657 
1658  break;
1659 
1660  case PITZER_TEMP_COMPLEX1:
1661  m_Beta0MX_ij[counterIJ] = beta0MX_coeff[0]
1662  + beta0MX_coeff[1]*tlin
1663  + beta0MX_coeff[2]*tquad
1664  + beta0MX_coeff[3]*tinv
1665  + beta0MX_coeff[4]*tln;
1666 
1667  m_Beta1MX_ij[counterIJ] = beta1MX_coeff[0]
1668  + beta1MX_coeff[1]*tlin
1669  + beta1MX_coeff[2]*tquad
1670  + beta1MX_coeff[3]*tinv
1671  + beta1MX_coeff[4]*tln;
1672 
1673  m_Beta2MX_ij[counterIJ] = beta2MX_coeff[0]
1674  + beta2MX_coeff[1]*tlin
1675  + beta2MX_coeff[2]*tquad
1676  + beta2MX_coeff[3]*tinv
1677  + beta2MX_coeff[4]*tln;
1678 
1679  m_CphiMX_ij[counterIJ] = CphiMX_coeff[0]
1680  + CphiMX_coeff[1]*tlin
1681  + CphiMX_coeff[2]*tquad
1682  + CphiMX_coeff[3]*tinv
1683  + CphiMX_coeff[4]*tln;
1684 
1685  m_Theta_ij[counterIJ] = Theta_coeff[0]
1686  + Theta_coeff[1]*tlin
1687  + Theta_coeff[2]*tquad
1688  + Theta_coeff[3]*tinv
1689  + Theta_coeff[4]*tln;
1690 
1691  m_Beta0MX_ij_L[counterIJ] = beta0MX_coeff[1]
1692  + beta0MX_coeff[2]*2.0*T
1693  - beta0MX_coeff[3]/(T*T)
1694  + beta0MX_coeff[4]/T;
1695 
1696  m_Beta1MX_ij_L[counterIJ] = beta1MX_coeff[1]
1697  + beta1MX_coeff[2]*2.0*T
1698  - beta1MX_coeff[3]/(T*T)
1699  + beta1MX_coeff[4]/T;
1700 
1701  m_Beta2MX_ij_L[counterIJ] = beta2MX_coeff[1]
1702  + beta2MX_coeff[2]*2.0*T
1703  - beta2MX_coeff[3]/(T*T)
1704  + beta2MX_coeff[4]/T;
1705 
1706  m_CphiMX_ij_L[counterIJ] = CphiMX_coeff[1]
1707  + CphiMX_coeff[2]*2.0*T
1708  - CphiMX_coeff[3]/(T*T)
1709  + CphiMX_coeff[4]/T;
1710 
1711  m_Theta_ij_L[counterIJ] = Theta_coeff[1]
1712  + Theta_coeff[2]*2.0*T
1713  - Theta_coeff[3]/(T*T)
1714  + Theta_coeff[4]/T;
1715 
1716  doDerivs = 2;
1717  if (doDerivs > 1) {
1718  m_Beta0MX_ij_LL[counterIJ] =
1719  + beta0MX_coeff[2]*2.0
1720  + 2.0*beta0MX_coeff[3]/(T*T*T)
1721  - beta0MX_coeff[4]/(T*T);
1722 
1723  m_Beta1MX_ij_LL[counterIJ] =
1724  + beta1MX_coeff[2]*2.0
1725  + 2.0*beta1MX_coeff[3]/(T*T*T)
1726  - beta1MX_coeff[4]/(T*T);
1727 
1728  m_Beta2MX_ij_LL[counterIJ] =
1729  + beta2MX_coeff[2]*2.0
1730  + 2.0*beta2MX_coeff[3]/(T*T*T)
1731  - beta2MX_coeff[4]/(T*T);
1732 
1733  m_CphiMX_ij_LL[counterIJ] =
1734  + CphiMX_coeff[2]*2.0
1735  + 2.0*CphiMX_coeff[3]/(T*T*T)
1736  - CphiMX_coeff[4]/(T*T);
1737 
1738  m_Theta_ij_LL[counterIJ] =
1739  + Theta_coeff[2]*2.0
1740  + 2.0*Theta_coeff[3]/(T*T*T)
1741  - Theta_coeff[4]/(T*T);
1742  }
1743 
1744 #ifdef DEBUG_HKM
1745  /*
1746  * Turn terms off for debugging
1747  */
1748  //m_Beta0MX_ij_L[counterIJ] = 0;
1749  //m_Beta0MX_ij_LL[counterIJ] = 0;
1750  //m_Beta1MX_ij_L[counterIJ] = 0;
1751  //m_Beta1MX_ij_LL[counterIJ] = 0;
1752  //m_CphiMX_ij_L[counterIJ] = 0;
1753  //m_CphiMX_ij_LL[counterIJ] = 0;
1754 #endif
1755  break;
1756  }
1757  }
1758  }
1759 
1760  // Lambda interactions and Mu_nnn
1761  // i must be neutral for this term to be nonzero. We take advantage of this
1762  // here to lower the operation count.
1763  for (i = 1; i < m_kk; i++) {
1764  if (charge(i) == 0.0) {
1765  for (j = 1; j < m_kk; j++) {
1766  n = i * m_kk + j;
1767  const double* Lambda_coeff = m_Lambda_nj_coeff.ptrColumn(n);
1768  switch (m_formPitzerTemp) {
1769  case PITZER_TEMP_CONSTANT:
1770  m_Lambda_nj(i,j) = Lambda_coeff[0];
1771  break;
1772  case PITZER_TEMP_LINEAR:
1773  m_Lambda_nj(i,j) = Lambda_coeff[0] + Lambda_coeff[1]*tlin;
1774  m_Lambda_nj_L(i,j) = Lambda_coeff[1];
1775  m_Lambda_nj_LL(i,j) = 0.0;
1776  break;
1777  case PITZER_TEMP_COMPLEX1:
1778  m_Lambda_nj(i,j) = Lambda_coeff[0]
1779  + Lambda_coeff[1]*tlin
1780  + Lambda_coeff[2]*tquad
1781  + Lambda_coeff[3]*tinv
1782  + Lambda_coeff[4]*tln;
1783 
1784  m_Lambda_nj_L(i,j) = Lambda_coeff[1]
1785  + Lambda_coeff[2]*2.0*T
1786  - Lambda_coeff[3]/(T*T)
1787  + Lambda_coeff[4]/T;
1788 
1789  m_Lambda_nj_LL(i,j) =
1790  Lambda_coeff[2]*2.0
1791  + 2.0*Lambda_coeff[3]/(T*T*T)
1792  - Lambda_coeff[4]/(T*T);
1793  }
1794 
1795  if (j == i) {
1796  const double* Mu_coeff = m_Mu_nnn_coeff.ptrColumn(i);
1797  switch (m_formPitzerTemp) {
1798  case PITZER_TEMP_CONSTANT:
1799  m_Mu_nnn[i] = Mu_coeff[0];
1800  break;
1801  case PITZER_TEMP_LINEAR:
1802  m_Mu_nnn[i] = Mu_coeff[0] + Mu_coeff[1]*tlin;
1803  m_Mu_nnn_L[i] = Mu_coeff[1];
1804  m_Mu_nnn_LL[i] = 0.0;
1805  break;
1806  case PITZER_TEMP_COMPLEX1:
1807  m_Mu_nnn[i] = Mu_coeff[0]
1808  + Mu_coeff[1]*tlin
1809  + Mu_coeff[2]*tquad
1810  + Mu_coeff[3]*tinv
1811  + Mu_coeff[4]*tln;
1812 
1813  m_Mu_nnn_L[i] = Mu_coeff[1]
1814  + Mu_coeff[2]*2.0*T
1815  - Mu_coeff[3]/(T*T)
1816  + Mu_coeff[4]/T;
1817 
1818  m_Mu_nnn_LL[i] =
1819  Mu_coeff[2]*2.0
1820  + 2.0*Mu_coeff[3]/(T*T*T)
1821  - Mu_coeff[4]/(T*T);
1822  }
1823  }
1824  }
1825  }
1826  }
1827 
1828 
1829  for (i = 1; i < m_kk; i++) {
1830  for (j = 1; j < m_kk; j++) {
1831  for (size_t k = 1; k < m_kk; k++) {
1832  n = i * m_kk *m_kk + j * m_kk + k ;
1833  const double* Psi_coeff = m_Psi_ijk_coeff.ptrColumn(n);
1834  switch (m_formPitzerTemp) {
1835  case PITZER_TEMP_CONSTANT:
1836  m_Psi_ijk[n] = Psi_coeff[0];
1837  break;
1838  case PITZER_TEMP_LINEAR:
1839  m_Psi_ijk[n] = Psi_coeff[0] + Psi_coeff[1]*tlin;
1840  m_Psi_ijk_L[n] = Psi_coeff[1];
1841  m_Psi_ijk_LL[n] = 0.0;
1842  break;
1843  case PITZER_TEMP_COMPLEX1:
1844  m_Psi_ijk[n] = Psi_coeff[0]
1845  + Psi_coeff[1]*tlin
1846  + Psi_coeff[2]*tquad
1847  + Psi_coeff[3]*tinv
1848  + Psi_coeff[4]*tln;
1849 
1850  m_Psi_ijk_L[n] = Psi_coeff[1]
1851  + Psi_coeff[2]*2.0*T
1852  - Psi_coeff[3]/(T*T)
1853  + Psi_coeff[4]/T;
1854 
1855  m_Psi_ijk_LL[n] =
1856  Psi_coeff[2]*2.0
1857  + 2.0*Psi_coeff[3]/(T*T*T)
1858  - Psi_coeff[4]/(T*T);
1859  }
1860  }
1861  }
1862  }
1863 
1864 }
1865 
1866 void HMWSoln::
1868 {
1869  /*
1870  * HKM -> Assumption is made that the solvent is
1871  * species 0.
1872  */
1873  if (m_indexSolvent != 0) {
1874  printf("Wrong index solvent value!\n");
1875  exit(EXIT_FAILURE);
1876  }
1877 
1878 #ifdef DEBUG_MODE
1879  int printE = 0;
1880  if (temperature() == 323.15) {
1881  printE = 0;
1882  }
1883 #endif
1884  std::string sni, snj, snk;
1885 
1886  /*
1887  * Use the CROPPED molality of the species in solution.
1888  */
1889  const double* molality = DATA_PTR(m_molalitiesCropped);
1890 
1891  /*
1892  * These are data inputs about the Pitzer correlation. They come
1893  * from the input file for the Pitzer model.
1894  */
1895  const double* beta0MX = DATA_PTR(m_Beta0MX_ij);
1896  const double* beta1MX = DATA_PTR(m_Beta1MX_ij);
1897  const double* beta2MX = DATA_PTR(m_Beta2MX_ij);
1898  const double* CphiMX = DATA_PTR(m_CphiMX_ij);
1899  const double* thetaij = DATA_PTR(m_Theta_ij);
1900  const double* alpha1MX = DATA_PTR(m_Alpha1MX_ij);
1901  const double* alpha2MX = DATA_PTR(m_Alpha2MX_ij);
1902 
1903  const double* psi_ijk = DATA_PTR(m_Psi_ijk);
1904  //n = k + j * m_kk + i * m_kk * m_kk;
1905 
1906 
1907  double* gamma_Unscaled = DATA_PTR(m_gamma_tmp);
1908  /*
1909  * Local variables defined by Coltrin
1910  */
1911  double etheta[5][5], etheta_prime[5][5], sqrtIs;
1912  /*
1913  * Molality based ionic strength of the solution
1914  */
1915  double Is = 0.0;
1916  /*
1917  * Molarcharge of the solution: In Pitzer's notation,
1918  * this is his variable called "Z".
1919  */
1920  double molarcharge = 0.0;
1921  /*
1922  * molalitysum is the sum of the molalities over all solutes,
1923  * even those with zero charge.
1924  */
1925  double molalitysumUncropped = 0.0;
1926 
1927  double* gfunc = DATA_PTR(m_gfunc_IJ);
1928  double* g2func = DATA_PTR(m_g2func_IJ);
1929  double* hfunc = DATA_PTR(m_hfunc_IJ);
1930  double* h2func = DATA_PTR(m_h2func_IJ);
1931  double* BMX = DATA_PTR(m_BMX_IJ);
1932  double* BprimeMX = DATA_PTR(m_BprimeMX_IJ);
1933  double* BphiMX = DATA_PTR(m_BphiMX_IJ);
1934  double* Phi = DATA_PTR(m_Phi_IJ);
1935  double* Phiprime = DATA_PTR(m_Phiprime_IJ);
1936  double* Phiphi = DATA_PTR(m_PhiPhi_IJ);
1937  double* CMX = DATA_PTR(m_CMX_IJ);
1938 
1939 
1940  double x1, x2;
1941  double Aphi, F, zsqF;
1942  double sum1, sum2, sum3, sum4, sum5, term1;
1943  double sum_m_phi_minus_1, osmotic_coef, lnwateract;
1944 
1945  int z1, z2;
1946  size_t n, i, j, m, counterIJ, counterIJ2;
1947 
1948 #ifdef DEBUG_MODE
1949  if (m_debugCalc) {
1950  printf("\n Debugging information from hmw_act \n");
1951  }
1952 #endif
1953  /*
1954  * Make sure the counter variables are setup
1955  */
1956  counterIJ_setup();
1957 
1958  /*
1959  * ---------- Calculate common sums over solutes ---------------------
1960  */
1961  for (n = 1; n < m_kk; n++) {
1962  // ionic strength
1963  Is += charge(n) * charge(n) * molality[n];
1964  // total molar charge
1965  molarcharge += fabs(charge(n)) * molality[n];
1966  molalitysumUncropped += m_molalities[n];
1967  }
1968  Is *= 0.5;
1969 
1970  /*
1971  * Store the ionic molality in the object for reference.
1972  */
1973  m_IionicMolality = Is;
1974  sqrtIs = sqrt(Is);
1975 #ifdef DEBUG_MODE
1976  if (m_debugCalc) {
1977  printf(" Step 1: \n");
1978  printf(" ionic strenth = %14.7le \n total molar "
1979  "charge = %14.7le \n", Is, molarcharge);
1980  }
1981 #endif
1982 
1983  /*
1984  * The following call to calc_lambdas() calculates all 16 elements
1985  * of the elambda and elambda1 arrays, given the value of the
1986  * ionic strength (Is)
1987  */
1988  calc_lambdas(Is);
1989 
1990  /*
1991  * ----- Step 2: Find the coefficients E-theta and -------------------
1992  * E-thetaprime for all combinations of positive
1993  * unlike charges up to 4
1994  */
1995 #ifdef DEBUG_MODE
1996  if (m_debugCalc) {
1997  printf(" Step 2: \n");
1998  }
1999 #endif
2000  for (z1 = 1; z1 <=4; z1++) {
2001  for (z2 =1; z2 <=4; z2++) {
2002  calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
2003 #ifdef DEBUG_MODE
2004  if (m_debugCalc) {
2005  printf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n",
2006  z1, z2, etheta[z1][z2], etheta_prime[z1][z2]);
2007  }
2008 #endif
2009  }
2010  }
2011 
2012 #ifdef DEBUG_MODE
2013  if (m_debugCalc) {
2014  printf(" Step 3: \n");
2015  printf(" Species Species g(x) "
2016  " hfunc(x) \n");
2017  }
2018 #endif
2019 
2020  /*
2021  * calculate g(x) and hfunc(x) for each cation-anion pair MX
2022  * In the original literature, hfunc, was called gprime. However,
2023  * it's not the derivative of g(x), so I renamed it.
2024  */
2025  for (i = 1; i < (m_kk - 1); i++) {
2026  for (j = (i+1); j < m_kk; j++) {
2027  /*
2028  * Find the counterIJ for the symmetric binary interaction
2029  */
2030  n = m_kk*i + j;
2031  counterIJ = m_CounterIJ[n];
2032 
2033  /*
2034  * Only loop over oppositely charge species
2035  */
2036  if (charge(i)*charge(j) < 0) {
2037  /*
2038  * x is a reduced function variable
2039  */
2040  x1 = sqrtIs * alpha1MX[counterIJ];
2041  if (x1 > 1.0E-100) {
2042  gfunc[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 * x1);
2043  hfunc[counterIJ] = -2.0 *
2044  (1.0-(1.0 + x1 + 0.5 * x1 * x1) * exp(-x1)) / (x1 * x1);
2045  } else {
2046  gfunc[counterIJ] = 0.0;
2047  hfunc[counterIJ] = 0.0;
2048  }
2049 
2050  if (beta2MX[counterIJ] != 0.0) {
2051  x2 = sqrtIs * alpha2MX[counterIJ];
2052  if (x2 > 1.0E-100) {
2053  g2func[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
2054  h2func[counterIJ] = -2.0 *
2055  (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
2056  } else {
2057  g2func[counterIJ] = 0.0;
2058  h2func[counterIJ] = 0.0;
2059  }
2060  }
2061  } else {
2062  gfunc[counterIJ] = 0.0;
2063  hfunc[counterIJ] = 0.0;
2064  }
2065 #ifdef DEBUG_MODE
2066  if (m_debugCalc) {
2067  sni = speciesName(i);
2068  snj = speciesName(j);
2069  printf(" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(),
2070  gfunc[counterIJ], hfunc[counterIJ]);
2071  }
2072 #endif
2073  }
2074  }
2075 
2076  /*
2077  * --------- SUBSECTION TO CALCULATE BMX, BprimeMX, BphiMX ----------
2078  * --------- Agrees with Pitzer, Eq. (49), (51), (55)
2079  */
2080 #ifdef DEBUG_MODE
2081  if (m_debugCalc) {
2082  printf(" Step 4: \n");
2083  printf(" Species Species BMX "
2084  "BprimeMX BphiMX \n");
2085  }
2086 #endif
2087 
2088  for (i = 1; i < m_kk - 1; i++) {
2089  for (j = i+1; j < m_kk; j++) {
2090  /*
2091  * Find the counterIJ for the symmetric binary interaction
2092  */
2093  n = m_kk*i + j;
2094  counterIJ = m_CounterIJ[n];
2095 
2096 #ifdef DEBUG_MODE
2097  if (printE) {
2098  if (counterIJ == 2) {
2099  printf("%s %s\n", speciesName(i).c_str(),
2100  speciesName(j).c_str());
2101  printf("beta0MX[%d] = %g\n", (int) counterIJ, beta0MX[counterIJ]);
2102  printf("beta1MX[%d] = %g\n", (int) counterIJ, beta1MX[counterIJ]);
2103  printf("beta2MX[%d] = %g\n", (int) counterIJ, beta2MX[counterIJ]);
2104  }
2105  }
2106 #endif
2107  /*
2108  * both species have a non-zero charge, and one is positive
2109  * and the other is negative
2110  */
2111  if (charge(i)*charge(j) < 0.0) {
2112  BMX[counterIJ] = beta0MX[counterIJ]
2113  + beta1MX[counterIJ] * gfunc[counterIJ]
2114  + beta2MX[counterIJ] * g2func[counterIJ];
2115 #ifdef DEBUG_MODE
2116  if (m_debugCalc) {
2117  printf("%d %g: %g %g %g %g\n",
2118  (int) counterIJ, BMX[counterIJ], beta0MX[counterIJ],
2119  beta1MX[counterIJ], beta2MX[counterIJ], gfunc[counterIJ]);
2120  }
2121 #endif
2122  if (Is > 1.0E-150) {
2123  BprimeMX[counterIJ] = (beta1MX[counterIJ] * hfunc[counterIJ]/Is +
2124  beta2MX[counterIJ] * h2func[counterIJ]/Is);
2125  } else {
2126  BprimeMX[counterIJ] = 0.0;
2127  }
2128  BphiMX[counterIJ] = BMX[counterIJ] + Is*BprimeMX[counterIJ];
2129  } else {
2130  BMX[counterIJ] = 0.0;
2131  BprimeMX[counterIJ] = 0.0;
2132  BphiMX[counterIJ] = 0.0;
2133  }
2134 #ifdef DEBUG_MODE
2135  if (m_debugCalc) {
2136  sni = speciesName(i);
2137  snj = speciesName(j);
2138  printf(" %-16s %-16s %11.7f %11.7f %11.7f \n",
2139  sni.c_str(), snj.c_str(),
2140  BMX[counterIJ], BprimeMX[counterIJ], BphiMX[counterIJ]);
2141  }
2142 #endif
2143  }
2144  }
2145 
2146  /*
2147  * --------- SUBSECTION TO CALCULATE CMX ----------
2148  * --------- Agrees with Pitzer, Eq. (53).
2149  */
2150 #ifdef DEBUG_MODE
2151  if (m_debugCalc) {
2152  printf(" Step 5: \n");
2153  printf(" Species Species CMX \n");
2154  }
2155 #endif
2156  for (i = 1; i < m_kk-1; i++) {
2157  for (j = i+1; j < m_kk; j++) {
2158  /*
2159  * Find the counterIJ for the symmetric binary interaction
2160  */
2161  n = m_kk*i + j;
2162  counterIJ = m_CounterIJ[n];
2163  /*
2164  * both species have a non-zero charge, and one is positive
2165  * and the other is negative
2166  */
2167  if (charge(i)*charge(j) < 0.0) {
2168  CMX[counterIJ] = CphiMX[counterIJ]/
2169  (2.0* sqrt(fabs(charge(i)*charge(j))));
2170  } else {
2171  CMX[counterIJ] = 0.0;
2172  }
2173 #ifdef DEBUG_MODE
2174  if (printE) {
2175  if (counterIJ == 2) {
2176  printf("%s %s\n", speciesName(i).c_str(),
2177  speciesName(j).c_str());
2178  printf("CphiMX[%d] = %g\n", (int) counterIJ, CphiMX[counterIJ]);
2179  }
2180  }
2181 #endif
2182 #ifdef DEBUG_MODE
2183  if (m_debugCalc) {
2184  sni = speciesName(i);
2185  snj = speciesName(j);
2186  printf(" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(),
2187  CMX[counterIJ]);
2188  }
2189 #endif
2190  }
2191  }
2192 
2193  /*
2194  * ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ----------
2195  * --------- Agrees with Pitzer, Eq. 72, 73, 74
2196  */
2197 #ifdef DEBUG_MODE
2198  if (m_debugCalc) {
2199  printf(" Step 6: \n");
2200  printf(" Species Species Phi_ij "
2201  " Phiprime_ij Phi^phi_ij \n");
2202  }
2203 #endif
2204  for (i = 1; i < m_kk-1; i++) {
2205  for (j = i+1; j < m_kk; j++) {
2206  /*
2207  * Find the counterIJ for the symmetric binary interaction
2208  */
2209  n = m_kk*i + j;
2210  counterIJ = m_CounterIJ[n];
2211  /*
2212  * both species have a non-zero charge, and one is positive
2213  * and the other is negative
2214  */
2215  if (charge(i)*charge(j) > 0) {
2216  z1 = (int) fabs(charge(i));
2217  z2 = (int) fabs(charge(j));
2218  Phi[counterIJ] = thetaij[counterIJ] + etheta[z1][z2];
2219  Phiprime[counterIJ] = etheta_prime[z1][z2];
2220  Phiphi[counterIJ] = Phi[counterIJ] + Is * Phiprime[counterIJ];
2221  } else {
2222  Phi[counterIJ] = 0.0;
2223  Phiprime[counterIJ] = 0.0;
2224  Phiphi[counterIJ] = 0.0;
2225  }
2226 #ifdef DEBUG_MODE
2227  if (m_debugCalc) {
2228  sni = speciesName(i);
2229  snj = speciesName(j);
2230  printf(" %-16s %-16s %10.6f %10.6f %10.6f \n",
2231  sni.c_str(), snj.c_str(),
2232  Phi[counterIJ], Phiprime[counterIJ], Phiphi[counterIJ]);
2233  }
2234 #endif
2235  }
2236  }
2237 
2238  /*
2239  * ------------- SUBSECTION FOR CALCULATION OF F ----------------------
2240  * ------------ Agrees with Pitzer Eqn. (65) --------------------------
2241  */
2242 #ifdef DEBUG_MODE
2243  if (m_debugCalc) {
2244  printf(" Step 7: \n");
2245  }
2246 #endif
2247  // A_Debye_Huckel = 0.5092; (units = sqrt(kg/gmol))
2248  // A_Debye_Huckel = 0.5107; <- This value is used to match GWB data
2249  // ( A * ln(10) = 1.17593)
2250  // Aphi = A_Debye_Huckel * 2.30258509 / 3.0;
2251  Aphi = m_A_Debye / 3.0;
2252  F = -Aphi * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
2253  + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
2254 #ifdef DEBUG_MODE
2255  if (printE) {
2256  printf("Aphi = %20.13g\n", Aphi);
2257  }
2258 #endif
2259 #ifdef DEBUG_MODE
2260  if (m_debugCalc) {
2261  printf(" initial value of F = %10.6f \n", F);
2262  }
2263 #endif
2264  for (i = 1; i < m_kk-1; i++) {
2265  for (j = i+1; j < m_kk; j++) {
2266  /*
2267  * Find the counterIJ for the symmetric binary interaction
2268  */
2269  n = m_kk*i + j;
2270  counterIJ = m_CounterIJ[n];
2271  /*
2272  * both species have a non-zero charge, and one is positive
2273  * and the other is negative
2274  */
2275  if (charge(i)*charge(j) < 0) {
2276  F = F + molality[i]*molality[j] * BprimeMX[counterIJ];
2277  }
2278  /*
2279  * Both species have a non-zero charge, and they
2280  * have the same sign
2281  */
2282  if (charge(i)*charge(j) > 0) {
2283  F = F + molality[i]*molality[j] * Phiprime[counterIJ];
2284  }
2285 #ifdef DEBUG_MODE
2286  if (m_debugCalc) {
2287  printf(" F = %10.6f \n", F);
2288  }
2289 #endif
2290  }
2291  }
2292 #ifdef DEBUG_MODE
2293  if (m_debugCalc) {
2294  printf(" Step 8: Summing in All Contributions to Activity Coefficients \n");
2295  }
2296 #endif
2297 
2298  for (i = 1; i < m_kk; i++) {
2299 
2300  /*
2301  * -------- SUBSECTION FOR CALCULATING THE ACTCOEFF FOR CATIONS -----
2302  * -------- -> equations agree with my notes, Eqn. (118).
2303  * -> Equations agree with Pitzer, eqn.(63)
2304  */
2305  if (charge(i) > 0.0) {
2306 
2307 #ifdef DEBUG_MODE
2308  if (m_debugCalc) {
2309  sni = speciesName(i);
2310  printf(" Contributions to ln(ActCoeff_%s):\n", sni.c_str());
2311  }
2312 #endif
2313  // species i is the cation (positive) to calc the actcoeff
2314  zsqF = charge(i)*charge(i)*F;
2315 #ifdef DEBUG_MODE
2316  if (m_debugCalc) {
2317  printf(" Unary term: z*z*F = %10.5f\n", zsqF);
2318  }
2319 #endif
2320  sum1 = 0.0;
2321  sum2 = 0.0;
2322  sum3 = 0.0;
2323  sum4 = 0.0;
2324  sum5 = 0.0;
2325  for (j = 1; j < m_kk; j++) {
2326  /*
2327  * Find the counterIJ for the symmetric binary interaction
2328  */
2329  n = m_kk*i + j;
2330  counterIJ = m_CounterIJ[n];
2331 
2332  if (charge(j) < 0.0) {
2333  // sum over all anions
2334  sum1 = sum1 + molality[j]*
2335  (2.0*BMX[counterIJ] + molarcharge*CMX[counterIJ]);
2336 #ifdef DEBUG_MODE
2337  if (m_debugCalc) {
2338  snj = speciesName(j) + ":";
2339  printf(" Bin term with %-13s 2 m_j BMX = %10.5f\n", snj.c_str(),
2340  molality[j]*2.0*BMX[counterIJ]);
2341  printf(" m_j Z CMX = %10.5f\n",
2342  molality[j]* molarcharge*CMX[counterIJ]);
2343  }
2344 #endif
2345  if (j < m_kk-1) {
2346  /*
2347  * This term is the ternary interaction involving the
2348  * non-duplicate sum over double anions, j, k, with
2349  * respect to the cation, i.
2350  */
2351  for (size_t k = j+1; k < m_kk; k++) {
2352  // an inner sum over all anions
2353  if (charge(k) < 0.0) {
2354  n = k + j * m_kk + i * m_kk * m_kk;
2355  sum3 = sum3 + molality[j]*molality[k]*psi_ijk[n];
2356 #ifdef DEBUG_MODE
2357  if (m_debugCalc) {
2358  if (psi_ijk[n] != 0.0) {
2359  snj = speciesName(j) + "," + speciesName(k) + ":";
2360  printf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(),
2361  molality[j]*molality[k]*psi_ijk[n]);
2362  }
2363  }
2364 #endif
2365  }
2366  }
2367  }
2368  }
2369 
2370 
2371  if (charge(j) > 0.0) {
2372  // sum over all cations
2373  if (j != i) {
2374  sum2 = sum2 + molality[j]*(2.0*Phi[counterIJ]);
2375 #ifdef DEBUG_MODE
2376  if (m_debugCalc) {
2377  if ((molality[j] * Phi[counterIJ])!= 0.0) {
2378  snj = speciesName(j) + ":";
2379  printf(" Phi term with %-12s 2 m_j Phi_cc = %10.5f\n", snj.c_str(),
2380  molality[j]*(2.0*Phi[counterIJ]));
2381  }
2382  }
2383 #endif
2384  }
2385  for (size_t k = 1; k < m_kk; k++) {
2386  if (charge(k) < 0.0) {
2387  // two inner sums over anions
2388 
2389  n = k + j * m_kk + i * m_kk * m_kk;
2390  sum2 = sum2 + molality[j]*molality[k]*psi_ijk[n];
2391 #ifdef DEBUG_MODE
2392  if (m_debugCalc) {
2393  if (psi_ijk[n] != 0.0) {
2394  snj = speciesName(j) + "," + speciesName(k) + ":";
2395  printf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(),
2396  molality[j]*molality[k]*psi_ijk[n]);
2397  }
2398  }
2399 #endif
2400  /*
2401  * Find the counterIJ for the j,k interaction
2402  */
2403  n = m_kk*j + k;
2404  counterIJ2 = m_CounterIJ[n];
2405  sum4 = sum4 + (fabs(charge(i))*
2406  molality[j]*molality[k]*CMX[counterIJ2]);
2407 #ifdef DEBUG_MODE
2408  if (m_debugCalc) {
2409  if ((molality[j]*molality[k]*CMX[counterIJ2]) != 0.0) {
2410  snj = speciesName(j) + "," + speciesName(k) + ":";
2411  printf(" Tern CMX term on %-16s abs(z_i) m_j m_k CMX = %10.5f\n", snj.c_str(),
2412  fabs(charge(i))* molality[j]*molality[k]*CMX[counterIJ2]);
2413  }
2414  }
2415 #endif
2416  }
2417  }
2418  }
2419 
2420  /*
2421  * Handle neutral j species
2422  */
2423  if (charge(j) == 0) {
2424  sum5 = sum5 + molality[j]*2.0*m_Lambda_nj(j,i);
2425 #ifdef DEBUG_MODE
2426  if (m_debugCalc) {
2427  if ((molality[j]*2.0*m_Lambda_nj(j,i)) != 0.0) {
2428  snj = speciesName(j) + ":";
2429  printf(" Lambda term with %-12s 2 m_j lam_ji = %10.5f\n", snj.c_str(),
2430  molality[j]*2.0*m_Lambda_nj(j,i));
2431  }
2432  }
2433 #endif
2434  /*
2435  * Zeta interaction term
2436  */
2437  for (size_t k = 1; k < m_kk; k++) {
2438  if (charge(k) < 0.0) {
2439  size_t izeta = j;
2440  size_t jzeta = i;
2441  n = izeta * m_kk * m_kk + jzeta * m_kk + k;
2442  double zeta = psi_ijk[n];
2443  if (zeta != 0.0) {
2444  sum5 = sum5 + molality[j]*molality[k]*zeta;
2445 #ifdef DEBUG_MODE
2446  if (m_debugCalc) {
2447  snj = speciesName(j) + "," + speciesName(k) + ":";
2448  printf(" Zeta term on %-16s m_n m_a zeta_nMa = %10.5f\n", snj.c_str(),
2449  molality[j]*molality[k]*psi_ijk[n]);
2450  }
2451 #endif
2452  }
2453  }
2454  }
2455  }
2456  }
2457  /*
2458  * Add all of the contributions up to yield the log of the
2459  * solute activity coefficients (molality scale)
2460  */
2461  m_lnActCoeffMolal_Unscaled[i] = zsqF + sum1 + sum2 + sum3 + sum4 + sum5;
2462  gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]);
2463 #ifdef DEBUG_MODE
2464  if (m_debugCalc) {
2465  sni = speciesName(i);
2466  printf(" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f \n",
2467  sni.c_str(), m_lnActCoeffMolal_Unscaled[i], gamma_Unscaled[i]);
2468  }
2469 #endif
2470  }
2471 
2472  /*
2473  * -------- SUBSECTION FOR CALCULATING THE ACTCOEFF FOR ANIONS ------
2474  * -------- -> equations agree with my notes, Eqn. (119).
2475  * -> Equations agree with Pitzer, eqn.(64)
2476  */
2477  if (charge(i) < 0) {
2478 
2479 #ifdef DEBUG_MODE
2480  if (m_debugCalc) {
2481  sni = speciesName(i);
2482  printf(" Contributions to ln(ActCoeff_%s):\n", sni.c_str());
2483  }
2484 #endif
2485 
2486  // species i is an anion (negative)
2487  zsqF = charge(i)*charge(i)*F;
2488 #ifdef DEBUG_MODE
2489  if (m_debugCalc) {
2490  printf(" Unary term: z*z*F = %10.5f\n", zsqF);
2491  }
2492 #endif
2493  sum1 = 0.0;
2494  sum2 = 0.0;
2495  sum3 = 0.0;
2496  sum4 = 0.0;
2497  sum5 = 0.0;
2498  for (j = 1; j < m_kk; j++) {
2499  /*
2500  * Find the counterIJ for the symmetric binary interaction
2501  */
2502  n = m_kk*i + j;
2503  counterIJ = m_CounterIJ[n];
2504 
2505  /*
2506  * For Anions, do the cation interactions.
2507  */
2508  if (charge(j) > 0) {
2509  sum1 = sum1 + molality[j]*
2510  (2.0*BMX[counterIJ]+molarcharge*CMX[counterIJ]);
2511 #ifdef DEBUG_MODE
2512  if (m_debugCalc) {
2513  snj = speciesName(j) + ":";
2514  printf(" Bin term with %-13s 2 m_j BMX = %10.5f\n", snj.c_str(),
2515  molality[j]*2.0*BMX[counterIJ]);
2516  printf(" m_j Z CMX = %10.5f\n",
2517  molality[j]* molarcharge*CMX[counterIJ]);
2518  }
2519 #endif
2520  if (j < m_kk-1) {
2521  for (size_t k = j+1; k < m_kk; k++) {
2522  // an inner sum over all cations
2523  if (charge(k) > 0) {
2524  n = k + j * m_kk + i * m_kk * m_kk;
2525  sum3 = sum3 + molality[j]*molality[k]*psi_ijk[n];
2526 #ifdef DEBUG_MODE
2527  if (m_debugCalc) {
2528  if (psi_ijk[n] != 0.0) {
2529  snj = speciesName(j) + "," + speciesName(k) + ":";
2530  printf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(),
2531  molality[j]*molality[k]*psi_ijk[n]);
2532  }
2533  }
2534 #endif
2535  }
2536  }
2537  }
2538  }
2539 
2540  /*
2541  * For Anions, do the other anion interactions.
2542  */
2543  if (charge(j) < 0.0) {
2544  // sum over all anions
2545  if (j != i) {
2546  sum2 = sum2 + molality[j]*(2.0*Phi[counterIJ]);
2547 #ifdef DEBUG_MODE
2548  if (m_debugCalc) {
2549  if ((molality[j] * Phi[counterIJ])!= 0.0) {
2550  snj = speciesName(j) + ":";
2551  printf(" Phi term with %-12s 2 m_j Phi_aa = %10.5f\n", snj.c_str(),
2552  molality[j]*(2.0*Phi[counterIJ]));
2553  }
2554  }
2555 #endif
2556  }
2557  for (size_t k = 1; k < m_kk; k++) {
2558  if (charge(k) > 0.0) {
2559  // two inner sums over cations
2560  n = k + j * m_kk + i * m_kk * m_kk;
2561  sum2 = sum2 + molality[j]*molality[k]*psi_ijk[n];
2562 #ifdef DEBUG_MODE
2563  if (m_debugCalc) {
2564  if (psi_ijk[n] != 0.0) {
2565  snj = speciesName(j) + "," + speciesName(k) + ":";
2566  printf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(),
2567  molality[j]*molality[k]*psi_ijk[n]);
2568  }
2569  }
2570 #endif
2571  /*
2572  * Find the counterIJ for the symmetric binary interaction
2573  */
2574  n = m_kk*j + k;
2575  counterIJ2 = m_CounterIJ[n];
2576  sum4 = sum4 +
2577  (fabs(charge(i))*
2578  molality[j]*molality[k]*CMX[counterIJ2]);
2579 #ifdef DEBUG_MODE
2580  if (m_debugCalc) {
2581  if ((molality[j]*molality[k]*CMX[counterIJ2]) != 0.0) {
2582  snj = speciesName(j) + "," + speciesName(k) + ":";
2583  printf(" Tern CMX term on %-16s abs(z_i) m_j m_k CMX = %10.5f\n", snj.c_str(),
2584  fabs(charge(i))* molality[j]*molality[k]*CMX[counterIJ2]);
2585  }
2586  }
2587 #endif
2588  }
2589  }
2590  }
2591 
2592  /*
2593  * for Anions, do the neutral species interaction
2594  */
2595  if (charge(j) == 0.0) {
2596  sum5 = sum5 + molality[j]*2.0*m_Lambda_nj(j,i);
2597 #ifdef DEBUG_MODE
2598  if (m_debugCalc) {
2599  if ((molality[j]*2.0*m_Lambda_nj(j,i)) != 0.0) {
2600  snj = speciesName(j) + ":";
2601  printf(" Lambda term with %-12s 2 m_j lam_ji = %10.5f\n", snj.c_str(),
2602  molality[j]*2.0*m_Lambda_nj(j,i));
2603  }
2604  }
2605 #endif
2606  /*
2607  * Zeta interaction term
2608  */
2609  for (size_t k = 1; k < m_kk; k++) {
2610  if (charge(k) > 0.0) {
2611  size_t izeta = j;
2612  size_t jzeta = k;
2613  size_t kzeta = i;
2614  n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
2615  double zeta = psi_ijk[n];
2616  if (zeta != 0.0) {
2617  sum5 = sum5 + molality[j]*molality[k]*zeta;
2618 #ifdef DEBUG_MODE
2619  if (m_debugCalc) {
2620  snj = speciesName(j) + "," + speciesName(k) + ":";
2621  printf(" Zeta term on %-16s m_n m_c zeta_ncX = %10.5f\n", snj.c_str(),
2622  molality[j]*molality[k]*psi_ijk[n]);
2623  }
2624 #endif
2625  }
2626  }
2627  }
2628  }
2629  }
2630  m_lnActCoeffMolal_Unscaled[i] = zsqF + sum1 + sum2 + sum3 + sum4 + sum5;
2631  gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]);
2632 #ifdef DEBUG_MODE
2633  if (m_debugCalc) {
2634  sni = speciesName(i);
2635  printf(" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f\n",
2636  sni.c_str(), m_lnActCoeffMolal_Unscaled[i], gamma_Unscaled[i]);
2637  }
2638 #endif
2639  }
2640  /*
2641  * ------ SUBSECTION FOR CALCULATING NEUTRAL SOLUTE ACT COEFF -------
2642  * ------ -> equations agree with my notes,
2643  * -> Equations agree with Pitzer,
2644  */
2645  if (charge(i) == 0.0) {
2646 #ifdef DEBUG_MODE
2647  if (m_debugCalc) {
2648  sni = speciesName(i);
2649  printf(" Contributions to ln(ActCoeff_%s):\n", sni.c_str());
2650  }
2651 #endif
2652  sum1 = 0.0;
2653  sum3 = 0.0;
2654  for (j = 1; j < m_kk; j++) {
2655  sum1 = sum1 + molality[j]*2.0*m_Lambda_nj(i,j);
2656 #ifdef DEBUG_MODE
2657  if (m_debugCalc) {
2658  if (m_Lambda_nj(i,j) != 0.0) {
2659  snj = speciesName(j) + ":";
2660  printf(" Lambda_n term on %-16s 2 m_j lambda_n_j = %10.5f\n", snj.c_str(),
2661  molality[j]*2.0*m_Lambda_nj(i,j));
2662  }
2663  }
2664 #endif
2665 
2666  /*
2667  * Zeta term -> we piggyback on the psi term
2668  */
2669  if (charge(j) > 0.0) {
2670  for (size_t k = 1; k < m_kk; k++) {
2671  if (charge(k) < 0.0) {
2672  n = k + j * m_kk + i * m_kk * m_kk;
2673  sum3 = sum3 + molality[j]*molality[k]*psi_ijk[n];
2674 #ifdef DEBUG_MODE
2675  if (m_debugCalc) {
2676  if (psi_ijk[n] != 0.0) {
2677  snj = speciesName(j) + "," + speciesName(k) + ":";
2678  printf(" Zeta term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(),
2679  molality[j]*molality[k]*psi_ijk[n]);
2680  }
2681  }
2682 #endif
2683  }
2684  }
2685  }
2686  }
2687  sum2 = 3.0 * molality[i]* molality[i] * m_Mu_nnn[i];
2688 #ifdef DEBUG_MODE
2689  if (m_debugCalc) {
2690  if (m_Mu_nnn[i] != 0.0) {
2691  printf(" Mu_nnn term 3 m_n m_n Mu_n_n = %10.5f\n",
2692  3.0 * molality[i]* molality[i] * m_Mu_nnn[i]);
2693  }
2694  }
2695 #endif
2696 
2697  m_lnActCoeffMolal_Unscaled[i] = sum1 + sum2 + sum3;
2698  gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]);
2699 #ifdef DEBUG_MODE
2700  if (m_debugCalc) {
2701  sni = speciesName(i);
2702  printf(" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f\n",
2703  sni.c_str(), m_lnActCoeffMolal_Unscaled[i], gamma_Unscaled[i]);
2704  }
2705 #endif
2706  }
2707 
2708  }
2709 #ifdef DEBUG_MODE
2710  if (m_debugCalc) {
2711  printf(" Step 9: \n");
2712  }
2713 #endif
2714  /*
2715  * -------- SUBSECTION FOR CALCULATING THE OSMOTIC COEFF ---------
2716  * -------- -> equations agree with my notes, Eqn. (117).
2717  * -> Equations agree with Pitzer, eqn.(62)
2718  */
2719  sum1 = 0.0;
2720  sum2 = 0.0;
2721  sum3 = 0.0;
2722  sum4 = 0.0;
2723  sum5 = 0.0;
2724  double sum6 = 0.0;
2725  double sum7 = 0.0;
2726  /*
2727  * term1 is the DH term in the osmotic coefficient expression
2728  * b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer
2729  * implementations.
2730  * Is = Ionic strength on the molality scale (units of (gmol/kg))
2731  * Aphi = A_Debye / 3 (units of sqrt(kg/gmol))
2732  */
2733  term1 = -Aphi * pow(Is,1.5) / (1.0 + 1.2 * sqrt(Is));
2734 
2735  for (j = 1; j < m_kk; j++) {
2736  /*
2737  * Loop Over Cations
2738  */
2739  if (charge(j) > 0.0) {
2740  for (size_t k = 1; k < m_kk; k++) {
2741  if (charge(k) < 0.0) {
2742  /*
2743  * Find the counterIJ for the symmetric j,k binary interaction
2744  */
2745  n = m_kk*j + k;
2746  counterIJ = m_CounterIJ[n];
2747 
2748  sum1 = sum1 + molality[j]*molality[k]*
2749  (BphiMX[counterIJ] + molarcharge*CMX[counterIJ]);
2750  }
2751  }
2752 
2753  for (size_t k = j+1; k < m_kk; k++) {
2754  if (j == (m_kk-1)) {
2755  // we should never reach this step
2756  printf("logic error 1 in Step 9 of hmw_act");
2757  exit(EXIT_FAILURE);
2758  }
2759  if (charge(k) > 0.0) {
2760  /*
2761  * Find the counterIJ for the symmetric j,k binary interaction
2762  * between 2 cations.
2763  */
2764  n = m_kk*j + k;
2765  counterIJ = m_CounterIJ[n];
2766  sum2 = sum2 + molality[j]*molality[k]*Phiphi[counterIJ];
2767  for (m = 1; m < m_kk; m++) {
2768  if (charge(m) < 0.0) {
2769  // species m is an anion
2770  n = m + k * m_kk + j * m_kk * m_kk;
2771  sum2 = sum2 +
2772  molality[j]*molality[k]*molality[m]*psi_ijk[n];
2773  }
2774  }
2775  }
2776  }
2777  }
2778 
2779  /*
2780  * Loop Over Anions
2781  */
2782  if (charge(j) < 0) {
2783  for (size_t k = j+1; k < m_kk; k++) {
2784  if (j == m_kk-1) {
2785  // we should never reach this step
2786  printf("logic error 2 in Step 9 of hmw_act");
2787  exit(EXIT_FAILURE);
2788  }
2789  if (charge(k) < 0) {
2790  /*
2791  * Find the counterIJ for the symmetric j,k binary interaction
2792  * between two anions
2793  */
2794  n = m_kk*j + k;
2795  counterIJ = m_CounterIJ[n];
2796 
2797  sum3 = sum3 + molality[j]*molality[k]*Phiphi[counterIJ];
2798  for (m = 1; m < m_kk; m++) {
2799  if (charge(m) > 0.0) {
2800  n = m + k * m_kk + j * m_kk * m_kk;
2801  sum3 = sum3 +
2802  molality[j]*molality[k]*molality[m]*psi_ijk[n];
2803  }
2804  }
2805  }
2806  }
2807  }
2808 
2809  /*
2810  * Loop Over Neutral Species
2811  */
2812  if (charge(j) == 0) {
2813  for (size_t k = 1; k < m_kk; k++) {
2814  if (charge(k) < 0.0) {
2815  sum4 = sum4 + molality[j]*molality[k]*m_Lambda_nj(j,k);
2816  }
2817  if (charge(k) > 0.0) {
2818  sum5 = sum5 + molality[j]*molality[k]*m_Lambda_nj(j,k);
2819  }
2820  if (charge(k) == 0.0) {
2821  if (k > j) {
2822  sum6 = sum6 + molality[j]*molality[k]*m_Lambda_nj(j,k);
2823  } else if (k == j) {
2824  sum6 = sum6 + 0.5 * molality[j]*molality[k]*m_Lambda_nj(j,k);
2825  }
2826  }
2827  if (charge(k) < 0.0) {
2828  size_t izeta = j;
2829  for (m = 1; m < m_kk; m++) {
2830  if (charge(m) > 0.0) {
2831  size_t jzeta = m;
2832  n = k + jzeta * m_kk + izeta * m_kk * m_kk;
2833  double zeta = psi_ijk[n];
2834  if (zeta != 0.0) {
2835  sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta;
2836  }
2837  }
2838  }
2839  }
2840  }
2841  sum7 += molality[j]*molality[j]*molality[j]*m_Mu_nnn[j];
2842  }
2843  }
2844  sum_m_phi_minus_1 = 2.0 *
2845  (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
2846  /*
2847  * Calculate the osmotic coefficient from
2848  * osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i)
2849  */
2850  if (molalitysumUncropped > 1.0E-150) {
2851  osmotic_coef = 1.0 + (sum_m_phi_minus_1 / molalitysumUncropped);
2852  } else {
2853  osmotic_coef = 1.0;
2854  }
2855 #ifdef DEBUG_MODE
2856  if (printE) {
2857  printf("OsmCoef - 1 = %20.13g\n", osmotic_coef - 1.0);
2858  }
2859 #endif
2860 #ifdef DEBUG_MODE
2861  if (m_debugCalc) {
2862  printf(" term1=%10.6f sum1=%10.6f sum2=%10.6f "
2863  "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n",
2864  term1, sum1, sum2, sum3, sum4, sum5);
2865  printf(" sum_m_phi_minus_1=%10.6f osmotic_coef=%10.6f\n",
2866  sum_m_phi_minus_1, osmotic_coef);
2867  }
2868 
2869  if (m_debugCalc) {
2870  printf(" Step 10: \n");
2871  }
2872 #endif
2873  lnwateract = -(m_weightSolvent/1000.0) * molalitysumUncropped * osmotic_coef;
2874 
2875  /*
2876  * In Cantera, we define the activity coefficient of the solvent as
2877  *
2878  * act_0 = actcoeff_0 * Xmol_0
2879  *
2880  * We have just computed act_0. However, this routine returns
2881  * ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
2882  */
2883  double xmolSolvent = moleFraction(m_indexSolvent);
2884  double xx = std::max(m_xmolSolventMIN, xmolSolvent);
2885  m_lnActCoeffMolal_Unscaled[0] = lnwateract - log(xx);
2886 #ifdef DEBUG_MODE
2887  if (m_debugCalc) {
2888  double wateract = exp(lnwateract);
2889  printf(" Weight of Solvent = %16.7g\n", m_weightSolvent);
2890  printf(" molalitySumUncropped = %16.7g\n", molalitysumUncropped);
2891  printf(" ln_a_water=%10.6f a_water=%10.6f\n\n",
2892  lnwateract, wateract);
2893  }
2894 #endif
2895 }
2896 
2898 {
2899  /*
2900  * Zero the unscaled 2nd derivatives
2901  */
2902  m_dlnActCoeffMolaldT_Unscaled.assign(m_kk, 0.0);
2903  /*
2904  * Do the actual calculation of the unscaled temperature derivatives
2905  */
2907 
2908  //double xmolSolvent = moleFraction(m_indexSolvent);
2909  //double xx = MAX(m_xmolSolventMIN, xmolSolvent);
2910  // double lnxs = log(xx);
2911 
2912  for (size_t k = 1; k < m_kk; k++) {
2913  if (CROP_speciesCropped_[k] == 2) {
2915  }
2916  }
2917 
2918  if (CROP_speciesCropped_[0]) {
2920  }
2921 
2922  /*
2923  * Do the pH scaling to the derivatives
2924  */
2926 
2927 
2928 }
2929 
2931 {
2932  /*
2933  * It may be assumed that the Pitzer activity coefficient routine is
2934  * called immediately preceding the calling of this routine. Therefore,
2935  * some quantities do not need to be recalculated in this routine.
2936  */
2937 
2938  /*
2939  * HKM -> Assumption is made that the solvent is
2940  * species 0.
2941  */
2942 #ifdef DEBUG_MODE
2943  m_debugCalc = 0;
2944 #endif
2945  if (m_indexSolvent != 0) {
2946  printf("Wrong index solvent value!\n");
2947  exit(EXIT_FAILURE);
2948  }
2949 
2950  std::string sni, snj, snk;
2951 
2952  const double* molality = DATA_PTR(m_molalitiesCropped);
2953  const double* beta0MX_L = DATA_PTR(m_Beta0MX_ij_L);
2954  const double* beta1MX_L = DATA_PTR(m_Beta1MX_ij_L);
2955  const double* beta2MX_L = DATA_PTR(m_Beta2MX_ij_L);
2956  const double* CphiMX_L = DATA_PTR(m_CphiMX_ij_L);
2957  const double* thetaij_L = DATA_PTR(m_Theta_ij_L);
2958  const double* alpha1MX = DATA_PTR(m_Alpha1MX_ij);
2959  const double* alpha2MX = DATA_PTR(m_Alpha2MX_ij);
2960  const double* psi_ijk_L = DATA_PTR(m_Psi_ijk_L);
2961  double* d_gamma_dT_Unscaled = DATA_PTR(m_gamma_tmp);
2962  /*
2963  * Local variables defined by Coltrin
2964  */
2965  double etheta[5][5], etheta_prime[5][5], sqrtIs;
2966  /*
2967  * Molality based ionic strength of the solution
2968  */
2969  double Is = 0.0;
2970  /*
2971  * Molarcharge of the solution: In Pitzer's notation,
2972  * this is his variable called "Z".
2973  */
2974  double molarcharge = 0.0;
2975  /*
2976  * molalitysum is the sum of the molalities over all solutes,
2977  * even those with zero charge.
2978  */
2979  double molalitysum = 0.0;
2980 
2981  double* gfunc = DATA_PTR(m_gfunc_IJ);
2982  double* g2func = DATA_PTR(m_g2func_IJ);
2983  double* hfunc = DATA_PTR(m_hfunc_IJ);
2984  double* h2func = DATA_PTR(m_h2func_IJ);
2985  double* BMX_L = DATA_PTR(m_BMX_IJ_L);
2986  double* BprimeMX_L= DATA_PTR(m_BprimeMX_IJ_L);
2987  double* BphiMX_L = DATA_PTR(m_BphiMX_IJ_L);
2988  double* Phi_L = DATA_PTR(m_Phi_IJ_L);
2989  double* Phiprime = DATA_PTR(m_Phiprime_IJ);
2990  double* Phiphi_L = DATA_PTR(m_PhiPhi_IJ_L);
2991  double* CMX_L = DATA_PTR(m_CMX_IJ_L);
2992 
2993  double x1, x2;
2994  double dFdT, zsqdFdT;
2995  double sum1, sum2, sum3, sum4, sum5, term1;
2996  double sum_m_phi_minus_1, d_osmotic_coef_dT, d_lnwateract_dT;
2997 
2998  int z1, z2;
2999  size_t n, i, j, m, counterIJ, counterIJ2;
3000 
3001 #ifdef DEBUG_MODE
3002  if (m_debugCalc) {
3003  printf("\n Debugging information from "
3004  "s_Pitzer_dlnMolalityActCoeff_dT()\n");
3005  }
3006 #endif
3007  /*
3008  * Make sure the counter variables are setup
3009  */
3010  counterIJ_setup();
3011 
3012  /*
3013  * ---------- Calculate common sums over solutes ---------------------
3014  */
3015  for (n = 1; n < m_kk; n++) {
3016  // ionic strength
3017  Is += charge(n) * charge(n) * molality[n];
3018  // total molar charge
3019  molarcharge += fabs(charge(n)) * molality[n];
3020  molalitysum += molality[n];
3021  }
3022  Is *= 0.5;
3023 
3024  /*
3025  * Store the ionic molality in the object for reference.
3026  */
3027  m_IionicMolality = Is;
3028  sqrtIs = sqrt(Is);
3029 #ifdef DEBUG_MODE
3030  if (m_debugCalc) {
3031  printf(" Step 1: \n");
3032  printf(" ionic strenth = %14.7le \n total molar "
3033  "charge = %14.7le \n", Is, molarcharge);
3034  }
3035 #endif
3036 
3037  /*
3038  * The following call to calc_lambdas() calculates all 16 elements
3039  * of the elambda and elambda1 arrays, given the value of the
3040  * ionic strength (Is)
3041  */
3042  calc_lambdas(Is);
3043 
3044  /*
3045  * ----- Step 2: Find the coefficients E-theta and -------------------
3046  * E-thetaprime for all combinations of positive
3047  * unlike charges up to 4
3048  */
3049 #ifdef DEBUG_MODE
3050  if (m_debugCalc) {
3051  printf(" Step 2: \n");
3052  }
3053 #endif
3054  for (z1 = 1; z1 <=4; z1++) {
3055  for (z2 =1; z2 <=4; z2++) {
3056  calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
3057 #ifdef DEBUG_MODE
3058  if (m_debugCalc) {
3059  printf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n",
3060  z1, z2, etheta[z1][z2], etheta_prime[z1][z2]);
3061  }
3062 #endif
3063  }
3064  }
3065 
3066 #ifdef DEBUG_MODE
3067  if (m_debugCalc) {
3068  printf(" Step 3: \n");
3069  printf(" Species Species g(x) "
3070  " hfunc(x) \n");
3071  }
3072 #endif
3073 
3074  /*
3075  * calculate g(x) and hfunc(x) for each cation-anion pair MX
3076  * In the original literature, hfunc, was called gprime. However,
3077  * it's not the derivative of g(x), so I renamed it.
3078  */
3079  for (i = 1; i < (m_kk - 1); i++) {
3080  for (j = (i+1); j < m_kk; j++) {
3081  /*
3082  * Find the counterIJ for the symmetric binary interaction
3083  */
3084  n = m_kk*i + j;
3085  counterIJ = m_CounterIJ[n];
3086  /*
3087  * Only loop over oppositely charge species
3088  */
3089  if (charge(i)*charge(j) < 0) {
3090  /*
3091  * x is a reduced function variable
3092  */
3093  x1 = sqrtIs * alpha1MX[counterIJ];
3094  if (x1 > 1.0E-100) {
3095  gfunc[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 * x1);
3096  hfunc[counterIJ] = -2.0 *
3097  (1.0-(1.0 + x1 + 0.5 * x1 *x1) * exp(-x1)) / (x1 * x1);
3098  } else {
3099  gfunc[counterIJ] = 0.0;
3100  hfunc[counterIJ] = 0.0;
3101  }
3102 
3103  if (beta2MX_L[counterIJ] != 0.0) {
3104  x2 = sqrtIs * alpha2MX[counterIJ];
3105  if (x2 > 1.0E-100) {
3106  g2func[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
3107  h2func[counterIJ] = -2.0 *
3108  (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
3109  } else {
3110  g2func[counterIJ] = 0.0;
3111  h2func[counterIJ] = 0.0;
3112  }
3113  }
3114  } else {
3115  gfunc[counterIJ] = 0.0;
3116  hfunc[counterIJ] = 0.0;
3117  }
3118 #ifdef DEBUG_MODE
3119  if (m_debugCalc) {
3120  sni = speciesName(i);
3121  snj = speciesName(j);
3122  printf(" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(),
3123  gfunc[counterIJ], hfunc[counterIJ]);
3124  }
3125 #endif
3126  }
3127  }
3128 
3129  /*
3130  * ------- SUBSECTION TO CALCULATE BMX_L, BprimeMX_L, BphiMX_L ----------
3131  * ------- These are now temperature derivatives of the
3132  * previously calculated quantities.
3133  */
3134 #ifdef DEBUG_MODE
3135  if (m_debugCalc) {
3136  printf(" Step 4: \n");
3137  printf(" Species Species BMX "
3138  "BprimeMX BphiMX \n");
3139  }
3140 #endif
3141 
3142  for (i = 1; i < m_kk - 1; i++) {
3143  for (j = i+1; j < m_kk; j++) {
3144  /*
3145  * Find the counterIJ for the symmetric binary interaction
3146  */
3147  n = m_kk*i + j;
3148  counterIJ = m_CounterIJ[n];
3149  /*
3150  * both species have a non-zero charge, and one is positive
3151  * and the other is negative
3152  */
3153  if (charge(i)*charge(j) < 0.0) {
3154  BMX_L[counterIJ] = beta0MX_L[counterIJ]
3155  + beta1MX_L[counterIJ] * gfunc[counterIJ]
3156  + beta2MX_L[counterIJ] * gfunc[counterIJ];
3157 #ifdef DEBUG_MODE
3158  if (m_debugCalc) {
3159  printf("%d %g: %g %g %g %g\n",
3160  (int) counterIJ, BMX_L[counterIJ], beta0MX_L[counterIJ],
3161  beta1MX_L[counterIJ], beta2MX_L[counterIJ], gfunc[counterIJ]);
3162  }
3163 #endif
3164  if (Is > 1.0E-150) {
3165  BprimeMX_L[counterIJ] = (beta1MX_L[counterIJ] * hfunc[counterIJ]/Is +
3166  beta2MX_L[counterIJ] * h2func[counterIJ]/Is);
3167  } else {
3168  BprimeMX_L[counterIJ] = 0.0;
3169  }
3170  BphiMX_L[counterIJ] = BMX_L[counterIJ] + Is*BprimeMX_L[counterIJ];
3171  } else {
3172  BMX_L[counterIJ] = 0.0;
3173  BprimeMX_L[counterIJ] = 0.0;
3174  BphiMX_L[counterIJ] = 0.0;
3175  }
3176 #ifdef DEBUG_MODE
3177  if (m_debugCalc) {
3178  sni = speciesName(i);
3179  snj = speciesName(j);
3180  printf(" %-16s %-16s %11.7f %11.7f %11.7f \n",
3181  sni.c_str(), snj.c_str(),
3182  BMX_L[counterIJ], BprimeMX_L[counterIJ], BphiMX_L[counterIJ]);
3183  }
3184 #endif
3185  }
3186  }
3187 
3188  /*
3189  * --------- SUBSECTION TO CALCULATE CMX_L ----------
3190  * ---------
3191  */
3192 #ifdef DEBUG_MODE
3193  if (m_debugCalc) {
3194  printf(" Step 5: \n");
3195  printf(" Species Species CMX \n");
3196  }
3197 #endif
3198  for (i = 1; i < m_kk-1; i++) {
3199  for (j = i+1; j < m_kk; j++) {
3200  /*
3201  * Find the counterIJ for the symmetric binary interaction
3202  */
3203  n = m_kk*i + j;
3204  counterIJ = m_CounterIJ[n];
3205  /*
3206  * both species have a non-zero charge, and one is positive
3207  * and the other is negative
3208  */
3209  if (charge(i)*charge(j) < 0.0) {
3210  CMX_L[counterIJ] = CphiMX_L[counterIJ]/
3211  (2.0* sqrt(fabs(charge(i)*charge(j))));
3212  } else {
3213  CMX_L[counterIJ] = 0.0;
3214  }
3215 #ifdef DEBUG_MODE
3216  if (m_debugCalc) {
3217  sni = speciesName(i);
3218  snj = speciesName(j);
3219  printf(" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(),
3220  CMX_L[counterIJ]);
3221  }
3222 #endif
3223  }
3224  }
3225 
3226  /*
3227  * ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ----------
3228  * --------
3229  */
3230 #ifdef DEBUG_MODE
3231  if (m_debugCalc) {
3232  printf(" Step 6: \n");
3233  printf(" Species Species Phi_ij "
3234  " Phiprime_ij Phi^phi_ij \n");
3235  }
3236 #endif
3237  for (i = 1; i < m_kk-1; i++) {
3238  for (j = i+1; j < m_kk; j++) {
3239  /*
3240  * Find the counterIJ for the symmetric binary interaction
3241  */
3242  n = m_kk*i + j;
3243  counterIJ = m_CounterIJ[n];
3244  /*
3245  * both species have a non-zero charge, and one is positive
3246  * and the other is negative
3247  */
3248  if (charge(i)*charge(j) > 0) {
3249  z1 = (int) fabs(charge(i));
3250  z2 = (int) fabs(charge(j));
3251  //Phi[counterIJ] = thetaij_L[counterIJ] + etheta[z1][z2];
3252  Phi_L[counterIJ] = thetaij_L[counterIJ];
3253  //Phiprime[counterIJ] = etheta_prime[z1][z2];
3254  Phiprime[counterIJ] = 0.0;
3255  Phiphi_L[counterIJ] = Phi_L[counterIJ] + Is * Phiprime[counterIJ];
3256  } else {
3257  Phi_L[counterIJ] = 0.0;
3258  Phiprime[counterIJ] = 0.0;
3259  Phiphi_L[counterIJ] = 0.0;
3260  }
3261 #ifdef DEBUG_MODE
3262  if (m_debugCalc) {
3263  sni = speciesName(i);
3264  snj = speciesName(j);
3265  printf(" %-16s %-16s %10.6f %10.6f %10.6f \n",
3266  sni.c_str(), snj.c_str(),
3267  Phi_L[counterIJ], Phiprime[counterIJ], Phiphi_L[counterIJ]);
3268  }
3269 #endif
3270  }
3271  }
3272 
3273  /*
3274  * ----------- SUBSECTION FOR CALCULATION OF dFdT ---------------------
3275  */
3276 #ifdef DEBUG_MODE
3277  if (m_debugCalc) {
3278  printf(" Step 7: \n");
3279  }
3280 #endif
3281  // A_Debye_Huckel = 0.5092; (units = sqrt(kg/gmol))
3282  // A_Debye_Huckel = 0.5107; <- This value is used to match GWB data
3283  // ( A * ln(10) = 1.17593)
3284  // Aphi = A_Debye_Huckel * 2.30258509 / 3.0;
3285 
3286  double dA_DebyedT = dA_DebyedT_TP();
3287  double dAphidT = dA_DebyedT /3.0;
3288 #ifdef DEBUG_HKM
3289  //dAphidT = 0.0;
3290 #endif
3291  //F = -Aphi * ( sqrt(Is) / (1.0 + 1.2*sqrt(Is))
3292  // + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
3293  //dAphidT = Al / (4.0 * GasConstant * T * T);
3294  dFdT = -dAphidT * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
3295  + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
3296 #ifdef DEBUG_MODE
3297  if (m_debugCalc) {
3298  printf(" initial value of dFdT = %10.6f \n", dFdT);
3299  }
3300 #endif
3301  for (i = 1; i < m_kk-1; i++) {
3302  for (j = i+1; j < m_kk; j++) {
3303  /*
3304  * Find the counterIJ for the symmetric binary interaction
3305  */
3306  n = m_kk*i + j;
3307  counterIJ = m_CounterIJ[n];
3308  /*
3309  * both species have a non-zero charge, and one is positive
3310  * and the other is negative
3311  */
3312  if (charge(i)*charge(j) < 0) {
3313  dFdT = dFdT + molality[i]*molality[j] * BprimeMX_L[counterIJ];
3314  }
3315  /*
3316  * Both species have a non-zero charge, and they
3317  * have the same sign, e.g., both positive or both negative.
3318  */
3319  if (charge(i)*charge(j) > 0) {
3320  dFdT = dFdT + molality[i]*molality[j] * Phiprime[counterIJ];
3321  }
3322 #ifdef DEBUG_MODE
3323  if (m_debugCalc) {
3324  printf(" dFdT = %10.6f \n", dFdT);
3325  }
3326 #endif
3327  }
3328  }
3329 #ifdef DEBUG_MODE
3330  if (m_debugCalc) {
3331  printf(" Step 8: \n");
3332  }
3333 #endif
3334 
3335  for (i = 1; i < m_kk; i++) {
3336 
3337  /*
3338  * -------- SUBSECTION FOR CALCULATING THE dACTCOEFFdT FOR CATIONS -----
3339  * --
3340  */
3341  if (charge(i) > 0) {
3342  // species i is the cation (positive) to calc the actcoeff
3343  zsqdFdT = charge(i)*charge(i)*dFdT;
3344  sum1 = 0.0;
3345  sum2 = 0.0;
3346  sum3 = 0.0;
3347  sum4 = 0.0;
3348  sum5 = 0.0;
3349  for (j = 1; j < m_kk; j++) {
3350  /*
3351  * Find the counterIJ for the symmetric binary interaction
3352  */
3353  n = m_kk*i + j;
3354  counterIJ = m_CounterIJ[n];
3355 
3356  if (charge(j) < 0.0) {
3357  // sum over all anions
3358  sum1 = sum1 + molality[j]*
3359  (2.0*BMX_L[counterIJ] + molarcharge*CMX_L[counterIJ]);
3360  if (j < m_kk-1) {
3361  /*
3362  * This term is the ternary interaction involving the
3363  * non-duplicate sum over double anions, j, k, with
3364  * respect to the cation, i.
3365  */
3366  for (size_t k = j+1; k < m_kk; k++) {
3367  // an inner sum over all anions
3368  if (charge(k) < 0.0) {
3369  n = k + j * m_kk + i * m_kk * m_kk;
3370  sum3 = sum3 + molality[j]*molality[k]*psi_ijk_L[n];
3371  }
3372  }
3373  }
3374  }
3375 
3376 
3377  if (charge(j) > 0.0) {
3378  // sum over all cations
3379  if (j != i) {
3380  sum2 = sum2 + molality[j]*(2.0*Phi_L[counterIJ]);
3381  }
3382  for (size_t k = 1; k < m_kk; k++) {
3383  if (charge(k) < 0.0) {
3384  // two inner sums over anions
3385 
3386  n = k + j * m_kk + i * m_kk * m_kk;
3387  sum2 = sum2 + molality[j]*molality[k]*psi_ijk_L[n];
3388  /*
3389  * Find the counterIJ for the j,k interaction
3390  */
3391  n = m_kk*j + k;
3392  counterIJ2 = m_CounterIJ[n];
3393  sum4 = sum4 + (fabs(charge(i))*
3394  molality[j]*molality[k]*CMX_L[counterIJ2]);
3395  }
3396  }
3397  }
3398 
3399  /*
3400  * Handle neutral j species
3401  */
3402  if (charge(j) == 0) {
3403  sum5 = sum5 + molality[j]*2.0*m_Lambda_nj_L(j,i);
3404  }
3405  /*
3406  * Zeta interaction term
3407  */
3408  for (size_t k = 1; k < m_kk; k++) {
3409  if (charge(k) < 0.0) {
3410  size_t izeta = j;
3411  size_t jzeta = i;
3412  n = izeta * m_kk * m_kk + jzeta * m_kk + k;
3413  double zeta_L = psi_ijk_L[n];
3414  if (zeta_L != 0.0) {
3415  sum5 = sum5 + molality[j]*molality[k]*zeta_L;
3416  }
3417  }
3418  }
3419  }
3420  /*
3421  * Add all of the contributions up to yield the log of the
3422  * solute activity coefficients (molality scale)
3423  */
3425  zsqdFdT + sum1 + sum2 + sum3 + sum4 + sum5;
3426  d_gamma_dT_Unscaled[i] = exp(m_dlnActCoeffMolaldT_Unscaled[i]);
3427 #ifdef DEBUG_MODE
3428  if (m_debugCalc) {
3429  sni = speciesName(i);
3430  printf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f \n",
3431  sni.c_str(), m_dlnActCoeffMolaldT_Unscaled[i], d_gamma_dT_Unscaled[i]);
3432  printf(" %12g %12g %12g %12g %12g %12g\n",
3433  zsqdFdT, sum1, sum2, sum3, sum4, sum5);
3434  }
3435 #endif
3436  }
3437 
3438  /*
3439  * ------ SUBSECTION FOR CALCULATING THE dACTCOEFFdT FOR ANIONS ------
3440  *
3441  */
3442  if (charge(i) < 0) {
3443  // species i is an anion (negative)
3444  zsqdFdT = charge(i)*charge(i)*dFdT;
3445  sum1 = 0.0;
3446  sum2 = 0.0;
3447  sum3 = 0.0;
3448  sum4 = 0.0;
3449  sum5 = 0.0;
3450  for (j = 1; j < m_kk; j++) {
3451  /*
3452  * Find the counterIJ for the symmetric binary interaction
3453  */
3454  n = m_kk*i + j;
3455  counterIJ = m_CounterIJ[n];
3456 
3457  /*
3458  * For Anions, do the cation interactions.
3459  */
3460  if (charge(j) > 0) {
3461  sum1 = sum1 + molality[j]*
3462  (2.0*BMX_L[counterIJ] + molarcharge*CMX_L[counterIJ]);
3463  if (j < m_kk-1) {
3464  for (size_t k = j+1; k < m_kk; k++) {
3465  // an inner sum over all cations
3466  if (charge(k) > 0) {
3467  n = k + j * m_kk + i * m_kk * m_kk;
3468  sum3 = sum3 + molality[j]*molality[k]*psi_ijk_L[n];
3469  }
3470  }
3471  }
3472  }
3473 
3474  /*
3475  * For Anions, do the other anion interactions.
3476  */
3477  if (charge(j) < 0.0) {
3478  // sum over all anions
3479  if (j != i) {
3480  sum2 = sum2 + molality[j]*(2.0*Phi_L[counterIJ]);
3481  }
3482  for (size_t k = 1; k < m_kk; k++) {
3483  if (charge(k) > 0.0) {
3484  // two inner sums over cations
3485  n = k + j * m_kk + i * m_kk * m_kk;
3486  sum2 = sum2 + molality[j]*molality[k]*psi_ijk_L[n];
3487  /*
3488  * Find the counterIJ for the symmetric binary interaction
3489  */
3490  n = m_kk*j + k;
3491  counterIJ2 = m_CounterIJ[n];
3492  sum4 = sum4 +
3493  (fabs(charge(i))*
3494  molality[j]*molality[k]*CMX_L[counterIJ2]);
3495  }
3496  }
3497  }
3498 
3499  /*
3500  * for Anions, do the neutral species interaction
3501  */
3502  if (charge(j) == 0.0) {
3503  sum5 = sum5 + molality[j]*2.0*m_Lambda_nj_L(j,i);
3504  for (size_t k = 1; k < m_kk; k++) {
3505  if (charge(k) > 0.0) {
3506  size_t izeta = j;
3507  size_t jzeta = k;
3508  size_t kzeta = i;
3509  n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
3510  double zeta_L = psi_ijk_L[n];
3511  if (zeta_L != 0.0) {
3512  sum5 = sum5 + molality[j]*molality[k]*zeta_L;
3513  }
3514  }
3515  }
3516  }
3517  }
3519  zsqdFdT + sum1 + sum2 + sum3 + sum4 + sum5;
3520  d_gamma_dT_Unscaled[i] = exp(m_dlnActCoeffMolaldT_Unscaled[i]);
3521 #ifdef DEBUG_MODE
3522  if (m_debugCalc) {
3523  sni = speciesName(i);
3524  printf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f\n",
3525  sni.c_str(), m_dlnActCoeffMolaldT_Unscaled[i], d_gamma_dT_Unscaled[i]);
3526  printf(" %12g %12g %12g %12g %12g %12g\n",
3527  zsqdFdT, sum1, sum2, sum3, sum4, sum5);
3528  }
3529 #endif
3530  }
3531  /*
3532  * ------ SUBSECTION FOR CALCULATING NEUTRAL SOLUTE ACT COEFF -------
3533  * ------ -> equations agree with my notes,
3534  * -> Equations agree with Pitzer,
3535  */
3536  if (charge(i) == 0.0) {
3537  sum1 = 0.0;
3538  sum3 = 0.0;
3539  for (j = 1; j < m_kk; j++) {
3540  sum1 = sum1 + molality[j]*2.0*m_Lambda_nj_L(i,j);
3541  /*
3542  * Zeta term -> we piggyback on the psi term
3543  */
3544  if (charge(j) > 0.0) {
3545  for (size_t k = 1; k < m_kk; k++) {
3546  if (charge(k) < 0.0) {
3547  n = k + j * m_kk + i * m_kk * m_kk;
3548  sum3 = sum3 + molality[j]*molality[k]*psi_ijk_L[n];
3549  }
3550  }
3551  }
3552  }
3553  sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_L[i];
3554  m_dlnActCoeffMolaldT_Unscaled[i] = sum1 + sum2 + sum3;
3555  d_gamma_dT_Unscaled[i] = exp(m_dlnActCoeffMolaldT_Unscaled[i]);
3556 #ifdef DEBUG_MODE
3557  if (m_debugCalc) {
3558  sni = speciesName(i);
3559  printf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f \n",
3560  sni.c_str(), m_dlnActCoeffMolaldT_Unscaled[i], d_gamma_dT_Unscaled[i]);
3561  }
3562 #endif
3563  }
3564 
3565  }
3566 #ifdef DEBUG_MODE
3567  if (m_debugCalc) {
3568  printf(" Step 9: \n");
3569  }
3570 #endif
3571  /*
3572  * ------ SUBSECTION FOR CALCULATING THE d OSMOTIC COEFF dT ---------
3573  *
3574  */
3575  sum1 = 0.0;
3576  sum2 = 0.0;
3577  sum3 = 0.0;
3578  sum4 = 0.0;
3579  sum5 = 0.0;
3580  double sum6 = 0.0;
3581  double sum7 = 0.0;
3582  /*
3583  * term1 is the temperature derivative of the
3584  * DH term in the osmotic coefficient expression
3585  * b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer
3586  * implementations.
3587  * Is = Ionic strength on the molality scale (units of (gmol/kg))
3588  * Aphi = A_Debye / 3 (units of sqrt(kg/gmol))
3589  */
3590  term1 = -dAphidT * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
3591 
3592  for (j = 1; j < m_kk; j++) {
3593  /*
3594  * Loop Over Cations
3595  */
3596  if (charge(j) > 0.0) {
3597  for (size_t k = 1; k < m_kk; k++) {
3598  if (charge(k) < 0.0) {
3599  /*
3600  * Find the counterIJ for the symmetric j,k binary interaction
3601  */
3602  n = m_kk*j + k;
3603  counterIJ = m_CounterIJ[n];
3604 
3605  sum1 = sum1 + molality[j]*molality[k]*
3606  (BphiMX_L[counterIJ] + molarcharge*CMX_L[counterIJ]);
3607  }
3608  }
3609 
3610  for (size_t k = j+1; k < m_kk; k++) {
3611  if (j == (m_kk-1)) {
3612  // we should never reach this step
3613  printf("logic error 1 in Step 9 of hmw_act");
3614  exit(EXIT_FAILURE);
3615  }
3616  if (charge(k) > 0.0) {
3617  /*
3618  * Find the counterIJ for the symmetric j,k binary interaction
3619  * between 2 cations.
3620  */
3621  n = m_kk*j + k;
3622  counterIJ = m_CounterIJ[n];
3623  sum2 = sum2 + molality[j]*molality[k]*Phiphi_L[counterIJ];
3624  for (m = 1; m < m_kk; m++) {
3625  if (charge(m) < 0.0) {
3626  // species m is an anion
3627  n = m + k * m_kk + j * m_kk * m_kk;
3628  sum2 = sum2 +
3629  molality[j]*molality[k]*molality[m]*psi_ijk_L[n];
3630  }
3631  }
3632  }
3633  }
3634  }
3635 
3636  /*
3637  * Loop Over Anions
3638  */
3639  if (charge(j) < 0) {
3640  for (size_t k = j+1; k < m_kk; k++) {
3641  if (j == m_kk-1) {
3642  // we should never reach this step
3643  printf("logic error 2 in Step 9 of hmw_act");
3644  exit(EXIT_FAILURE);
3645  }
3646  if (charge(k) < 0) {
3647  /*
3648  * Find the counterIJ for the symmetric j,k binary interaction
3649  * between two anions
3650  */
3651  n = m_kk*j + k;
3652  counterIJ = m_CounterIJ[n];
3653 
3654  sum3 = sum3 + molality[j]*molality[k]*Phiphi_L[counterIJ];
3655  for (m = 1; m < m_kk; m++) {
3656  if (charge(m) > 0.0) {
3657  n = m + k * m_kk + j * m_kk * m_kk;
3658  sum3 = sum3 +
3659  molality[j]*molality[k]*molality[m]*psi_ijk_L[n];
3660  }
3661  }
3662  }
3663  }
3664  }
3665 
3666  /*
3667  * Loop Over Neutral Species
3668  */
3669  if (charge(j) == 0) {
3670  for (size_t k = 1; k < m_kk; k++) {
3671  if (charge(k) < 0.0) {
3672  sum4 = sum4 + molality[j]*molality[k]*m_Lambda_nj_L(j,k);
3673  }
3674  if (charge(k) > 0.0) {
3675  sum5 = sum5 + molality[j]*molality[k]*m_Lambda_nj_L(j,k);
3676  }
3677  if (charge(k) == 0.0) {
3678  if (k > j) {
3679  sum6 = sum6 + molality[j]*molality[k]*m_Lambda_nj_L(j,k);
3680  } else if (k == j) {
3681  sum6 = sum6 + 0.5 * molality[j]*molality[k]*m_Lambda_nj_L(j,k);
3682  }
3683  }
3684  if (charge(k) < 0.0) {
3685  size_t izeta = j;
3686  for (m = 1; m < m_kk; m++) {
3687  if (charge(m) > 0.0) {
3688  size_t jzeta = m;
3689  n = k + jzeta * m_kk + izeta * m_kk * m_kk;
3690  double zeta_L = psi_ijk_L[n];
3691  if (zeta_L != 0.0) {
3692  sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_L;
3693  }
3694  }
3695  }
3696  }
3697  }
3698  sum7 += molality[j]*molality[j]*molality[j]*m_Mu_nnn_L[j];
3699  }
3700  }
3701  sum_m_phi_minus_1 = 2.0 *
3702  (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
3703  /*
3704  * Calculate the osmotic coefficient from
3705  * osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i)
3706  */
3707  if (molalitysum > 1.0E-150) {
3708  d_osmotic_coef_dT = 0.0 + (sum_m_phi_minus_1 / molalitysum);
3709  } else {
3710  d_osmotic_coef_dT = 0.0;
3711  }
3712 #ifdef DEBUG_MODE
3713  if (m_debugCalc) {
3714  printf(" term1=%10.6f sum1=%10.6f sum2=%10.6f "
3715  "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n",
3716  term1, sum1, sum2, sum3, sum4, sum5);
3717  printf(" sum_m_phi_minus_1=%10.6f d_osmotic_coef_dT =%10.6f\n",
3718  sum_m_phi_minus_1, d_osmotic_coef_dT);
3719  }
3720 
3721  if (m_debugCalc) {
3722  printf(" Step 10: \n");
3723  }
3724 #endif
3725  d_lnwateract_dT = -(m_weightSolvent/1000.0) * molalitysum * d_osmotic_coef_dT;
3726 
3727  /*
3728  * In Cantera, we define the activity coefficient of the solvent as
3729  *
3730  * act_0 = actcoeff_0 * Xmol_0
3731  *
3732  * We have just computed act_0. However, this routine returns
3733  * ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
3734  */
3735  //double xmolSolvent = moleFraction(m_indexSolvent);
3736  m_dlnActCoeffMolaldT_Unscaled[0] = d_lnwateract_dT;
3737 #ifdef DEBUG_MODE
3738  if (m_debugCalc) {
3739  double d_wateract_dT = exp(d_lnwateract_dT);
3740  printf(" d_ln_a_water_dT = %10.6f d_a_water_dT=%10.6f\n\n",
3741  d_lnwateract_dT, d_wateract_dT);
3742  }
3743 #endif
3744 }
3745 
3747 {
3748  /*
3749  * Zero the unscaled 2nd derivatives
3750  */
3751  m_d2lnActCoeffMolaldT2_Unscaled.assign(m_kk, 0.0);
3752  /*
3753  * Calculate the unscaled 2nd derivatives
3754  */
3756 
3757  //double xmolSolvent = moleFraction(m_indexSolvent);
3758  //double xx = MAX(m_xmolSolventMIN, xmolSolvent);
3759  //double lnxs = log(xx);
3760 
3761  for (size_t k = 1; k < m_kk; k++) {
3762  if (CROP_speciesCropped_[k] == 2) {
3764  }
3765  }
3766 
3767  if (CROP_speciesCropped_[0]) {
3769  }
3770 
3771  /*
3772  * Scale the 2nd derivatives
3773  */
3775 }
3776 
3778 {
3779  /*
3780  * HKM -> Assumption is made that the solvent is species 0.
3781  */
3782 #ifdef DEBUG_MODE
3783  m_debugCalc = 0;
3784 #endif
3785  if (m_indexSolvent != 0) {
3786  printf("Wrong index solvent value!\n");
3787  exit(EXIT_FAILURE);
3788  }
3789 
3790  std::string sni, snj, snk;
3791 
3792  const double* molality = DATA_PTR(m_molalitiesCropped);
3793  const double* beta0MX_LL= DATA_PTR(m_Beta0MX_ij_LL);
3794  const double* beta1MX_LL= DATA_PTR(m_Beta1MX_ij_LL);
3795  const double* beta2MX_LL= DATA_PTR(m_Beta2MX_ij_LL);
3796  const double* CphiMX_LL = DATA_PTR(m_CphiMX_ij_LL);
3797  const double* thetaij_LL= DATA_PTR(m_Theta_ij_LL);
3798  const double* alpha1MX = DATA_PTR(m_Alpha1MX_ij);
3799  const double* alpha2MX = DATA_PTR(m_Alpha2MX_ij);
3800  const double* psi_ijk_LL= DATA_PTR(m_Psi_ijk_LL);
3801 
3802  /*
3803  * Local variables defined by Coltrin
3804  */
3805  double etheta[5][5], etheta_prime[5][5], sqrtIs;
3806  /*
3807  * Molality based ionic strength of the solution
3808  */
3809  double Is = 0.0;
3810  /*
3811  * Molarcharge of the solution: In Pitzer's notation,
3812  * this is his variable called "Z".
3813  */
3814  double molarcharge = 0.0;
3815  /*
3816  * molalitysum is the sum of the molalities over all solutes,
3817  * even those with zero charge.
3818  */
3819  double molalitysum = 0.0;
3820 
3821  double* gfunc = DATA_PTR(m_gfunc_IJ);
3822  double* g2func = DATA_PTR(m_g2func_IJ);
3823  double* hfunc = DATA_PTR(m_hfunc_IJ);
3824  double* h2func = DATA_PTR(m_h2func_IJ);
3825  double* BMX_LL = DATA_PTR(m_BMX_IJ_LL);
3826  double* BprimeMX_LL=DATA_PTR(m_BprimeMX_IJ_LL);
3827  double* BphiMX_LL= DATA_PTR(m_BphiMX_IJ_LL);
3828  double* Phi_LL = DATA_PTR(m_Phi_IJ_LL);
3829  double* Phiprime = DATA_PTR(m_Phiprime_IJ);
3830  double* Phiphi_LL= DATA_PTR(m_PhiPhi_IJ_LL);
3831  double* CMX_LL = DATA_PTR(m_CMX_IJ_LL);
3832 
3833 
3834  double x1, x2;
3835  double d2FdT2, zsqd2FdT2;
3836  double sum1, sum2, sum3, sum4, sum5, term1;
3837  double sum_m_phi_minus_1, d2_osmotic_coef_dT2, d2_lnwateract_dT2;
3838 
3839  int z1, z2;
3840  size_t n, i, j, m, counterIJ, counterIJ2;
3841 
3842 #ifdef DEBUG_MODE
3843  if (m_debugCalc) {
3844  printf("\n Debugging information from "
3845  "s_Pitzer_d2lnMolalityActCoeff_dT2()\n");
3846  }
3847 #endif
3848  /*
3849  * Make sure the counter variables are setup
3850  */
3851  counterIJ_setup();
3852 
3853 
3854  /*
3855  * ---------- Calculate common sums over solutes ---------------------
3856  */
3857  for (n = 1; n < m_kk; n++) {
3858  // ionic strength
3859  Is += charge(n) * charge(n) * molality[n];
3860  // total molar charge
3861  molarcharge += fabs(charge(n)) * molality[n];
3862  molalitysum += molality[n];
3863  }
3864  Is *= 0.5;
3865 
3866  /*
3867  * Store the ionic molality in the object for reference.
3868  */
3869  m_IionicMolality = Is;
3870  sqrtIs = sqrt(Is);
3871 #ifdef DEBUG_MODE
3872  if (m_debugCalc) {
3873  printf(" Step 1: \n");
3874  printf(" ionic strenth = %14.7le \n total molar "
3875  "charge = %14.7le \n", Is, molarcharge);
3876  }
3877 #endif
3878 
3879  /*
3880  * The following call to calc_lambdas() calculates all 16 elements
3881  * of the elambda and elambda1 arrays, given the value of the
3882  * ionic strength (Is)
3883  */
3884  calc_lambdas(Is);
3885 
3886  /*
3887  * ----- Step 2: Find the coefficients E-theta and -------------------
3888  * E-thetaprime for all combinations of positive
3889  * unlike charges up to 4
3890  */
3891 #ifdef DEBUG_MODE
3892  if (m_debugCalc) {
3893  printf(" Step 2: \n");
3894  }
3895 #endif
3896  for (z1 = 1; z1 <=4; z1++) {
3897  for (z2 =1; z2 <=4; z2++) {
3898  calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
3899 #ifdef DEBUG_MODE
3900  if (m_debugCalc) {
3901  printf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n",
3902  z1, z2, etheta[z1][z2], etheta_prime[z1][z2]);
3903  }
3904 #endif
3905  }
3906  }
3907 
3908 #ifdef DEBUG_MODE
3909  if (m_debugCalc) {
3910  printf(" Step 3: \n");
3911  printf(" Species Species g(x) "
3912  " hfunc(x) \n");
3913  }
3914 #endif
3915 
3916  /*
3917  *
3918  * calculate gfunc(x) and hfunc(x) for each cation-anion pair MX
3919  * In the original literature, hfunc, was called gprime. However,
3920  * it's not the derivative of gfunc(x), so I renamed it.
3921  */
3922  for (i = 1; i < (m_kk - 1); i++) {
3923  for (j = (i+1); j < m_kk; j++) {
3924  /*
3925  * Find the counterIJ for the symmetric binary interaction
3926  */
3927  n = m_kk*i + j;
3928  counterIJ = m_CounterIJ[n];
3929  /*
3930  * Only loop over oppositely charge species
3931  */
3932  if (charge(i)*charge(j) < 0) {
3933  /*
3934  * x is a reduced function variable
3935  */
3936  x1 = sqrtIs * alpha1MX[counterIJ];
3937  if (x1 > 1.0E-100) {
3938  gfunc[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 *x1);
3939  hfunc[counterIJ] = -2.0*
3940  (1.0-(1.0 + x1 + 0.5*x1 * x1) * exp(-x1)) / (x1 * x1);
3941  } else {
3942  gfunc[counterIJ] = 0.0;
3943  hfunc[counterIJ] = 0.0;
3944  }
3945 
3946  if (beta2MX_LL[counterIJ] != 0.0) {
3947  x2 = sqrtIs * alpha2MX[counterIJ];
3948  if (x2 > 1.0E-100) {
3949  g2func[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
3950  h2func[counterIJ] = -2.0 *
3951  (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
3952  } else {
3953  g2func[counterIJ] = 0.0;
3954  h2func[counterIJ] = 0.0;
3955  }
3956  }
3957  } else {
3958  gfunc[counterIJ] = 0.0;
3959  hfunc[counterIJ] = 0.0;
3960  }
3961 #ifdef DEBUG_MODE
3962  if (m_debugCalc) {
3963  sni = speciesName(i);
3964  snj = speciesName(j);
3965  printf(" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(),
3966  gfunc[counterIJ], hfunc[counterIJ]);
3967  }
3968 #endif
3969  }
3970  }
3971  /*
3972  * ------- SUBSECTION TO CALCULATE BMX_L, BprimeMX_LL, BphiMX_L ----------
3973  * ------- These are now temperature derivatives of the
3974  * previously calculated quantities.
3975  */
3976 #ifdef DEBUG_MODE
3977  if (m_debugCalc) {
3978  printf(" Step 4: \n");
3979  printf(" Species Species BMX "
3980  "BprimeMX BphiMX \n");
3981  }
3982 #endif
3983 
3984  for (i = 1; i < m_kk - 1; i++) {
3985  for (j = i+1; j < m_kk; j++) {
3986  /*
3987  * Find the counterIJ for the symmetric binary interaction
3988  */
3989  n = m_kk*i + j;
3990  counterIJ = m_CounterIJ[n];
3991  /*
3992  * both species have a non-zero charge, and one is positive
3993  * and the other is negative
3994  */
3995  if (charge(i)*charge(j) < 0.0) {
3996  BMX_LL[counterIJ] = beta0MX_LL[counterIJ]
3997  + beta1MX_LL[counterIJ] * gfunc[counterIJ]
3998  + beta2MX_LL[counterIJ] * g2func[counterIJ];
3999 #ifdef DEBUG_MODE
4000  if (m_debugCalc) {
4001  printf("%d %g: %g %g %g %g\n",
4002  (int) counterIJ, BMX_LL[counterIJ], beta0MX_LL[counterIJ],
4003  beta1MX_LL[counterIJ], beta2MX_LL[counterIJ], gfunc[counterIJ]);
4004  }
4005 #endif
4006  if (Is > 1.0E-150) {
4007  BprimeMX_LL[counterIJ] = (beta1MX_LL[counterIJ] * hfunc[counterIJ]/Is +
4008  beta2MX_LL[counterIJ] * h2func[counterIJ]/Is);
4009  } else {
4010  BprimeMX_LL[counterIJ] = 0.0;
4011  }
4012  BphiMX_LL[counterIJ] = BMX_LL[counterIJ] + Is*BprimeMX_LL[counterIJ];
4013  } else {
4014  BMX_LL[counterIJ] = 0.0;
4015  BprimeMX_LL[counterIJ] = 0.0;
4016  BphiMX_LL[counterIJ] = 0.0;
4017  }
4018 #ifdef DEBUG_MODE
4019  if (m_debugCalc) {
4020  sni = speciesName(i);
4021  snj = speciesName(j);
4022  printf(" %-16s %-16s %11.7f %11.7f %11.7f \n",
4023  sni.c_str(), snj.c_str(),
4024  BMX_LL[counterIJ], BprimeMX_LL[counterIJ], BphiMX_LL[counterIJ]);
4025  }
4026 #endif
4027  }
4028  }
4029 
4030  /*
4031  * --------- SUBSECTION TO CALCULATE CMX_LL ----------
4032  * ---------
4033  */
4034 #ifdef DEBUG_MODE
4035  if (m_debugCalc) {
4036  printf(" Step 5: \n");
4037  printf(" Species Species CMX \n");
4038  }
4039 #endif
4040  for (i = 1; i < m_kk-1; i++) {
4041  for (j = i+1; j < m_kk; j++) {
4042  /*
4043  * Find the counterIJ for the symmetric binary interaction
4044  */
4045  n = m_kk*i + j;
4046  counterIJ = m_CounterIJ[n];
4047  /*
4048  * both species have a non-zero charge, and one is positive
4049  * and the other is negative
4050  */
4051  if (charge(i)*charge(j) < 0.0) {
4052  CMX_LL[counterIJ] = CphiMX_LL[counterIJ]/
4053  (2.0* sqrt(fabs(charge(i)*charge(j))));
4054  } else {
4055  CMX_LL[counterIJ] = 0.0;
4056  }
4057 #ifdef DEBUG_MODE
4058  if (m_debugCalc) {
4059  sni = speciesName(i);
4060  snj = speciesName(j);
4061  printf(" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(),
4062  CMX_LL[counterIJ]);
4063  }
4064 #endif
4065  }
4066  }
4067 
4068  /*
4069  * ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ----------
4070  * --------
4071  */
4072 #ifdef DEBUG_MODE
4073  if (m_debugCalc) {
4074  printf(" Step 6: \n");
4075  printf(" Species Species Phi_ij "
4076  " Phiprime_ij Phi^phi_ij \n");
4077  }
4078 #endif
4079  for (i = 1; i < m_kk-1; i++) {
4080  for (j = i+1; j < m_kk; j++) {
4081  /*
4082  * Find the counterIJ for the symmetric binary interaction
4083  */
4084  n = m_kk*i + j;
4085  counterIJ = m_CounterIJ[n];
4086  /*
4087  * both species have a non-zero charge, and one is positive
4088  * and the other is negative
4089  */
4090  if (charge(i)*charge(j) > 0) {
4091  z1 = (int) fabs(charge(i));
4092  z2 = (int) fabs(charge(j));
4093  //Phi[counterIJ] = thetaij[counterIJ] + etheta[z1][z2];
4094  //Phi_L[counterIJ] = thetaij_L[counterIJ];
4095  Phi_LL[counterIJ] = thetaij_LL[counterIJ];
4096  //Phiprime[counterIJ] = etheta_prime[z1][z2];
4097  Phiprime[counterIJ] = 0.0;
4098  //Phiphi[counterIJ] = Phi[counterIJ] + Is * Phiprime[counterIJ];
4099  //Phiphi_L[counterIJ] = Phi_L[counterIJ] + Is * Phiprime[counterIJ];
4100  Phiphi_LL[counterIJ] = Phi_LL[counterIJ];
4101  } else {
4102  Phi_LL[counterIJ] = 0.0;
4103  Phiprime[counterIJ] = 0.0;
4104  Phiphi_LL[counterIJ] = 0.0;
4105  }
4106 #ifdef DEBUG_MODE
4107  if (m_debugCalc) {
4108  sni = speciesName(i);
4109  snj = speciesName(j);
4110  //printf(" %-16s %-16s %10.6f %10.6f %10.6f \n",
4111  // sni.c_str(), snj.c_str(),
4112  // Phi_L[counterIJ], Phiprime[counterIJ], Phiphi_L[counterIJ] );
4113  printf(" %-16s %-16s %10.6f %10.6f %10.6f \n",
4114  sni.c_str(), snj.c_str(),
4115  Phi_LL[counterIJ], Phiprime[counterIJ], Phiphi_LL[counterIJ]);
4116  }
4117 #endif
4118  }
4119  }
4120 
4121  /*
4122  * ----------- SUBSECTION FOR CALCULATION OF d2FdT2 ---------------------
4123  */
4124 #ifdef DEBUG_MODE
4125  if (m_debugCalc) {
4126  printf(" Step 7: \n");
4127  }
4128 #endif
4129  // A_Debye_Huckel = 0.5092; (units = sqrt(kg/gmol))
4130  // A_Debye_Huckel = 0.5107; <- This value is used to match GWB data
4131  // ( A * ln(10) = 1.17593)
4132  // Aphi = A_Debye_Huckel * 2.30258509 / 3.0;
4133  // Aphi = m_A_Debye / 3.0;
4134 
4135  //double dA_DebyedT = dA_DebyedT_TP();
4136  //double dAphidT = dA_DebyedT /3.0;
4137  double d2AphidT2 = d2A_DebyedT2_TP() / 3.0;
4138 #ifdef DEBUG_HKM
4139  //d2AphidT2 = 0.0;
4140 #endif
4141  //F = -Aphi * ( sqrt(Is) / (1.0 + 1.2*sqrt(Is))
4142  // + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
4143  //dAphidT = Al / (4.0 * GasConstant * T * T);
4144  //dFdT = -dAphidT * ( sqrt(Is) / (1.0 + 1.2*sqrt(Is))
4145  // + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
4146  d2FdT2 = -d2AphidT2 * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
4147  + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
4148 #ifdef DEBUG_MODE
4149  if (m_debugCalc) {
4150  printf(" initial value of d2FdT2 = %10.6f \n", d2FdT2);
4151  }
4152 #endif
4153  for (i = 1; i < m_kk-1; i++) {
4154  for (j = i+1; j < m_kk; j++) {
4155  /*
4156  * Find the counterIJ for the symmetric binary interaction
4157  */
4158  n = m_kk*i + j;
4159  counterIJ = m_CounterIJ[n];
4160  /*
4161  * both species have a non-zero charge, and one is positive
4162  * and the other is negative
4163  */
4164  if (charge(i)*charge(j) < 0) {
4165  d2FdT2 = d2FdT2 + molality[i]*molality[j] * BprimeMX_LL[counterIJ];
4166  }
4167  /*
4168  * Both species have a non-zero charge, and they
4169  * have the same sign, e.g., both positive or both negative.
4170  */
4171  if (charge(i)*charge(j) > 0) {
4172  d2FdT2 = d2FdT2 + molality[i]*molality[j] * Phiprime[counterIJ];
4173  }
4174 #ifdef DEBUG_MODE
4175  if (m_debugCalc) {
4176  printf(" d2FdT2 = %10.6f \n", d2FdT2);
4177  }
4178 #endif
4179  }
4180  }
4181 #ifdef DEBUG_MODE
4182  if (m_debugCalc) {
4183  printf(" Step 8: \n");
4184  }
4185 #endif
4186 
4187  for (i = 1; i < m_kk; i++) {
4188 
4189  /*
4190  * -------- SUBSECTION FOR CALCULATING THE dACTCOEFFdT FOR CATIONS -----
4191  * --
4192  */
4193  if (charge(i) > 0) {
4194  // species i is the cation (positive) to calc the actcoeff
4195  zsqd2FdT2 = charge(i)*charge(i)*d2FdT2;
4196  sum1 = 0.0;
4197  sum2 = 0.0;
4198  sum3 = 0.0;
4199  sum4 = 0.0;
4200  sum5 = 0.0;
4201  for (j = 1; j < m_kk; j++) {
4202  /*
4203  * Find the counterIJ for the symmetric binary interaction
4204  */
4205  n = m_kk*i + j;
4206  counterIJ = m_CounterIJ[n];
4207 
4208  if (charge(j) < 0.0) {
4209  // sum over all anions
4210  sum1 = sum1 + molality[j]*
4211  (2.0*BMX_LL[counterIJ] + molarcharge*CMX_LL[counterIJ]);
4212  if (j < m_kk-1) {
4213  /*
4214  * This term is the ternary interaction involving the
4215  * non-duplicate sum over double anions, j, k, with
4216  * respect to the cation, i.
4217  */
4218  for (size_t k = j+1; k < m_kk; k++) {
4219  // an inner sum over all anions
4220  if (charge(k) < 0.0) {
4221  n = k + j * m_kk + i * m_kk * m_kk;
4222  sum3 = sum3 + molality[j]*molality[k]*psi_ijk_LL[n];
4223  }
4224  }
4225  }
4226  }
4227 
4228 
4229  if (charge(j) > 0.0) {
4230  // sum over all cations
4231  if (j != i) {
4232  sum2 = sum2 + molality[j]*(2.0*Phi_LL[counterIJ]);
4233  }
4234  for (size_t k = 1; k < m_kk; k++) {
4235  if (charge(k) < 0.0) {
4236  // two inner sums over anions
4237 
4238  n = k + j * m_kk + i * m_kk * m_kk;
4239  sum2 = sum2 + molality[j]*molality[k]*psi_ijk_LL[n];
4240  /*
4241  * Find the counterIJ for the j,k interaction
4242  */
4243  n = m_kk*j + k;
4244  counterIJ2 = m_CounterIJ[n];
4245  sum4 = sum4 + (fabs(charge(i))*
4246  molality[j]*molality[k]*CMX_LL[counterIJ2]);
4247  }
4248  }
4249  }
4250 
4251  /*
4252  * Handle neutral j species
4253  */
4254  if (charge(j) == 0) {
4255  sum5 = sum5 + molality[j]*2.0*m_Lambda_nj_LL(j,i);
4256  /*
4257  * Zeta interaction term
4258  */
4259  for (size_t k = 1; k < m_kk; k++) {
4260  if (charge(k) < 0.0) {
4261  size_t izeta = j;
4262  size_t jzeta = i;
4263  n = izeta * m_kk * m_kk + jzeta * m_kk + k;
4264  double zeta_LL = psi_ijk_LL[n];
4265  if (zeta_LL != 0.0) {
4266  sum5 = sum5 + molality[j]*molality[k]*zeta_LL;
4267  }
4268  }
4269  }
4270  }
4271  }
4272  /*
4273  * Add all of the contributions up to yield the log of the
4274  * solute activity coefficients (molality scale)
4275  */
4277  zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5;
4278 #ifdef DEBUG_MODE
4279  if (m_debugCalc) {
4280  sni = speciesName(i);
4281  printf(" %-16s d2lngammadT2[i]=%10.6f \n",
4282  sni.c_str(), m_d2lnActCoeffMolaldT2_Unscaled[i]);
4283  printf(" %12g %12g %12g %12g %12g %12g\n",
4284  zsqd2FdT2, sum1, sum2, sum3, sum4, sum5);
4285  }
4286 #endif
4287  }
4288 
4289 
4290  /*
4291  * ------ SUBSECTION FOR CALCULATING THE d2ACTCOEFFdT2 FOR ANIONS ------
4292  *
4293  */
4294  if (charge(i) < 0) {
4295  // species i is an anion (negative)
4296  zsqd2FdT2 = charge(i)*charge(i)*d2FdT2;
4297  sum1 = 0.0;
4298  sum2 = 0.0;
4299  sum3 = 0.0;
4300  sum4 = 0.0;
4301  sum5 = 0.0;
4302  for (j = 1; j < m_kk; j++) {
4303  /*
4304  * Find the counterIJ for the symmetric binary interaction
4305  */
4306  n = m_kk*i + j;
4307  counterIJ = m_CounterIJ[n];
4308 
4309  /*
4310  * For Anions, do the cation interactions.
4311  */
4312  if (charge(j) > 0) {
4313  sum1 = sum1 + molality[j]*
4314  (2.0*BMX_LL[counterIJ] + molarcharge*CMX_LL[counterIJ]);
4315  if (j < m_kk-1) {
4316  for (size_t k = j+1; k < m_kk; k++) {
4317  // an inner sum over all cations
4318  if (charge(k) > 0) {
4319  n = k + j * m_kk + i * m_kk * m_kk;
4320  sum3 = sum3 + molality[j]*molality[k]*psi_ijk_LL[n];
4321  }
4322  }
4323  }
4324  }
4325 
4326  /*
4327  * For Anions, do the other anion interactions.
4328  */
4329  if (charge(j) < 0.0) {
4330  // sum over all anions
4331  if (j != i) {
4332  sum2 = sum2 + molality[j]*(2.0*Phi_LL[counterIJ]);
4333  }
4334  for (size_t k = 1; k < m_kk; k++) {
4335  if (charge(k) > 0.0) {
4336  // two inner sums over cations
4337  n = k + j * m_kk + i * m_kk * m_kk;
4338  sum2 = sum2 + molality[j]*molality[k]*psi_ijk_LL[n];
4339  /*
4340  * Find the counterIJ for the symmetric binary interaction
4341  */
4342  n = m_kk*j + k;
4343  counterIJ2 = m_CounterIJ[n];
4344  sum4 = sum4 +
4345  (fabs(charge(i))*
4346  molality[j]*molality[k]*CMX_LL[counterIJ2]);
4347  }
4348  }
4349  }
4350 
4351  /*
4352  * for Anions, do the neutral species interaction
4353  */
4354  if (charge(j) == 0.0) {
4355  sum5 = sum5 + molality[j]*2.0*m_Lambda_nj_LL(j,i);
4356  /*
4357  * Zeta interaction term
4358  */
4359  for (size_t k = 1; k < m_kk; k++) {
4360  if (charge(k) > 0.0) {
4361  size_t izeta = j;
4362  size_t jzeta = k;
4363  size_t kzeta = i;
4364  n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
4365  double zeta_LL = psi_ijk_LL[n];
4366  if (zeta_LL != 0.0) {
4367  sum5 = sum5 + molality[j]*molality[k]*zeta_LL;
4368  }
4369  }
4370  }
4371  }
4372  }
4374  zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5;
4375 #ifdef DEBUG_MODE
4376  if (m_debugCalc) {
4377  sni = speciesName(i);
4378  printf(" %-16s d2lngammadT2[i]=%10.6f\n",
4379  sni.c_str(), m_d2lnActCoeffMolaldT2_Unscaled[i]);
4380  printf(" %12g %12g %12g %12g %12g %12g\n",
4381  zsqd2FdT2, sum1, sum2, sum3, sum4, sum5);
4382  }
4383 #endif
4384  }
4385  /*
4386  * ------ SUBSECTION FOR CALCULATING NEUTRAL SOLUTE ACT COEFF -------
4387  * ------ -> equations agree with my notes,
4388  * -> Equations agree with Pitzer,
4389  */
4390  if (charge(i) == 0.0) {
4391  sum1 = 0.0;
4392  sum3 = 0.0;
4393  for (j = 1; j < m_kk; j++) {
4394  sum1 = sum1 + molality[j]*2.0*m_Lambda_nj_LL(i,j);
4395  /*
4396  * Zeta term -> we piggyback on the psi term
4397  */
4398  if (charge(j) > 0.0) {
4399  for (size_t k = 1; k < m_kk; k++) {
4400  if (charge(k) < 0.0) {
4401  n = k + j * m_kk + i * m_kk * m_kk;
4402  sum3 = sum3 + molality[j]*molality[k]*psi_ijk_LL[n];
4403  }
4404  }
4405  }
4406  }
4407  sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_LL[i];
4408  m_d2lnActCoeffMolaldT2_Unscaled[i] = sum1 + sum2 + sum3;
4409 #ifdef DEBUG_MODE
4410  if (m_debugCalc) {
4411  sni = speciesName(i);
4412  printf(" %-16s d2lngammadT2[i]=%10.6f \n",
4413  sni.c_str(), m_d2lnActCoeffMolaldT2_Unscaled[i]);
4414  }
4415 #endif
4416  }
4417 
4418  }
4419 #ifdef DEBUG_MODE
4420  if (m_debugCalc) {
4421  printf(" Step 9: \n");
4422  }
4423 #endif
4424 
4425  /*
4426  * ------ SUBSECTION FOR CALCULATING THE d2 OSMOTIC COEFF dT2 ---------
4427  */
4428  sum1 = 0.0;
4429  sum2 = 0.0;
4430  sum3 = 0.0;
4431  sum4 = 0.0;
4432  sum5 = 0.0;
4433  double sum6 = 0.0;
4434  double sum7 = 0.0;
4435  /*
4436  * term1 is the temperature derivative of the
4437  * DH term in the osmotic coefficient expression
4438  * b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer
4439  * implementations.
4440  * Is = Ionic strength on the molality scale (units of (gmol/kg))
4441  * Aphi = A_Debye / 3 (units of sqrt(kg/gmol))
4442  */
4443  term1 = -d2AphidT2 * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
4444 
4445  for (j = 1; j < m_kk; j++) {
4446  /*
4447  * Loop Over Cations
4448  */
4449  if (charge(j) > 0.0) {
4450  for (size_t k = 1; k < m_kk; k++) {
4451  if (charge(k) < 0.0) {
4452  /*
4453  * Find the counterIJ for the symmetric j,k binary interaction
4454  */
4455  n = m_kk*j + k;
4456  counterIJ = m_CounterIJ[n];
4457 
4458  sum1 = sum1 + molality[j]*molality[k]*
4459  (BphiMX_LL[counterIJ] + molarcharge*CMX_LL[counterIJ]);
4460  }
4461  }
4462 
4463  for (size_t k = j+1; k < m_kk; k++) {
4464  if (j == (m_kk-1)) {
4465  // we should never reach this step
4466  printf("logic error 1 in Step 9 of hmw_act");
4467  exit(EXIT_FAILURE);
4468  }
4469  if (charge(k) > 0.0) {
4470  /*
4471  * Find the counterIJ for the symmetric j,k binary interaction
4472  * between 2 cations.
4473  */
4474  n = m_kk*j + k;
4475  counterIJ = m_CounterIJ[n];
4476  sum2 = sum2 + molality[j]*molality[k]*Phiphi_LL[counterIJ];
4477  for (m = 1; m < m_kk; m++) {
4478  if (charge(m) < 0.0) {
4479  // species m is an anion
4480  n = m + k * m_kk + j * m_kk * m_kk;
4481  sum2 = sum2 +
4482  molality[j]*molality[k]*molality[m]*psi_ijk_LL[n];
4483  }
4484  }
4485  }
4486  }
4487  }
4488 
4489  /*
4490  * Loop Over Anions
4491  */
4492  if (charge(j) < 0) {
4493  for (size_t k = j+1; k < m_kk; k++) {
4494  if (j == m_kk-1) {
4495  // we should never reach this step
4496  printf("logic error 2 in Step 9 of hmw_act");
4497  exit(EXIT_FAILURE);
4498  }
4499  if (charge(k) < 0) {
4500  /*
4501  * Find the counterIJ for the symmetric j,k binary interaction
4502  * between two anions
4503  */
4504  n = m_kk*j + k;
4505  counterIJ = m_CounterIJ[n];
4506 
4507  sum3 = sum3 + molality[j]*molality[k]*Phiphi_LL[counterIJ];
4508  for (m = 1; m < m_kk; m++) {
4509  if (charge(m) > 0.0) {
4510  n = m + k * m_kk + j * m_kk * m_kk;
4511  sum3 = sum3 +
4512  molality[j]*molality[k]*molality[m]*psi_ijk_LL[n];
4513  }
4514  }
4515  }
4516  }
4517  }
4518 
4519  /*
4520  * Loop Over Neutral Species
4521  */
4522  if (charge(j) == 0) {
4523  for (size_t k = 1; k < m_kk; k++) {
4524  if (charge(k) < 0.0) {
4525  sum4 = sum4 + molality[j]*molality[k]*m_Lambda_nj_LL(j,k);
4526  }
4527  if (charge(k) > 0.0) {
4528  sum5 = sum5 + molality[j]*molality[k]*m_Lambda_nj_LL(j,k);
4529  }
4530  if (charge(k) == 0.0) {
4531  if (k > j) {
4532  sum6 = sum6 + molality[j]*molality[k]*m_Lambda_nj_LL(j,k);
4533  } else if (k == j) {
4534  sum6 = sum6 + 0.5 * molality[j]*molality[k]*m_Lambda_nj_LL(j,k);
4535  }
4536  }
4537  if (charge(k) < 0.0) {
4538  size_t izeta = j;
4539  for (m = 1; m < m_kk; m++) {
4540  if (charge(m) > 0.0) {
4541  size_t jzeta = m;
4542  n = k + jzeta * m_kk + izeta * m_kk * m_kk;
4543  double zeta_LL = psi_ijk_LL[n];
4544  if (zeta_LL != 0.0) {
4545  sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_LL;
4546  }
4547  }
4548  }
4549  }
4550  }
4551 
4552  sum7 += molality[j] * molality[j] * molality[j] * m_Mu_nnn_LL[j];
4553  }
4554  }
4555  sum_m_phi_minus_1 = 2.0 *
4556  (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
4557  /*
4558  * Calculate the osmotic coefficient from
4559  * osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i)
4560  */
4561  if (molalitysum > 1.0E-150) {
4562  d2_osmotic_coef_dT2 = 0.0 + (sum_m_phi_minus_1 / molalitysum);
4563  } else {
4564  d2_osmotic_coef_dT2 = 0.0;
4565  }
4566 #ifdef DEBUG_MODE
4567  if (m_debugCalc) {
4568  printf(" term1=%10.6f sum1=%10.6f sum2=%10.6f "
4569  "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n",
4570  term1, sum1, sum2, sum3, sum4, sum5);
4571  printf(" sum_m_phi_minus_1=%10.6f d2_osmotic_coef_dT2=%10.6f\n",
4572  sum_m_phi_minus_1, d2_osmotic_coef_dT2);
4573  }
4574 
4575  if (m_debugCalc) {
4576  printf(" Step 10: \n");
4577  }
4578 #endif
4579  d2_lnwateract_dT2 = -(m_weightSolvent/1000.0) * molalitysum * d2_osmotic_coef_dT2;
4580 
4581  /*
4582  * In Cantera, we define the activity coefficient of the solvent as
4583  *
4584  * act_0 = actcoeff_0 * Xmol_0
4585  *
4586  * We have just computed act_0. However, this routine returns
4587  * ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
4588  */
4589  m_d2lnActCoeffMolaldT2_Unscaled[0] = d2_lnwateract_dT2;
4590 
4591 #ifdef DEBUG_MODE
4592  if (m_debugCalc) {
4593  double d2_wateract_dT2 = exp(d2_lnwateract_dT2);
4594  printf(" d2_ln_a_water_dT2 = %10.6f d2_a_water_dT2=%10.6f\n\n",
4595  d2_lnwateract_dT2, d2_wateract_dT2);
4596  }
4597 #endif
4598 }
4599 
4601 {
4602  m_dlnActCoeffMolaldP_Unscaled.assign(m_kk, 0.0);
4604 
4605  for (size_t k = 1; k < m_kk; k++) {
4606  if (CROP_speciesCropped_[k] == 2) {
4608  }
4609  }
4610 
4611  if (CROP_speciesCropped_[0]) {
4613  }
4614 
4616 }
4617 
4619 {
4620  /*
4621  * HKM -> Assumption is made that the solvent is species 0.
4622  */
4623 #ifdef DEBUG_MODE
4624  m_debugCalc = 0;
4625 #endif
4626  if (m_indexSolvent != 0) {
4627  printf("Wrong index solvent value!\n");
4628  exit(EXIT_FAILURE);
4629  }
4630 
4631  std::string sni, snj, snk;
4632 
4633  const double* molality = DATA_PTR(m_molalitiesCropped);
4634  const double* beta0MX_P = DATA_PTR(m_Beta0MX_ij_P);
4635  const double* beta1MX_P = DATA_PTR(m_Beta1MX_ij_P);
4636  const double* beta2MX_P = DATA_PTR(m_Beta2MX_ij_P);
4637  const double* CphiMX_P = DATA_PTR(m_CphiMX_ij_P);
4638  const double* thetaij_P = DATA_PTR(m_Theta_ij_P);
4639  const double* alpha1MX = DATA_PTR(m_Alpha1MX_ij);
4640  const double* alpha2MX = DATA_PTR(m_Alpha2MX_ij);
4641  const double* psi_ijk_P = DATA_PTR(m_Psi_ijk_P);
4642 
4643  /*
4644  * Local variables defined by Coltrin
4645  */
4646  double etheta[5][5], etheta_prime[5][5], sqrtIs;
4647  /*
4648  * Molality based ionic strength of the solution
4649  */
4650  double Is = 0.0;
4651  /*
4652  * Molarcharge of the solution: In Pitzer's notation,
4653  * this is his variable called "Z".
4654  */
4655  double molarcharge = 0.0;
4656  /*
4657  * molalitysum is the sum of the molalities over all solutes,
4658  * even those with zero charge.
4659  */
4660  double molalitysum = 0.0;
4661 
4662  double* gfunc = DATA_PTR(m_gfunc_IJ);
4663  double* g2func = DATA_PTR(m_g2func_IJ);
4664  double* hfunc = DATA_PTR(m_hfunc_IJ);
4665  double* h2func = DATA_PTR(m_h2func_IJ);
4666  double* BMX_P = DATA_PTR(m_BMX_IJ_P);
4667  double* BprimeMX_P= DATA_PTR(m_BprimeMX_IJ_P);
4668  double* BphiMX_P = DATA_PTR(m_BphiMX_IJ_P);
4669  double* Phi_P = DATA_PTR(m_Phi_IJ_P);
4670  double* Phiprime = DATA_PTR(m_Phiprime_IJ);
4671  double* Phiphi_P = DATA_PTR(m_PhiPhi_IJ_P);
4672  double* CMX_P = DATA_PTR(m_CMX_IJ_P);
4673 
4674  double x1, x2;
4675  double dFdP, zsqdFdP;
4676  double sum1, sum2, sum3, sum4, sum5, term1;
4677  double sum_m_phi_minus_1, d_osmotic_coef_dP, d_lnwateract_dP;
4678 
4679  int z1, z2;
4680  size_t n, i, j, m, counterIJ, counterIJ2;
4681 
4682  double currTemp = temperature();
4683  double currPres = pressure();
4684 
4685 #ifdef DEBUG_MODE
4686  if (m_debugCalc) {
4687  printf("\n Debugging information from "
4688  "s_Pitzer_dlnMolalityActCoeff_dP()\n");
4689  }
4690 #endif
4691  /*
4692  * Make sure the counter variables are setup
4693  */
4694  counterIJ_setup();
4695 
4696  /*
4697  * ---------- Calculate common sums over solutes ---------------------
4698  */
4699  for (n = 1; n < m_kk; n++) {
4700  // ionic strength
4701  Is += charge(n) * charge(n) * molality[n];
4702  // total molar charge
4703  molarcharge += fabs(charge(n)) * molality[n];
4704  molalitysum += molality[n];
4705  }
4706  Is *= 0.5;
4707 
4708  /*
4709  * Store the ionic molality in the object for reference.
4710  */
4711  m_IionicMolality = Is;
4712  sqrtIs = sqrt(Is);
4713 #ifdef DEBUG_MODE
4714  if (m_debugCalc) {
4715  printf(" Step 1: \n");
4716  printf(" ionic strenth = %14.7le \n total molar "
4717  "charge = %14.7le \n", Is, molarcharge);
4718  }
4719 #endif
4720 
4721  /*
4722  * The following call to calc_lambdas() calculates all 16 elements
4723  * of the elambda and elambda1 arrays, given the value of the
4724  * ionic strength (Is)
4725  */
4726  calc_lambdas(Is);
4727 
4728 
4729  /*
4730  * ----- Step 2: Find the coefficients E-theta and -------------------
4731  * E-thetaprime for all combinations of positive
4732  * unlike charges up to 4
4733  */
4734 #ifdef DEBUG_MODE
4735  if (m_debugCalc) {
4736  printf(" Step 2: \n");
4737  }
4738 #endif
4739  for (z1 = 1; z1 <=4; z1++) {
4740  for (z2 =1; z2 <=4; z2++) {
4741  calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
4742 #ifdef DEBUG_MODE
4743  if (m_debugCalc) {
4744  printf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n",
4745  z1, z2, etheta[z1][z2], etheta_prime[z1][z2]);
4746  }
4747 #endif
4748  }
4749  }
4750 
4751 #ifdef DEBUG_MODE
4752  if (m_debugCalc) {
4753  printf(" Step 3: \n");
4754  printf(" Species Species g(x) "
4755  " hfunc(x) \n");
4756  }
4757 #endif
4758 
4759  /*
4760  *
4761  * calculate g(x) and hfunc(x) for each cation-anion pair MX
4762  * In the original literature, hfunc, was called gprime. However,
4763  * it's not the derivative of g(x), so I renamed it.
4764  */
4765  for (i = 1; i < (m_kk - 1); i++) {
4766  for (j = (i+1); j < m_kk; j++) {
4767  /*
4768  * Find the counterIJ for the symmetric binary interaction
4769  */
4770  n = m_kk*i + j;
4771  counterIJ = m_CounterIJ[n];
4772  /*
4773  * Only loop over oppositely charge species
4774  */
4775  if (charge(i)*charge(j) < 0) {
4776  /*
4777  * x is a reduced function variable
4778  */
4779  x1 = sqrtIs * alpha1MX[counterIJ];
4780  if (x1 > 1.0E-100) {
4781  gfunc[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 * x1);
4782  hfunc[counterIJ] = -2.0*
4783  (1.0-(1.0 + x1 + 0.5 * x1 * x1) * exp(-x1)) / (x1 * x1);
4784  } else {
4785  gfunc[counterIJ] = 0.0;
4786  hfunc[counterIJ] = 0.0;
4787  }
4788 
4789  if (beta2MX_P[counterIJ] != 0.0) {
4790  x2 = sqrtIs * alpha2MX[counterIJ];
4791  if (x2 > 1.0E-100) {
4792  g2func[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
4793  h2func[counterIJ] = -2.0 *
4794  (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
4795  } else {
4796  g2func[counterIJ] = 0.0;
4797  h2func[counterIJ] = 0.0;
4798  }
4799  }
4800  } else {
4801  gfunc[counterIJ] = 0.0;
4802  hfunc[counterIJ] = 0.0;
4803  }
4804 #ifdef DEBUG_MODE
4805  if (m_debugCalc) {
4806  sni = speciesName(i);
4807  snj = speciesName(j);
4808  printf(" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(),
4809  gfunc[counterIJ], hfunc[counterIJ]);
4810  }
4811 #endif
4812  }
4813  }
4814 
4815  /*
4816  * ------- SUBSECTION TO CALCULATE BMX_P, BprimeMX_P, BphiMX_P ----------
4817  * ------- These are now temperature derivatives of the
4818  * previously calculated quantities.
4819  */
4820 #ifdef DEBUG_MODE
4821  if (m_debugCalc) {
4822  printf(" Step 4: \n");
4823  printf(" Species Species BMX "
4824  "BprimeMX BphiMX \n");
4825  }
4826 #endif
4827 
4828  for (i = 1; i < m_kk - 1; i++) {
4829  for (j = i+1; j < m_kk; j++) {
4830  /*
4831  * Find the counterIJ for the symmetric binary interaction
4832  */
4833  n = m_kk*i + j;
4834  counterIJ = m_CounterIJ[n];
4835  /*
4836  * both species have a non-zero charge, and one is positive
4837  * and the other is negative
4838  */
4839  if (charge(i)*charge(j) < 0.0) {
4840  BMX_P[counterIJ] = beta0MX_P[counterIJ]
4841  + beta1MX_P[counterIJ] * gfunc[counterIJ]
4842  + beta2MX_P[counterIJ] * g2func[counterIJ];
4843 #ifdef DEBUG_MODE
4844  if (m_debugCalc) {
4845  printf("%d %g: %g %g %g %g\n",
4846  (int) counterIJ, BMX_P[counterIJ], beta0MX_P[counterIJ],
4847  beta1MX_P[counterIJ], beta2MX_P[counterIJ], gfunc[counterIJ]);
4848  }
4849 #endif
4850  if (Is > 1.0E-150) {
4851  BprimeMX_P[counterIJ] = (beta1MX_P[counterIJ] * hfunc[counterIJ]/Is +
4852  beta2MX_P[counterIJ] * h2func[counterIJ]/Is);
4853  } else {
4854  BprimeMX_P[counterIJ] = 0.0;
4855  }
4856  BphiMX_P[counterIJ] = BMX_P[counterIJ] + Is*BprimeMX_P[counterIJ];
4857  } else {
4858  BMX_P[counterIJ] = 0.0;
4859  BprimeMX_P[counterIJ] = 0.0;
4860  BphiMX_P[counterIJ] = 0.0;
4861  }
4862 #ifdef DEBUG_MODE
4863  if (m_debugCalc) {
4864  sni = speciesName(i);
4865  snj = speciesName(j);
4866  printf(" %-16s %-16s %11.7f %11.7f %11.7f \n",
4867  sni.c_str(), snj.c_str(),
4868  BMX_P[counterIJ], BprimeMX_P[counterIJ], BphiMX_P[counterIJ]);
4869  }
4870 #endif
4871  }
4872  }
4873 
4874  /*
4875  * --------- SUBSECTION TO CALCULATE CMX_P ----------
4876  */
4877 #ifdef DEBUG_MODE
4878  if (m_debugCalc) {
4879  printf(" Step 5: \n");
4880  printf(" Species Species CMX \n");
4881  }
4882 #endif
4883  for (i = 1; i < m_kk-1; i++) {
4884  for (j = i+1; j < m_kk; j++) {
4885  /*
4886  * Find the counterIJ for the symmetric binary interaction
4887  */
4888  n = m_kk*i + j;
4889  counterIJ = m_CounterIJ[n];
4890  /*
4891  * both species have a non-zero charge, and one is positive
4892  * and the other is negative
4893  */
4894  if (charge(i)*charge(j) < 0.0) {
4895  CMX_P[counterIJ] = CphiMX_P[counterIJ]/
4896  (2.0* sqrt(fabs(charge(i)*charge(j))));
4897  } else {
4898  CMX_P[counterIJ] = 0.0;
4899  }
4900 #ifdef DEBUG_MODE
4901  if (m_debugCalc) {
4902  sni = speciesName(i);
4903  snj = speciesName(j);
4904  printf(" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(),
4905  CMX_P[counterIJ]);
4906  }
4907 #endif
4908  }
4909  }
4910 
4911  /*
4912  * ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ----------
4913  * --------
4914  */
4915 #ifdef DEBUG_MODE
4916  if (m_debugCalc) {
4917  printf(" Step 6: \n");
4918  printf(" Species Species Phi_ij "
4919  " Phiprime_ij Phi^phi_ij \n");
4920  }
4921 #endif
4922  for (i = 1; i < m_kk-1; i++) {
4923  for (j = i+1; j < m_kk; j++) {
4924  /*
4925  * Find the counterIJ for the symmetric binary interaction
4926  */
4927  n = m_kk*i + j;
4928  counterIJ = m_CounterIJ[n];
4929  /*
4930  * both species have a non-zero charge, and one is positive
4931  * and the other is negative
4932  */
4933  if (charge(i)*charge(j) > 0) {
4934  z1 = (int) fabs(charge(i));
4935  z2 = (int) fabs(charge(j));
4936  //Phi[counterIJ] = thetaij_L[counterIJ] + etheta[z1][z2];
4937  Phi_P[counterIJ] = thetaij_P[counterIJ];
4938  //Phiprime[counterIJ] = etheta_prime[z1][z2];
4939  Phiprime[counterIJ] = 0.0;
4940  Phiphi_P[counterIJ] = Phi_P[counterIJ] + Is * Phiprime[counterIJ];
4941  } else {
4942  Phi_P[counterIJ] = 0.0;
4943  Phiprime[counterIJ] = 0.0;
4944  Phiphi_P[counterIJ] = 0.0;
4945  }
4946 #ifdef DEBUG_MODE
4947  if (m_debugCalc) {
4948  sni = speciesName(i);
4949  snj = speciesName(j);
4950  printf(" %-16s %-16s %10.6f %10.6f %10.6f \n",
4951  sni.c_str(), snj.c_str(),
4952  Phi_P[counterIJ], Phiprime[counterIJ], Phiphi_P[counterIJ]);
4953  }
4954 #endif
4955  }
4956  }
4957 
4958  /*
4959  * ----------- SUBSECTION FOR CALCULATION OF dFdT ---------------------
4960  */
4961 #ifdef DEBUG_MODE
4962  if (m_debugCalc) {
4963  printf(" Step 7: \n");
4964  }
4965 #endif
4966  // A_Debye_Huckel = 0.5092; (units = sqrt(kg/gmol))
4967  // A_Debye_Huckel = 0.5107; <- This value is used to match GWB data
4968  // ( A * ln(10) = 1.17593)
4969  // Aphi = A_Debye_Huckel * 2.30258509 / 3.0;
4970 
4971  double dA_DebyedP = dA_DebyedP_TP(currTemp, currPres);
4972  double dAphidP = dA_DebyedP /3.0;
4973 #ifdef DEBUG_MODE
4974  //dAphidT = 0.0;
4975 #endif
4976  //F = -Aphi * ( sqrt(Is) / (1.0 + 1.2*sqrt(Is))
4977  // + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
4978  //dAphidT = Al / (4.0 * GasConstant * T * T);
4979  dFdP = -dAphidP * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
4980  + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
4981 #ifdef DEBUG_MODE
4982  if (m_debugCalc) {
4983  printf(" initial value of dFdP = %10.6f \n", dFdP);
4984  }
4985 #endif
4986  for (i = 1; i < m_kk-1; i++) {
4987  for (j = i+1; j < m_kk; j++) {
4988  /*
4989  * Find the counterIJ for the symmetric binary interaction
4990  */
4991  n = m_kk*i + j;
4992  counterIJ = m_CounterIJ[n];
4993  /*
4994  * both species have a non-zero charge, and one is positive
4995  * and the other is negative
4996  */
4997  if (charge(i)*charge(j) < 0) {
4998  dFdP = dFdP + molality[i]*molality[j] * BprimeMX_P[counterIJ];
4999  }
5000  /*
5001  * Both species have a non-zero charge, and they
5002  * have the same sign, e.g., both positive or both negative.
5003  */
5004  if (charge(i)*charge(j) > 0) {
5005  dFdP = dFdP + molality[i]*molality[j] * Phiprime[counterIJ];
5006  }
5007 #ifdef DEBUG_MODE
5008  if (m_debugCalc) {
5009  printf(" dFdP = %10.6f \n", dFdP);
5010  }
5011 #endif
5012  }
5013  }
5014 #ifdef DEBUG_MODE
5015  if (m_debugCalc) {
5016  printf(" Step 8: \n");
5017  }
5018 #endif
5019 
5020 
5021  for (i = 1; i < m_kk; i++) {
5022 
5023  /*
5024  * -------- SUBSECTION FOR CALCULATING THE dACTCOEFFdP FOR CATIONS -----
5025  */
5026  if (charge(i) > 0) {
5027  // species i is the cation (positive) to calc the actcoeff
5028  zsqdFdP = charge(i)*charge(i)*dFdP;
5029  sum1 = 0.0;
5030  sum2 = 0.0;
5031  sum3 = 0.0;
5032  sum4 = 0.0;
5033  sum5 = 0.0;
5034  for (j = 1; j < m_kk; j++) {
5035  /*
5036  * Find the counterIJ for the symmetric binary interaction
5037  */
5038  n = m_kk*i + j;
5039  counterIJ = m_CounterIJ[n];
5040 
5041  if (charge(j) < 0.0) {
5042  // sum over all anions
5043  sum1 = sum1 + molality[j]*
5044  (2.0*BMX_P[counterIJ] + molarcharge*CMX_P[counterIJ]);
5045  if (j < m_kk-1) {
5046  /*
5047  * This term is the ternary interaction involving the
5048  * non-duplicate sum over double anions, j, k, with
5049  * respect to the cation, i.
5050  */
5051  for (size_t k = j+1; k < m_kk; k++) {
5052  // an inner sum over all anions
5053  if (charge(k) < 0.0) {
5054  n = k + j * m_kk + i * m_kk * m_kk;
5055  sum3 = sum3 + molality[j]*molality[k]*psi_ijk_P[n];
5056  }
5057  }
5058  }
5059  }
5060 
5061 
5062  if (charge(j) > 0.0) {
5063  // sum over all cations
5064  if (j != i) {
5065  sum2 = sum2 + molality[j]*(2.0*Phi_P[counterIJ]);
5066  }
5067  for (size_t k = 1; k < m_kk; k++) {
5068  if (charge(k) < 0.0) {
5069  // two inner sums over anions
5070 
5071  n = k + j * m_kk + i * m_kk * m_kk;
5072  sum2 = sum2 + molality[j]*molality[k]*psi_ijk_P[n];
5073  /*
5074  * Find the counterIJ for the j,k interaction
5075  */
5076  n = m_kk*j + k;
5077  counterIJ2 = m_CounterIJ[n];
5078  sum4 = sum4 + (fabs(charge(i))*
5079  molality[j]*molality[k]*CMX_P[counterIJ2]);
5080  }
5081  }
5082  }
5083 
5084  /*
5085  * for Anions, do the neutral species interaction
5086  */
5087  if (charge(j) == 0) {
5088  sum5 = sum5 + molality[j]*2.0*m_Lambda_nj_P(j,i);
5089  /*
5090  * Zeta interaction term
5091  */
5092  for (size_t k = 1; k < m_kk; k++) {
5093  if (charge(k) < 0.0) {
5094  size_t izeta = j;
5095  size_t jzeta = i;
5096  n = izeta * m_kk * m_kk + jzeta * m_kk + k;
5097  double zeta_P = psi_ijk_P[n];
5098  if (zeta_P != 0.0) {
5099  sum5 = sum5 + molality[j]*molality[k]*zeta_P;
5100  }
5101  }
5102  }
5103  }
5104  }
5105 
5106  /*
5107  * Add all of the contributions up to yield the log of the
5108  * solute activity coefficients (molality scale)
5109  */
5111  zsqdFdP + sum1 + sum2 + sum3 + sum4 + sum5;
5112 
5113 #ifdef DEBUG_MODE
5114  if (m_debugCalc) {
5115  sni = speciesName(i);
5116  printf(" %-16s lngamma[i]=%10.6f \n",
5117  sni.c_str(), m_dlnActCoeffMolaldP_Unscaled[i]);
5118  printf(" %12g %12g %12g %12g %12g %12g\n",
5119  zsqdFdP, sum1, sum2, sum3, sum4, sum5);
5120  }
5121 #endif
5122  }
5123 
5124 
5125  /*
5126  * ------ SUBSECTION FOR CALCULATING THE dACTCOEFFdP FOR ANIONS ------
5127  */
5128  if (charge(i) < 0) {
5129  // species i is an anion (negative)
5130  zsqdFdP = charge(i)*charge(i)*dFdP;
5131  sum1 = 0.0;
5132  sum2 = 0.0;
5133  sum3 = 0.0;
5134  sum4 = 0.0;
5135  sum5 = 0.0;
5136  for (j = 1; j < m_kk; j++) {
5137  /*
5138  * Find the counterIJ for the symmetric binary interaction
5139  */
5140  n = m_kk*i + j;
5141  counterIJ = m_CounterIJ[n];
5142 
5143  /*
5144  * For Anions, do the cation interactions.
5145  */
5146  if (charge(j) > 0) {
5147  sum1 = sum1 + molality[j]*
5148  (2.0*BMX_P[counterIJ] + molarcharge*CMX_P[counterIJ]);
5149  if (j < m_kk-1) {
5150  for (size_t k = j+1; k < m_kk; k++) {
5151  // an inner sum over all cations
5152  if (charge(k) > 0) {
5153  n = k + j * m_kk + i * m_kk * m_kk;
5154  sum3 = sum3 + molality[j]*molality[k]*psi_ijk_P[n];
5155  }
5156  }
5157  }
5158  }
5159 
5160  /*
5161  * For Anions, do the other anion interactions.
5162  */
5163  if (charge(j) < 0.0) {
5164  // sum over all anions
5165  if (j != i) {
5166  sum2 = sum2 + molality[j]*(2.0*Phi_P[counterIJ]);
5167  }
5168  for (size_t k = 1; k < m_kk; k++) {
5169  if (charge(k) > 0.0) {
5170  // two inner sums over cations
5171  n = k + j * m_kk + i * m_kk * m_kk;
5172  sum2 = sum2 + molality[j]*molality[k]*psi_ijk_P[n];
5173  /*
5174  * Find the counterIJ for the symmetric binary interaction
5175  */
5176  n = m_kk*j + k;
5177  counterIJ2 = m_CounterIJ[n];
5178  sum4 = sum4 +
5179  (fabs(charge(i))*
5180  molality[j]*molality[k]*CMX_P[counterIJ2]);
5181  }
5182  }
5183  }
5184 
5185  /*
5186  * for Anions, do the neutral species interaction
5187  */
5188  if (charge(j) == 0.0) {
5189  sum5 = sum5 + molality[j]*2.0*m_Lambda_nj_P(j,i);
5190  /*
5191  * Zeta interaction term
5192  */
5193  for (size_t k = 1; k < m_kk; k++) {
5194  if (charge(k) > 0.0) {
5195  size_t izeta = j;
5196  size_t jzeta = k;
5197  size_t kzeta = i;
5198  n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
5199  double zeta_P = psi_ijk_P[n];
5200  if (zeta_P != 0.0) {
5201  sum5 = sum5 + molality[j]*molality[k]*zeta_P;
5202  }
5203  }
5204  }
5205  }
5206  }
5208  zsqdFdP + sum1 + sum2 + sum3 + sum4 + sum5;
5209 #ifdef DEBUG_MODE
5210  if (m_debugCalc) {
5211  sni = speciesName(i);
5212  printf(" %-16s lndactcoeffmolaldP[i]=%10.6f \n",
5213  sni.c_str(), m_dlnActCoeffMolaldP_Unscaled[i]);
5214  printf(" %12g %12g %12g %12g %12g %12g\n",
5215  zsqdFdP, sum1, sum2, sum3, sum4, sum5);
5216  }
5217 #endif
5218  }
5219 
5220  /*
5221  * ------ SUBSECTION FOR CALCULATING d NEUTRAL SOLUTE ACT COEFF dP -------
5222  */
5223  if (charge(i) == 0.0) {
5224  sum1 = 0.0;
5225  sum3 = 0.0;
5226  for (j = 1; j < m_kk; j++) {
5227  sum1 += molality[j]*2.0*m_Lambda_nj_P(i,j);
5228  /*
5229  * Zeta term -> we piggyback on the psi term
5230  */
5231  if (charge(j) > 0.0) {
5232  for (size_t k = 1; k < m_kk; k++) {
5233  if (charge(k) < 0.0) {
5234  n = k + j * m_kk + i * m_kk * m_kk;
5235  sum3 = sum3 + molality[j]*molality[k]*psi_ijk_P[n];
5236  }
5237  }
5238  }
5239  }
5240  sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_P[i];
5241  m_dlnActCoeffMolaldP_Unscaled[i] = sum1 + sum2 + sum3;
5242 #ifdef DEBUG_MODE
5243  if (m_debugCalc) {
5244  sni = speciesName(i);
5245  printf(" %-16s dlnActCoeffMolaldP[i]=%10.6f \n",
5246  sni.c_str(), m_dlnActCoeffMolaldP_Unscaled[i]);
5247  }
5248 #endif
5249  }
5250 
5251  }
5252 #ifdef DEBUG_MODE
5253  if (m_debugCalc) {
5254  printf(" Step 9: \n");
5255  }
5256 #endif
5257 
5258  /*
5259  * ------ SUBSECTION FOR CALCULATING THE d OSMOTIC COEFF dP ---------
5260  */
5261  sum1 = 0.0;
5262  sum2 = 0.0;
5263  sum3 = 0.0;
5264  sum4 = 0.0;
5265  sum5 = 0.0;
5266  double sum6 = 0.0;
5267  double sum7 = 0.0;
5268  /*
5269  * term1 is the temperature derivative of the
5270  * DH term in the osmotic coefficient expression
5271  * b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer
5272  * implementations.
5273  * Is = Ionic strength on the molality scale (units of (gmol/kg))
5274  * Aphi = A_Debye / 3 (units of sqrt(kg/gmol))
5275  */
5276  term1 = -dAphidP * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
5277 
5278  for (j = 1; j < m_kk; j++) {
5279  /*
5280  * Loop Over Cations
5281  */
5282  if (charge(j) > 0.0) {
5283  for (size_t k = 1; k < m_kk; k++) {
5284  if (charge(k) < 0.0) {
5285  /*
5286  * Find the counterIJ for the symmetric j,k binary interaction
5287  */
5288  n = m_kk*j + k;
5289  counterIJ = m_CounterIJ[n];
5290 
5291  sum1 = sum1 + molality[j]*molality[k]*
5292  (BphiMX_P[counterIJ] + molarcharge*CMX_P[counterIJ]);
5293  }
5294  }
5295 
5296  for (size_t k = j+1; k < m_kk; k++) {
5297  if (j == (m_kk-1)) {
5298  // we should never reach this step
5299  printf("logic error 1 in Step 9 of hmw_act");
5300  exit(EXIT_FAILURE);
5301  }
5302  if (charge(k) > 0.0) {
5303  /*
5304  * Find the counterIJ for the symmetric j,k binary interaction
5305  * between 2 cations.
5306  */
5307  n = m_kk*j + k;
5308  counterIJ = m_CounterIJ[n];
5309  sum2 = sum2 + molality[j]*molality[k]*Phiphi_P[counterIJ];
5310  for (m = 1; m < m_kk; m++) {
5311  if (charge(m) < 0.0) {
5312  // species m is an anion
5313  n = m + k * m_kk + j * m_kk * m_kk;
5314  sum2 = sum2 +
5315  molality[j]*molality[k]*molality[m]*psi_ijk_P[n];
5316  }
5317  }
5318  }
5319  }
5320  }
5321 
5322 
5323  /*
5324  * Loop Over Anions
5325  */
5326  if (charge(j) < 0) {
5327  for (size_t k = j+1; k < m_kk; k++) {
5328  if (j == m_kk-1) {
5329  // we should never reach this step
5330  printf("logic error 2 in Step 9 of hmw_act");
5331  exit(EXIT_FAILURE);
5332  }
5333  if (charge(k) < 0) {
5334  /*
5335  * Find the counterIJ for the symmetric j,k binary interaction
5336  * between two anions
5337  */
5338  n = m_kk*j + k;
5339  counterIJ = m_CounterIJ[n];
5340 
5341  sum3 = sum3 + molality[j]*molality[k]*Phiphi_P[counterIJ];
5342  for (m = 1; m < m_kk; m++) {
5343  if (charge(m) > 0.0) {
5344  n = m + k * m_kk + j * m_kk * m_kk;
5345  sum3 = sum3 +
5346  molality[j]*molality[k]*molality[m]*psi_ijk_P[n];
5347  }
5348  }
5349  }
5350  }
5351  }
5352 
5353  /*
5354  * Loop Over Neutral Species
5355  */
5356  if (charge(j) == 0) {
5357  for (size_t k = 1; k < m_kk; k++) {
5358  if (charge(k) < 0.0) {
5359  sum4 = sum4 + molality[j]*molality[k]*m_Lambda_nj_P(j,k);
5360  }
5361  if (charge(k) > 0.0) {
5362  sum5 = sum5 + molality[j]*molality[k]*m_Lambda_nj_P(j,k);
5363  }
5364  if (charge(k) == 0.0) {
5365  if (k > j) {
5366  sum6 = sum6 + molality[j]*molality[k]*m_Lambda_nj_P(j,k);
5367  } else if (k == j) {
5368  sum6 = sum6 + 0.5 * molality[j]*molality[k]*m_Lambda_nj_P(j,k);
5369  }
5370  }
5371  if (charge(k) < 0.0) {
5372  size_t izeta = j;
5373  for (m = 1; m < m_kk; m++) {
5374  if (charge(m) > 0.0) {
5375  size_t jzeta = m;
5376  n = k + jzeta * m_kk + izeta * m_kk * m_kk;
5377  double zeta_P = psi_ijk_P[n];
5378  if (zeta_P != 0.0) {
5379  sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_P;
5380  }
5381  }
5382  }
5383  }
5384  }
5385 
5386  sum7 += molality[j] * molality[j] * molality[j] * m_Mu_nnn_P[j];
5387  }
5388  }
5389  sum_m_phi_minus_1 = 2.0 *
5390  (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
5391 
5392  /*
5393  * Calculate the osmotic coefficient from
5394  * osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i)
5395  */
5396  if (molalitysum > 1.0E-150) {
5397  d_osmotic_coef_dP = 0.0 + (sum_m_phi_minus_1 / molalitysum);
5398  } else {
5399  d_osmotic_coef_dP = 0.0;
5400  }
5401 #ifdef DEBUG_MODE
5402  if (m_debugCalc) {
5403  printf(" term1=%10.6f sum1=%10.6f sum2=%10.6f "
5404  "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n",
5405  term1, sum1, sum2, sum3, sum4, sum5);
5406  printf(" sum_m_phi_minus_1=%10.6f d_osmotic_coef_dP =%10.6f\n",
5407  sum_m_phi_minus_1, d_osmotic_coef_dP);
5408  }
5409 
5410  if (m_debugCalc) {
5411  printf(" Step 10: \n");
5412  }
5413 #endif
5414  d_lnwateract_dP = -(m_weightSolvent/1000.0) * molalitysum * d_osmotic_coef_dP;
5415 
5416 
5417  /*
5418  * In Cantera, we define the activity coefficient of the solvent as
5419  *
5420  * act_0 = actcoeff_0 * Xmol_0
5421  *
5422  * We have just computed act_0. However, this routine returns
5423  * ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
5424  */
5425  //double xmolSolvent = moleFraction(m_indexSolvent);
5426  m_dlnActCoeffMolaldP_Unscaled[0] = d_lnwateract_dP;
5427 #ifdef DEBUG_MODE
5428  if (m_debugCalc) {
5429  double d_wateract_dP = exp(d_lnwateract_dP);
5430  printf(" d_ln_a_water_dP = %10.6f d_a_water_dP=%10.6f\n\n",
5431  d_lnwateract_dP, d_wateract_dP);
5432  }
5433 #endif
5434 
5435 }
5436 
5437 void HMWSoln::calc_lambdas(double is) const
5438 {
5439  double aphi, dj, jfunc, jprime, t, x, zprod;
5440  int i, ij, j;
5441  /*
5442  * Coefficients c1-c4 are used to approximate
5443  * the integral function "J";
5444  * aphi is the Debye-Huckel constant at 25 C
5445  */
5446 
5447  double c1 = 4.581, c2 = 0.7237, c3 = 0.0120, c4 = 0.528;
5448 
5449  aphi = 0.392; /* Value at 25 C */
5450 #ifdef DEBUG_MODE
5451  if (m_debugCalc) {
5452  printf(" Is = %g\n", is);
5453  }
5454 #endif
5455  if (is < 1.0E-150) {
5456  for (i = 0; i < 17; i++) {
5457  elambda[i] = 0.0;
5458  elambda1[i] = 0.0;
5459  }
5460  return;
5461  }
5462  /*
5463  * Calculate E-lambda terms for charge combinations of like sign,
5464  * using method of Pitzer (1975). Charges up to 4 are calculated.
5465  */
5466 
5467  for (i=1; i<=4; i++) {
5468  for (j=i; j<=4; j++) {
5469  ij = i*j;
5470  /*
5471  * calculate the product of the charges
5472  */
5473  zprod = (double)ij;
5474  /*
5475  * calculate Xmn (A1) from Harvie, Weare (1980).
5476  */
5477  x = 6.0* zprod * aphi * sqrt(is); /* eqn 23 */
5478 
5479  jfunc = x / (4.0 + c1*pow(x,-c2)*exp(-c3*pow(x,c4))); /* eqn 47 */
5480 
5481  t = c3 * c4 * pow(x,c4);
5482  dj = c1* pow(x,(-c2-1.0)) * (c2+t) * exp(-c3*pow(x,c4));
5483  jprime = (jfunc/x)*(1.0 + jfunc*dj);
5484 
5485  elambda[ij] = zprod*jfunc / (4.0*is); /* eqn 14 */
5486  elambda1[ij] = (3.0*zprod*zprod*aphi*jprime/(4.0*sqrt(is))
5487  - elambda[ij])/is;
5488 #ifdef DEBUG_MODE
5489  if (m_debugCalc) {
5490  printf(" ij = %d, elambda = %g, elambda1 = %g\n",
5491  ij, elambda[ij], elambda1[ij]);
5492  }
5493 #endif
5494  }
5495  }
5496 }
5497 
5498 void HMWSoln::calc_thetas(int z1, int z2,
5499  double* etheta, double* etheta_prime) const
5500 {
5501  int i, j;
5502  double f1, f2;
5503 
5504  /*
5505  * Calculate E-theta(i) and E-theta'(I) using method of
5506  * Pitzer (1987)
5507  */
5508  i = abs(z1);
5509  j = abs(z2);
5510 
5511 #ifdef DEBUG_MODE
5512  if (i > 4 || j > 4) {
5513  printf("we shouldn't be here\n");
5514  exit(EXIT_FAILURE);
5515  }
5516 #endif
5517 
5518  if ((i == 0) || (j == 0)) {
5519  printf("ERROR calc_thetas called with one species being neutral\n");
5520  exit(EXIT_FAILURE);
5521  }
5522 
5523  /*
5524  * Check to see if the charges are of opposite sign. If they are of
5525  * opposite sign then their etheta interaction is zero.
5526  */
5527  if (z1*z2 < 0) {
5528  *etheta = 0.0;
5529  *etheta_prime = 0.0;
5530  }
5531  /*
5532  * Actually calculate the interaction.
5533  */
5534  else {
5535  f1 = (double)i / (2.0 * j);
5536  f2 = (double)j / (2.0 * i);
5537  *etheta = elambda[i*j] - f1*elambda[j*j] - f2*elambda[i*i];
5538  *etheta_prime = elambda1[i*j] - f1*elambda1[j*j] - f2*elambda1[i*i];
5539  }
5540 }
5541 
5543 {
5544  double tmp;
5545  /*
5546  * Calculate the molalities. Currently, the molalities
5547  * may not be current with respect to the contents of the
5548  * State objects' data.
5549  */
5550  calcMolalities();
5551 
5552  double xmolSolvent = moleFraction(m_indexSolvent);
5553  double xx = std::max(m_xmolSolventMIN, xmolSolvent);
5554  if (IMS_typeCutoff_ == 0) {
5555  for (size_t k = 1; k < m_kk; k++) {
5556  IMS_lnActCoeffMolal_[k]= 0.0;
5557  }
5558  IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
5559  return;
5560  } else if (IMS_typeCutoff_ == 1) {
5561  if (xmolSolvent > 3.0 * IMS_X_o_cutoff_/2.0) {
5562  for (size_t k = 1; k < m_kk; k++) {
5563  IMS_lnActCoeffMolal_[k]= 0.0;
5564  }
5565  IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
5566  return;
5567  } else if (xmolSolvent < IMS_X_o_cutoff_/2.0) {
5568  tmp = log(xx * IMS_gamma_k_min_);
5569  for (size_t k = 1; k < m_kk; k++) {
5570  IMS_lnActCoeffMolal_[k]= tmp;
5571  }
5573  return;
5574  } else {
5575  /*
5576  * If we are in the middle region, calculate the connecting polynomials
5577  */
5578  double xminus = xmolSolvent - IMS_X_o_cutoff_/2.0;
5579  double xminus2 = xminus * xminus;
5580  double xminus3 = xminus2 * xminus;
5581  double x_o_cut2 = IMS_X_o_cutoff_ * IMS_X_o_cutoff_;
5582  double x_o_cut3 = x_o_cut2 * IMS_X_o_cutoff_;
5583 
5584  double h2 = 3.5 * xminus2 / IMS_X_o_cutoff_ - 2.0 * xminus3 / x_o_cut2;
5585  double h2_prime = 7.0 * xminus / IMS_X_o_cutoff_ - 6.0 * xminus2 / x_o_cut2;
5586 
5587  double h1 = (1.0 - 3.0 * xminus2 / x_o_cut2 + 2.0 * xminus3/ x_o_cut3);
5588  double h1_prime = (- 6.0 * xminus / x_o_cut2 + 6.0 * xminus2/ x_o_cut3);
5589 
5590  double h1_g = h1 / IMS_gamma_o_min_;
5591  double h1_g_prime = h1_prime / IMS_gamma_o_min_;
5592 
5593  double alpha = 1.0 / (exp(1.0) * IMS_gamma_k_min_);
5594  double h1_f = h1 * alpha;
5595  double h1_f_prime = h1_prime * alpha;
5596 
5597  double f = h2 + h1_f;
5598  double f_prime = h2_prime + h1_f_prime;
5599 
5600  double g = h2 + h1_g;
5601  double g_prime = h2_prime + h1_g_prime;
5602 
5603  tmp = (xmolSolvent/ g * g_prime + (1.0-xmolSolvent) / f * f_prime);
5604  double lngammak = -1.0 - log(f) + tmp * xmolSolvent;
5605  double lngammao =-log(g) - tmp * (1.0-xmolSolvent);
5606 
5607  tmp = log(xmolSolvent) + lngammak;
5608  for (size_t k = 1; k < m_kk; k++) {
5609  IMS_lnActCoeffMolal_[k]= tmp;
5610  }
5612  }
5613  }
5614  // Exponentials - trial 2
5615  else if (IMS_typeCutoff_ == 2) {
5616  if (xmolSolvent > IMS_X_o_cutoff_) {
5617  for (size_t k = 1; k < m_kk; k++) {
5618  IMS_lnActCoeffMolal_[k]= 0.0;
5619  }
5620  IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
5621  return;
5622  } else {
5623 
5624  double xoverc = xmolSolvent/IMS_cCut_;
5625  double eterm = std::exp(-xoverc);
5626 
5627  double fptmp = IMS_bfCut_ - IMS_afCut_ / IMS_cCut_ - IMS_bfCut_*xoverc
5628  + 2.0*IMS_dfCut_*xmolSolvent - IMS_dfCut_*xmolSolvent*xoverc;
5629  double f_prime = 1.0 + eterm*fptmp;
5630  double f = xmolSolvent + IMS_efCut_
5631  + eterm * (IMS_afCut_ + xmolSolvent * (IMS_bfCut_ + IMS_dfCut_*xmolSolvent));
5632 
5633  double gptmp = IMS_bgCut_ - IMS_agCut_ / IMS_cCut_ - IMS_bgCut_*xoverc
5634  + 2.0*IMS_dgCut_*xmolSolvent - IMS_dgCut_*xmolSolvent*xoverc;
5635  double g_prime = 1.0 + eterm*gptmp;
5636  double g = xmolSolvent + IMS_egCut_
5637  + eterm * (IMS_agCut_ + xmolSolvent * (IMS_bgCut_ + IMS_dgCut_*xmolSolvent));
5638 
5639  tmp = (xmolSolvent / g * g_prime + (1.0 - xmolSolvent) / f * f_prime);
5640  double lngammak = -1.0 - log(f) + tmp * xmolSolvent;
5641  double lngammao =-log(g) - tmp * (1.0-xmolSolvent);
5642 
5643  tmp = log(xx) + lngammak;
5644  for (size_t k = 1; k < m_kk; k++) {
5645  IMS_lnActCoeffMolal_[k]= tmp;
5646  }
5648  }
5649  }
5650  return;
5651 }
5652 
5654 {
5655  size_t i, j, k;
5656  std::string sni, snj;
5657  calcMolalities();
5658  double* molality = DATA_PTR(m_molalitiesCropped);
5659  double* moleF = DATA_PTR(m_tmpV);
5660  /*
5661  * Update the coefficients wrt Temperature
5662  * Calculate the derivatives as well
5663  */
5665  getMoleFractions(moleF);
5666 
5667  printf("Index Name MoleF MolalityCropped Charge\n");
5668  for (k = 0; k < m_kk; k++) {
5669  sni = speciesName(k);
5670  printf("%2s %-16s %14.7le %14.7le %5.1f \n",
5671  int2str(k).c_str(), sni.c_str(), moleF[k], molality[k], charge(k));
5672  }
5673 
5674  printf("\n Species Species beta0MX "
5675  "beta1MX beta2MX CphiMX alphaMX thetaij \n");
5676  for (i = 1; i < m_kk - 1; i++) {
5677  sni = speciesName(i);
5678  for (j = i+1; j < m_kk; j++) {
5679  snj = speciesName(j);
5680  size_t n = i * m_kk + j;
5681  size_t ct = m_CounterIJ[n];
5682  printf(" %-16s %-16s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f \n",
5683  sni.c_str(), snj.c_str(),
5684  m_Beta0MX_ij[ct], m_Beta1MX_ij[ct],
5685  m_Beta2MX_ij[ct], m_CphiMX_ij[ct],
5686  m_Alpha1MX_ij[ct], m_Theta_ij[ct]);
5687 
5688 
5689  }
5690  }
5691 
5692  printf("\n Species Species Species "
5693  "psi \n");
5694  for (i = 1; i < m_kk; i++) {
5695  sni = speciesName(i);
5696  for (j = 1; j < m_kk; j++) {
5697  snj = speciesName(j);
5698  for (k = 1; k < m_kk; k++) {
5699  std::string snk = speciesName(k);
5700  size_t n = k + j * m_kk + i * m_kk * m_kk;
5701  if (m_Psi_ijk[n] != 0.0) {
5702  printf(" %-16s %-16s %-16s %9.5f \n",
5703  sni.c_str(), snj.c_str(),
5704  snk.c_str(), m_Psi_ijk[n]);
5705  }
5706  }
5707  }
5708  }
5709 }
5710 
5711 void HMWSoln::applyphScale(doublereal* acMolality) const
5712 {
5714  return;
5715  }
5717  doublereal lnGammaClMs2 = s_NBS_CLM_lnMolalityActCoeff();
5718  doublereal lnGammaCLMs1 = m_lnActCoeffMolal_Unscaled[m_indexCLM];
5719  doublereal afac = -1.0 *(lnGammaClMs2 - lnGammaCLMs1);
5720  for (size_t k = 0; k < m_kk; k++) {
5721  acMolality[k] *= exp(charge(k) * afac);
5722  }
5723 }
5724 
5726 {
5729  return;
5730  }
5732  doublereal lnGammaClMs2 = s_NBS_CLM_lnMolalityActCoeff();
5733  doublereal lnGammaCLMs1 = m_lnActCoeffMolal_Unscaled[m_indexCLM];
5734  doublereal afac = -1.0 *(lnGammaClMs2 - lnGammaCLMs1);
5735  for (size_t k = 0; k < m_kk; k++) {
5737  }
5738 }
5739 
5741 {
5744  return;
5745  }
5747  doublereal dlnGammaClM_dT_s2 = s_NBS_CLM_dlnMolalityActCoeff_dT();
5748  doublereal dlnGammaCLM_dT_s1 = m_dlnActCoeffMolaldT_Unscaled[m_indexCLM];
5749  doublereal afac = -1.0 *(dlnGammaClM_dT_s2 - dlnGammaCLM_dT_s1);
5750  for (size_t k = 0; k < m_kk; k++) {
5752  }
5753 }
5754 
5756 {
5759  return;
5760  }
5762  doublereal d2lnGammaClM_dT2_s2 = s_NBS_CLM_d2lnMolalityActCoeff_dT2();
5763  doublereal d2lnGammaCLM_dT2_s1 = m_d2lnActCoeffMolaldT2_Unscaled[m_indexCLM];
5764  doublereal afac = -1.0 *(d2lnGammaClM_dT2_s2 - d2lnGammaCLM_dT2_s1);
5765  for (size_t k = 0; k < m_kk; k++) {
5767  }
5768 }
5769 
5771 {
5774  return;
5775  }
5777  doublereal dlnGammaClM_dP_s2 = s_NBS_CLM_dlnMolalityActCoeff_dP();
5778  doublereal dlnGammaCLM_dP_s1 = m_dlnActCoeffMolaldP_Unscaled[m_indexCLM];
5779  doublereal afac = -1.0 *(dlnGammaClM_dP_s2 - dlnGammaCLM_dP_s1);
5780  for (size_t k = 0; k < m_kk; k++) {
5782  }
5783 }
5784 
5786 {
5787  doublereal sqrtIs = sqrt(m_IionicMolality);
5788  doublereal A = m_A_Debye;
5789  doublereal lnGammaClMs2 = - A * sqrtIs /(1.0 + 1.5 * sqrtIs);
5790  return lnGammaClMs2;
5791 }
5792 
5794 {
5795  doublereal sqrtIs = sqrt(m_IionicMolality);
5796  doublereal dAdT = dA_DebyedT_TP();
5797  return - dAdT * sqrtIs /(1.0 + 1.5 * sqrtIs);
5798 }
5799 
5801 {
5802  doublereal sqrtIs = sqrt(m_IionicMolality);
5803  doublereal d2AdT2 = d2A_DebyedT2_TP();
5804  return - d2AdT2 * sqrtIs /(1.0 + 1.5 * sqrtIs);
5805 }
5806 
5808 {
5809  doublereal sqrtIs = sqrt(m_IionicMolality);
5810  doublereal dAdP = dA_DebyedP_TP();
5811  return - dAdP * sqrtIs /(1.0 + 1.5 * sqrtIs);
5812 }
5813 
5815 {
5816 #ifdef DEBUG_MODE
5817  return m_debugCalc;
5818 #else
5819  return 0;
5820 #endif
5821 }
5822 
5823 }
vector_fp m_Theta_ij_P
Derivative of Theta_ij[i][j] wrt P.
Definition: HMWSoln.h:2651
vector_fp m_BMX_IJ_L
Derivative of BMX_IJ wrt T.
Definition: HMWSoln.h:2923
vector_fp m_PhiPhi_IJ_P
Derivative of m_PhiPhi_IJ wrt P.
Definition: HMWSoln.h:3037
void s_updateIMS_lnMolalityActCoeff() const
This function will be called to update the internally stored natural logarithm of the molality activi...
Definition: HMWSoln.cpp:5542
vector_fp m_CphiMX_ij_LL
Derivative of Cphi_ij[i][j] wrt TT.
Definition: HMWSoln.h:2604
void s_update_dlnMolalityActCoeff_dP() const
This function calculates the pressure derivative of the natural logarithm of the molality activity co...
Definition: HMWSoln.cpp:4600
virtual doublereal pressure() const
Pressure.
Definition: HMWSoln.cpp:671
virtual doublereal density() const
Density (kg/m^3).
Definition: Phase.h:534
vector_fp m_Theta_ij_L
Derivative of Theta_ij[i][j] wrt T.
Definition: HMWSoln.h:2639
std::string int2str(const int n, const std::string &fmt)
Convert an int to a string using a format converter.
Definition: stringUtils.cpp:40
vector_fp m_speciesCharge_Stoich
Stoichiometric species charge -> This is for calculations of the ionic strength which ignore ion-ion ...
Definition: HMWSoln.h:2443
virtual void getParameters(int &n, doublereal *const c) const
Get the equation of state parameters in a vector.
Definition: HMWSoln.cpp:1014
WaterProps * m_waterProps
Pointer to the water property calculator.
Definition: HMWSoln.h:2423
vector_fp m_CMX_IJ_P
Derivative of m_CMX_IJ wrt P.
Definition: HMWSoln.h:3061
vector_fp m_gfunc_IJ
Various temporary arrays used in the calculation of the Pitzer activity coefficients.
Definition: HMWSoln.h:2890
vector_fp m_Beta2MX_ij_L
Derivative of Beta2_ij[i][j] wrt T.
Definition: HMWSoln.h:2535
XML_Node * findXMLPhase(XML_Node *root, const std::string &idtarget)
Search an XML_Node tree for a named phase XML_Node.
Definition: xml.cpp:1104
vector_fp m_d2lnActCoeffMolaldT2_Unscaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt TT...
Definition: HMWSoln.h:2838
doublereal s_NBS_CLM_lnMolalityActCoeff() const
Calculate the Chlorine activity coefficient on the NBS scale.
Definition: HMWSoln.cpp:5785
doublereal MC_X_o_cutoff_
value of the solvent mole fraction that centers the cutoff polynomials for the cutoff =1 process; ...
Definition: HMWSoln.h:3119
vector_fp IMS_lnActCoeffMolal_
Logarithm of the molal activity coefficients.
Definition: HMWSoln.h:3073
void resize(size_t n, size_t m, doublereal v=0.0)
Resize the array, and fill the new entries with 'v'.
Definition: Array.h:128
Array2D m_Lambda_nj_coeff
Array of coefficients for Lambda_nj[i][j] in the Pitzer/HMW formulation.
Definition: HMWSoln.h:2751
double m_maxIionicStrength
Maximum value of the ionic strength allowed in the calculation of the activity coefficients.
Definition: HMWSoln.h:2344
vector_fp m_Beta1MX_ij_LL
Derivative of Beta1_ij[i][j] wrt TT.
Definition: HMWSoln.h:2503
doublereal m_Mnaught
This is the multiplication factor that goes inside log expressions involving the molalities of specie...
vector_fp m_Mu_nnn_L
Mu coefficient temperature derivative for the self-ternary neutral coefficient.
Definition: HMWSoln.h:2769
virtual doublereal cv_mole() const
Molar heat capacity at constant volume. Units: J/kmol/K.
Definition: HMWSoln.cpp:657
vector_fp m_dlnActCoeffMolaldP_Unscaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt P...
Definition: HMWSoln.h:2852
Array2D m_Beta0MX_ij_coeff
Array of coefficients for Beta0, a variable in Pitzer's papers.
Definition: HMWSoln.h:2480
double ADebye_J(double temperature=-1.0, double pressure=-1.0) const
Return Pitzer's definition of A_J.
Definition: HMWSoln.cpp:1134
virtual void getStandardVolumes(doublereal *vol) const
Get the molar volumes of each species in their standard states at the current T and P of the solution...
virtual doublereal logStandardConc(size_t k=0) const
Returns the natural logarithm of the standard concentration of the kth species.
Definition: HMWSoln.cpp:789
vector_fp m_Beta2MX_ij_P
Derivative of Beta2_ij[i][j] wrt P.
Definition: HMWSoln.h:2547
virtual void getPartialMolarCp(doublereal *cpbar) const
Return an array of partial molar heat capacities for the species in the mixture.
Definition: HMWSoln.cpp:978
vector_fp m_Beta2MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Definition: HMWSoln.h:2529
Array2D m_Lambda_nj_P
Derivative of Lambda_nj[i][j] wrt P.
Definition: HMWSoln.h:2733
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:173
vector_fp m_PhiPhi_IJ
Intermediate variable called PhiPhi in Pitzer's paper.
Definition: HMWSoln.h:3019
vector_fp m_Phi_IJ_L
Derivative of m_Phi_IJ wrt T.
Definition: HMWSoln.h:2995
virtual void getCp_R(doublereal *cpr) const
Get the nondimensional Heat Capacities at constant pressure for the standard state of the species at ...
int IMS_typeCutoff_
IMS Cutoff type.
Definition: HMWSoln.h:3076
Headers for the factory class that can create known ThermoPhase objects (see Thermodynamic Properties...
vector_fp m_dlnActCoeffMolaldT_Scaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt T...
Definition: HMWSoln.h:2817
virtual doublereal satPressure(doublereal T)
saturation pressure
Definition: PDSS.cpp:408
void setMolarDensity(const doublereal conc)
Set the internally stored molar density (kmol/m^3) for the phase.
Definition: HMWSoln.cpp:725
size_t m_indexCLM
Index of the phScale species.
vector_fp m_BMX_IJ_P
Derivative of BMX_IJ wrt P.
Definition: HMWSoln.h:2935
vector_fp m_Beta0MX_ij_L
Derivative of Beta0_ij[i][j] wrt T.
Definition: HMWSoln.h:2460
double elambda[17]
This is elambda, MEC.
Definition: HMWSoln.h:2874
Class XML_Node is a tree-based representation of the contents of an XML file.
Definition: xml.h:100
Implementation of a pressure dependent standard state virtual function for a Pure Water Phase (see Sp...
vector_fp m_dlnActCoeffMolaldP_Scaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt P...
Definition: HMWSoln.h:2845
int debugPrinting()
Return int specifying the amount of debug printing.
Definition: HMWSoln.cpp:5814
vector_fp m_Mu_nnn_P
Mu coefficient pressure derivative for the self-ternary neutral coefficient.
Definition: HMWSoln.h:2787
void s_updateScaling_pHScaling() const
Apply the current phScale to a set of activity Coefficients.
Definition: HMWSoln.cpp:5725
virtual void updateStandardStateThermo() const
Updates the standard state thermodynamic functions at the current T and P of the solution.
Array2D m_Lambda_nj_L
Derivative of Lambda_nj[i][j] wrt T. see m_Lambda_ij.
Definition: HMWSoln.h:2727
void setDensity(const doublereal rho)
Set the internally stored density (kg/m^3) of the phase.
Definition: HMWSoln.cpp:715
vector_fp m_BprimeMX_IJ_P
Derivative of BprimeMX wrt P.
Definition: HMWSoln.h:2959
double AionicRadius(int k=0) const
Reports the ionic radius of the kth species.
Definition: HMWSoln.cpp:1174
vector_fp m_CphiMX_ij_P
Derivative of Cphi_ij[i][j] wrt P.
Definition: HMWSoln.h:2610
void warn_deprecated(const std::string &method, const std::string &extra)
Print a warning indicating that method is deprecated.
Definition: global.cpp:76
vector_fp m_speciesSize
Vector of species sizes.
Definition: Phase.h:729
doublereal s_NBS_CLM_d2lnMolalityActCoeff_dT2() const
Calculate the second temperature derivative of the Chlorine activity coefficient on the NBS scale...
Definition: HMWSoln.cpp:5800
virtual doublereal relative_enthalpy() const
Excess molar enthalpy of the solution from the mixing process.
Definition: HMWSoln.cpp:578
vector_fp m_Psi_ijk_L
Derivative of Psi_ijk[n] wrt T.
Definition: HMWSoln.h:2684
Header for a class used to house several approximation routines for properties of water...
void s_update_d2lnMolalityActCoeff_dT2() const
This function calculates the temperature second derivative of the natural logarithm of the molality a...
Definition: HMWSoln.cpp:3746
vector_fp m_lnActCoeffMolal_Unscaled
Logarithm of the activity coefficients on the molality scale.
Definition: HMWSoln.h:2810
doublereal molarDensity() const
Molar density (kmol/m^3).
Definition: Phase.cpp:597
doublereal * ptrColumn(size_t j)
Return a pointer to the top of column j, columns are contiguous in memory.
Definition: Array.h:363
vector_fp m_Theta_ij
Array of 2D data for Theta_ij[i][j] in the Pitzer/HMW formulation.
Definition: HMWSoln.h:2633
vector_fp m_hfunc_IJ
hfunc, was called gprime in Pitzer's paper.
Definition: HMWSoln.h:2903
double elambda1[17]
This is elambda1, MEC.
Definition: HMWSoln.h:2877
virtual doublereal relative_molal_enthalpy() const
Excess molar enthalpy of the solution from the mixing process on a molality basis.
Definition: HMWSoln.cpp:591
vector_fp m_Beta1MX_ij_P
Derivative of Beta1_ij[i][j] wrt P.
Definition: HMWSoln.h:2509
vector_fp m_BprimeMX_IJ_LL
Derivative of BprimeMX wrt TT.
Definition: HMWSoln.h:2953
vector_fp m_PhiPhi_IJ_LL
Derivative of m_PhiPhi_IJ wrt TT.
Definition: HMWSoln.h:3031
ThermoPhase * duplMyselfAsThermoPhase() const
Duplicator from the ThermoPhase parent class.
Definition: HMWSoln.cpp:544
doublereal IMS_cCut_
Parameter in the polyExp cutoff treatment having to do with rate of exp decay.
Definition: HMWSoln.h:3089
virtual void setTemperature(const doublereal temp)
Set the temperature (K)
Definition: HMWSoln.cpp:731
void getMoleFractions(doublereal *const x) const
Get the species mole fraction vector.
Definition: Phase.cpp:519
vector_fp m_Psi_ijk_LL
Derivative of Psi_ijk[n] wrt TT.
Definition: HMWSoln.h:2690
double m_IionicMolality
Current value of the ionic strength on the molality scale Associated Salts, if present in the mechani...
Definition: HMWSoln.h:2338
vector_fp m_Phi_IJ
Intermediate variable called Phi in Pitzer's paper.
Definition: HMWSoln.h:2989
vector_fp m_BphiMX_IJ
Intermediate variable called BphiMX in Pitzer's paper.
Definition: HMWSoln.h:2965
Array2D m_CphiMX_ij_coeff
Array of coefficients for CphiMX, a parameter in the activity coefficient formulation.
Definition: HMWSoln.h:2619
virtual doublereal intEnergy_mole() const
Molar internal energy. Units: J/kmol.
Definition: HMWSoln.cpp:631
vector_fp m_Beta1MX_ij
Definition: HMWSoln.h:2491
virtual doublereal thermalExpansionCoeff() const
The thermal expansion coefficient.
Definition: HMWSoln.cpp:702
void s_updateScaling_pHScaling_dP() const
Apply the current phScale to a set of derivatives of the activity Coefficients wrt pressure...
Definition: HMWSoln.cpp:5770
vector_fp m_tmpV
vector of size m_kk, used as a temporary holding area.
Definition: HMWSoln.h:2429
virtual doublereal standardConcentration(size_t k=0) const
Return the standard concentration for the kth species.
Definition: HMWSoln.cpp:779
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:101
void getUnscaledMolalityActivityCoefficients(doublereal *acMolality) const
Get the array of unscaled non-dimensional molality based activity coefficients at the current solutio...
Definition: HMWSoln.cpp:845
HMWSoln()
Default Constructor.
Definition: HMWSoln.cpp:31
vector_fp m_lnActCoeffMolal_Scaled
Logarithm of the activity coefficients on the molality scale.
Definition: HMWSoln.h:2800
virtual void getActivities(doublereal *ac) const
Get the array of non-dimensional activities at the current solution temperature, pressure, and solution concentration.
Definition: HMWSoln.cpp:819
doublereal mean_X(const doublereal *const Q) const
Evaluate the mole-fraction-weighted mean of an array Q.
Definition: Phase.cpp:623
doublereal MC_X_o_min_
gamma_o value for the cutoff process at the zero solvent point
Definition: HMWSoln.h:3122
vector_fp m_Phi_IJ_P
Derivative of m_Phi_IJ wrt P.
Definition: HMWSoln.h:3007
PDSS * m_waterSS
Water standard state calculator.
Definition: HMWSoln.h:2414
size_t speciesIndex(const std::string &name) const
Returns the index of a species named 'name' within the Phase object.
Definition: Phase.cpp:229
vector_fp m_Alpha1MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Definition: HMWSoln.h:2568
bool importPhase(XML_Node &phase, ThermoPhase *th, SpeciesThermoFactory *spfactory)
Import a phase information into an empty thermophase object.
virtual void setPressure(doublereal p)
Set the internally stored pressure (Pa) at constant temperature and composition.
Definition: HMWSoln.cpp:676
vector_fp m_Phiprime_IJ
Intermediate variable called Phiprime in Pitzer's paper.
Definition: HMWSoln.h:3013
virtual void getPartialMolarEntropies(doublereal *sbar) const
Returns an array of partial molar entropies of the species in the solution.
Definition: HMWSoln.cpp:915
vector_fp m_h2func_IJ
hfunc2, was called gprime in Pitzer's paper.
Definition: HMWSoln.h:2910
vector_fp m_Psi_ijk_P
Derivative of Psi_ijk[n] wrt P.
Definition: HMWSoln.h:2696
void counterIJ_setup() const
Set up a counter variable for keeping track of symmetric binary interactions amongst the solute speci...
Definition: HMWSoln.cpp:1569
vector_fp m_CMX_IJ_L
Derivative of m_CMX_IJ wrt T.
Definition: HMWSoln.h:3049
virtual void getPartialMolarVolumes(doublereal *vbar) const
Return an array of partial molar volumes for the species in the mixture.
Definition: HMWSoln.cpp:960
bool m_molalitiesAreCropped
Boolean indicating whether the molalities are cropped or are modified.
Definition: HMWSoln.h:2863
vector_fp m_dlnActCoeffMolaldT_Unscaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt T...
Definition: HMWSoln.h:2824
Array2D m_Psi_ijk_coeff
Array of coefficients for Psi_ijk[n] in the Pitzer/HMW formulation.
Definition: