diff --git a/include/cantera/thermo/HMWSoln.h b/include/cantera/thermo/HMWSoln.h
index f6f8477c7..243ca2e54 100644
--- a/include/cantera/thermo/HMWSoln.h
+++ b/include/cantera/thermo/HMWSoln.h
@@ -168,93 +168,6 @@ class WaterProps;
* molality has defined units of gmol kg-1, and therefore the ionic strength has
* units of sqrt(gmol/kg).
*
- * In some instances, from some authors, a different formulation is used for the
- * ionic strength in the equations below. The different formulation is due to
- * the possibility of the existence of weak acids and how association wrt to the
- * weak acid equilibrium relation affects the calculation of the activity
- * coefficients via the assumed value of the ionic strength.
- *
- * If we are to assume that the association reaction doesn't have an effect on
- * the ionic strength, then we will want to consider the associated weak acid as
- * in effect being fully dissociated, when we calculate an effective value for
- * the ionic strength. We will call this calculated value, the stoichiometric
- * ionic strength, \f$ I_s \f$, putting a subscript s to denote it from the more
- * straightforward calculation of \f$ I \f$.
- *
- * \f[
- * I_s = \frac{1}{2} \sum_k{m_k^s z_k^2}
- * \f]
- *
- * Here, \f$ m_k^s \f$ is the value of the molalities calculated assuming that
- * all weak acid-base pairs are in their fully dissociated states. This
- * calculation may be simplified by considering that the weakly associated acid
- * may be made up of two charged species, k1 and k2, each with their own
- * charges, obeying the following relationship:
- *
- * \f[
- * z_k = z_{k1} + z_{k2}
- * \f]
- * Then, we may only need to specify one charge value, say, \f$ z_{k1}\f$, the
- * cation charge number, in order to get both numbers, since we have already
- * specified \f$ z_k \f$ in the definition of original species. Then, the
- * stoichiometric ionic strength may be calculated via the following formula.
- *
- * \f[
- * I_s = \frac{1}{2} \left(\sum_{k,ions}{m_k z_k^2}+
- * \sum_{k,weak_assoc}(m_k z_{k1}^2 + m_k z_{k2}^2) \right)
- * \f]
- *
- * The specification of which species are weakly associated acids is made in the
- * input file via the `stoichIsMods` XML block, where the charge for k1 is also
- * specified. An example is given below:
- *
- * @code
- *
- * NaCl(aq):-1.0
- *
- * @endcode
- *
- * Because we need the concept of a weakly associated acid in order to calculated
- * \f$ I_s \f$ we need to catalog all species in the phase. This is done using
- * the following categories:
- *
- * - `cEST_solvent` Solvent species (neutral)
- * - `cEST_chargedSpecies` Charged species (charged)
- * - `cEST_weakAcidAssociated` Species which can break apart into charged species.
- * It may or may not be charged. These may or
- * may not be be included in the
- * species solution vector.
- * - `cEST_strongAcidAssociated` Species which always breaks apart into charged species.
- * It may or may not be charged. Normally, these
- * aren't included in the speciation vector.
- * - `cEST_polarNeutral` Polar neutral species
- * - `cEST_nonpolarNeutral` Non polar neutral species
- *
- * Polar and non-polar neutral species are differentiated, because some
- * additions to the activity coefficient expressions distinguish between these
- * two types of solutes. This is the so-called salt-out effect.
- *
- * The type of species is specified in the `electrolyteSpeciesType` XML block.
- * Note, this is not considered a part of the specification of the standard
- * state for the species, at this time. Therefore, this information is put under
- * the `activityCoefficient` XML block. An example is given below
- *
- * @code
- *
- * H2L(L):solvent
- * H+:chargedSpecies
- * NaOH(aq):weakAcidAssociated
- * NaCl(aq):strongAcidAssociated
- * NH3(aq):polarNeutral
- * O2(aq):nonpolarNeutral
- *
- * @endcode
- *
- * Much of the species electrolyte type information is inferred from other
- * information in the input file. For example, as species which is charged is
- * given the "chargedSpecies" default category. A neutral solute species is put
- * into the "nonpolarNeutral" category by default.
- *
* ### Specification of the Excess Gibbs Free Energy
*
* Pitzer's formulation may best be represented as a specification of the excess
@@ -1785,18 +1698,6 @@ private:
*/
int m_formPitzerTemp;
- //! Vector containing the electrolyte species type
- /*!
- * The possible types are:
- * - solvent
- * - Charged Species
- * - weakAcidAssociated
- * - strongAcidAssociated
- * - polarNeutral
- * - nonpolarNeutral
- */
- vector_int m_electrolyteSpeciesType;
-
//! a_k = Size of the ionic species in the DH formulation. units = meters
vector_fp m_Aionic;
@@ -1812,11 +1713,6 @@ private:
//! Reference Temperature for the Pitzer formulations.
double m_TempPitzerRef;
- //! Stoichiometric ionic strength on the molality scale. This differs from
- //! m_IionicMolality in the sense that associated salts are treated as
- //! unassociated salts, when calculating the Ionic strength by this method.
- mutable double m_IionicMolalityStoich;
-
public:
/**
* Form of the constant outside the Debye-Huckel term called A. It's
@@ -1884,21 +1780,6 @@ private:
//! vector of size m_kk, used as a temporary holding area.
mutable vector_fp m_tmpV;
- /**
- * Stoichiometric species charge -> This is for calculations of the ionic
- * strength which ignore ion-ion pairing into neutral molecules. The
- * Stoichiometric species charge is the charge of one of the ion that would
- * occur if the species broke into two charged ion pairs.
- *
- * NaCl -> m_speciesCharge_Stoich = -1;
- * HSO4- -> H+ + SO42- = -2
- * -> The other charge is calculated.
- *
- * For species that aren't ion pairs, its equal to the m_speciesCharge[]
- * value.
- */
- vector_fp m_speciesCharge_Stoich;
-
/**
* Array of 2D data used in the Pitzer/HMW formulation. Beta0_ij[i][j] is
* the value of the Beta0 coefficient for the ij salt. It will be nonzero
@@ -2658,13 +2539,6 @@ private:
//! Calculate molality cut-off parameters
void calcMCCutoffParams_();
- //! Utility function to assign an integer value from a string for the
- //! ElectrolyteSpeciesType field.
- /*!
- * @param estString string name of the electrolyte species type
- */
- static int interp_est(const std::string& estString);
-
public:
//! Turn on copious debug printing when this is true
mutable int m_debugCalc;
diff --git a/src/thermo/HMWSoln.cpp b/src/thermo/HMWSoln.cpp
index 7628e33ad..9ac4e6074 100644
--- a/src/thermo/HMWSoln.cpp
+++ b/src/thermo/HMWSoln.cpp
@@ -29,7 +29,6 @@ HMWSoln::HMWSoln() :
m_IionicMolality(0.0),
m_maxIionicStrength(100.0),
m_TempPitzerRef(298.15),
- m_IionicMolalityStoich(0.0),
m_form_A_Debye(A_DEBYE_WATER),
m_A_Debye(1.172576), // units = sqrt(kg/gmol)
m_waterSS(0),
@@ -76,7 +75,6 @@ HMWSoln::HMWSoln(const std::string& inputFile, const std::string& id_) :
m_IionicMolality(0.0),
m_maxIionicStrength(100.0),
m_TempPitzerRef(298.15),
- m_IionicMolalityStoich(0.0),
m_form_A_Debye(A_DEBYE_WATER),
m_A_Debye(1.172576), // units = sqrt(kg/gmol)
m_waterSS(0),
@@ -120,7 +118,6 @@ HMWSoln::HMWSoln(XML_Node& phaseRoot, const std::string& id_) :
m_IionicMolality(0.0),
m_maxIionicStrength(100.0),
m_TempPitzerRef(298.15),
- m_IionicMolalityStoich(0.0),
m_form_A_Debye(A_DEBYE_WATER),
m_A_Debye(1.172576), // units = sqrt(kg/gmol)
m_waterSS(0),
@@ -611,10 +608,7 @@ double HMWSoln::AionicRadius(int k) const
void HMWSoln::initLengths()
{
- // Resize lengths equal to the number of species in the phase.
- m_electrolyteSpeciesType.resize(m_kk, cEST_polarNeutral);
m_speciesSize.resize(m_kk);
- m_speciesCharge_Stoich.resize(m_kk, 0.0);
m_Aionic.resize(m_kk, 0.0);
m_tmpV.resize(m_kk, 0.0);
m_molalitiesCropped.resize(m_kk, 0.0);
@@ -743,20 +737,6 @@ void HMWSoln::s_update_lnMolalityActCoeff() const
// coefficient calculations.
calcMolalitiesCropped();
- // Calculate the stoichiometric ionic charge. This isn't used in the Pitzer
- // formulation.
- m_IionicMolalityStoich = 0.0;
- for (size_t k = 0; k < m_kk; k++) {
- double z_k = charge(k);
- double zs_k1 = m_speciesCharge_Stoich[k];
- if (z_k == zs_k1) {
- m_IionicMolalityStoich += m_molalities[k] * z_k * z_k;
- } else {
- double zs_k2 = z_k - zs_k1;
- m_IionicMolalityStoich += m_molalities[k] * (zs_k1 * zs_k1 + zs_k2 * zs_k2);
- }
- }
-
// Update the temperature dependence of the pitzer coefficients and their
// derivatives
s_updatePitzer_CoeffWRTemp();
diff --git a/src/thermo/HMWSoln_input.cpp b/src/thermo/HMWSoln_input.cpp
index 45a0b74ae..1df8d5384 100644
--- a/src/thermo/HMWSoln_input.cpp
+++ b/src/thermo/HMWSoln_input.cpp
@@ -24,27 +24,6 @@ using namespace std;
namespace Cantera
{
-int HMWSoln::interp_est(const std::string& estString)
-{
- if (ba::iequals(estString, "solvent")) {
- return cEST_solvent;
- } else if (ba::iequals(estString, "chargedspecies")) {
- return cEST_chargedSpecies;
- } else if (ba::iequals(estString, "weakacidassociated")) {
- return cEST_weakAcidAssociated;
- } else if (ba::iequals(estString, "strongacidassociated")) {
- return cEST_strongAcidAssociated;
- } else if (ba::iequals(estString, "polarneutral")) {
- return cEST_polarNeutral;
- } else if (ba::iequals(estString, "nonpolarneutral")) {
- return cEST_nonpolarNeutral;
- }
- int retn, rval;
- if ((retn = sscanf(estString.c_str(), "%d", &rval)) != 1) {
- return -1;
- }
- return rval;
-}
void HMWSoln::readXMLBinarySalt(XML_Node& BinSalt)
{
@@ -1032,13 +1011,6 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
// water calculator.
m_waterProps.reset(new WaterProps(&dynamic_cast(*m_waterSS)));
- // Fill in parameters for the calculation of the stoichiometric Ionic
- // Strength. The default is that stoich charge is the same as the regular
- // charge.
- for (size_t k = 0; k < m_kk; k++) {
- m_speciesCharge_Stoich[k] = charge(k);
- }
-
// Go get all of the coefficients and factors in the activityCoefficients
// XML block
XML_Node* acNodePtr = 0;
@@ -1084,40 +1056,6 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
- // First look at the species database. Look for the subelement
- // "stoichIsMods" in each of the species SS databases.
- std::vector xspecies = speciesData();
- for (size_t k = 0; k < m_kk; k++) {
- size_t jmap = npos;
- string kname = speciesName(k);
- for (size_t j = 0; j < xspecies.size(); j++) {
- const XML_Node& sp = *xspecies[j];
- string jname = sp["name"];
- if (jname == kname) {
- jmap = j;
- break;
- }
- }
- if (jmap != npos) {
- const XML_Node& sp = *xspecies[jmap];
- getOptionalFloat(sp, "stoichIsMods", m_speciesCharge_Stoich[k]);
- }
- }
-
- // Now look at the activity coefficient database
- if (acNodePtr && acNodePtr->hasChild("stoichIsMods")) {
- XML_Node& sIsNode = acNodePtr->child("stoichIsMods");
- map msIs;
- getMap(sIsNode, msIs);
- for (const auto& b : msIs) {
- size_t kk = speciesIndex(b.first);
- if (kk != npos) {
- double val = fpValue(b.second);
- m_speciesCharge_Stoich[kk] = val;
- }
- }
- }
-
// Loop through the children getting multiple instances of parameters
if (acNodePtr) {
for (size_t i = 0; i < acNodePtr->nChildren(); i++) {
@@ -1149,55 +1087,6 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
readXMLCroppingCoefficients(acNode);
}
- // Fill in the vector specifying the electrolyte species type
- //
- // First fill in default values. Everything is either a charge species, a
- // nonpolar neutral, or the solvent.
- for (size_t k = 0; k < m_kk; k++) {
- if (fabs(charge(k)) > 0.0001) {
- m_electrolyteSpeciesType[k] = cEST_chargedSpecies;
- if (fabs(m_speciesCharge_Stoich[k] - charge(k)) > 0.0001) {
- m_electrolyteSpeciesType[k] = cEST_weakAcidAssociated;
- }
- } else if (fabs(m_speciesCharge_Stoich[k]) > 0.0001) {
- m_electrolyteSpeciesType[k] = cEST_weakAcidAssociated;
- } else {
- m_electrolyteSpeciesType[k] = cEST_nonpolarNeutral;
- }
- }
- m_electrolyteSpeciesType[0] = cEST_solvent;
-
- // First look at the species database. Look for the subelement
- // "stoichIsMods" in each of the species SS databases.
- std::vector xspecies = speciesData();
- for (size_t k = 0; k < m_kk; k++) {
- const XML_Node* spPtr = xspecies[k];
- if (spPtr && spPtr->hasChild("electrolyteSpeciesType")) {
- string est = getChildValue(*spPtr, "electrolyteSpeciesType");
- if ((m_electrolyteSpeciesType[k] = interp_est(est)) == -1) {
- throw CanteraError("HMWSoln::initThermoXML",
- "Bad electrolyte type: " + est);
- }
- }
- }
-
- // Then look at the phase thermo specification
- if (acNodePtr && acNodePtr->hasChild("electrolyteSpeciesType")) {
- XML_Node& ESTNode = acNodePtr->child("electrolyteSpeciesType");
- map msEST;
- getMap(ESTNode, msEST);
- for (const auto& b : msEST) {
- size_t kk = speciesIndex(b.first);
- if (kk != npos) {
- string est = b.second;
- if ((m_electrolyteSpeciesType[kk] = interp_est(est)) == -1) {
- throw CanteraError("HMWSoln::initThermoXML",
- "Bad electrolyte type: " + est);
- }
- }
- }
- }
-
IMS_typeCutoff_ = 2;
if (IMS_typeCutoff_ == 2) {
calcIMSCutoffParams_();