IdealMolalSoln.cpp

Go to the documentation of this file.

/**
 *  @file IdealMolalSoln.cpp
 *   ThermoPhase object for the ideal molal equation of
 * state (see \ref thermoprops 
 * and class \link Cantera::IdealMolalSoln IdealMolalSoln\endlink).
 *
 * Definition file for a derived class of ThermoPhase that handles
 * variable pressure standard state methods for calculating
 * thermodynamic properties that are further based upon 
 * activities on the molality scale. The Ideal molal
 * solution assumes that all molality-based activity
 * coefficients are equal to one. This turns out, actually, to be
 * highly nonlinear when the solvent densities get low.
 */
/*
 * Copywrite (2006) Sandia Corporation. Under the terms of 
 * Contract DE-AC04-94AL85000 with Sandia Corporation, the
 * U.S. Government retains certain rights in this software.
 */
/*
 *
 *  $Date: 2009-12-05 14:08:43 -0500 (Sat, 05 Dec 2009) $
 *  $Revision: 279 $
 */

#include "IdealMolalSoln.h"
#include "ThermoFactory.h"
#include <cmath>

//@{
#ifndef MAX
#define MAX(x,y)    (( (x) > (y) ) ? (x) : (y))
#endif
//@}

namespace Cantera {


  static double xxSmall = 1.0E-150;
  /**
   * Default constructor
   */
  IdealMolalSoln::IdealMolalSoln() :
    MolalityVPSSTP(),
    m_formGC(2),
    IMS_typeCutoff_(0),
    IMS_X_o_cutoff_(0.20),
    IMS_gamma_o_min_(0.00001),
    IMS_gamma_k_min_(10.0),
    IMS_cCut_(.05),
    IMS_slopefCut_(0.6),
    IMS_dfCut_(0.0),
    IMS_efCut_(0.0),
    IMS_afCut_(0.0),
    IMS_bfCut_(0.0),
    IMS_slopegCut_(0.0),
    IMS_dgCut_(0.0),
    IMS_egCut_(0.0),
    IMS_agCut_(0.0),
    IMS_bgCut_(0.0)
  {
  }

  /*
   * Copy Constructor:
   *
   *  Note this stuff will not work until the underlying phase
   *  has a working copy constructor
   */
  IdealMolalSoln::IdealMolalSoln(const IdealMolalSoln &b) :
    MolalityVPSSTP(b)
  {
    /*
     * Use the assignment operator to do the brunt
     * of the work for the copy construtor.
     */
    *this = b;
  }

  /*
   * operator=()
   *
   *  Note this stuff will not work until the underlying phase
   *  has a working assignment operator
   */
  IdealMolalSoln& IdealMolalSoln::
  operator=(const IdealMolalSoln &b) {
    if (&b != this) {
      MolalityVPSSTP::operator=(b);
      m_speciesMolarVolume  = b.m_speciesMolarVolume;
      m_formGC              = b.m_formGC;
      IMS_typeCutoff_       = b.IMS_typeCutoff_;
      IMS_X_o_cutoff_       = b.IMS_X_o_cutoff_;
      IMS_gamma_o_min_      = b.IMS_gamma_o_min_;
      IMS_gamma_k_min_      = b.IMS_gamma_k_min_;
      IMS_cCut_             = b.IMS_cCut_;
      IMS_slopefCut_        = b.IMS_slopefCut_;
      IMS_dfCut_            = b.IMS_dfCut_;
      IMS_efCut_            = b.IMS_efCut_;
      IMS_afCut_            = b.IMS_afCut_;
      IMS_bfCut_            = b.IMS_bfCut_;
      IMS_slopegCut_        = b.IMS_slopegCut_;
      IMS_dgCut_            = b.IMS_dgCut_;
      IMS_egCut_            = b.IMS_egCut_;
      IMS_agCut_            = b.IMS_agCut_;
      IMS_bgCut_            = b.IMS_bgCut_;
      m_expg0_RT            = b.m_expg0_RT;
      m_pe                  = b.m_pe;
      m_pp                  = b.m_pp;
      m_tmpV                = b.m_tmpV;
      IMS_lnActCoeffMolal_  = b.IMS_lnActCoeffMolal_;
    }
    return *this;
  }

  IdealMolalSoln::IdealMolalSoln(std::string inputFile, std::string id) :
    MolalityVPSSTP(),
    m_formGC(2),
    IMS_typeCutoff_(0),
    IMS_X_o_cutoff_(0.2),
    IMS_gamma_o_min_(0.00001),
    IMS_gamma_k_min_(10.0),
    IMS_cCut_(.05),
    IMS_slopefCut_(0.6),
    IMS_dfCut_(0.0),
    IMS_efCut_(0.0),
    IMS_afCut_(0.0),
    IMS_bfCut_(0.0),
    IMS_slopegCut_(0.0),
    IMS_dgCut_(0.0),
    IMS_egCut_(0.0),
    IMS_agCut_(0.0),
    IMS_bgCut_(0.0)
  {
    constructPhaseFile(inputFile, id);
  }

  IdealMolalSoln::IdealMolalSoln(XML_Node& root, std::string id) :
    MolalityVPSSTP(),
    m_formGC(2),
    IMS_typeCutoff_(0),
    IMS_X_o_cutoff_(0.2),
    IMS_gamma_o_min_(0.00001),
    IMS_gamma_k_min_(10.0),
    IMS_cCut_(.05),
    IMS_slopefCut_(0.6),
    IMS_dfCut_(0.0),
    IMS_efCut_(0.0),
    IMS_afCut_(0.0),
    IMS_bfCut_(0.0),
    IMS_slopegCut_(0.0),
    IMS_dgCut_(0.0),
    IMS_egCut_(0.0),
    IMS_agCut_(0.0),
    IMS_bgCut_(0.0)
  {
    constructPhaseXML(root, id);
  }

  /*
   *
   * ~IdealMolalSoln():   (virtual)
   *
   * Destructor: does nothing:
   *
   */
  IdealMolalSoln::~IdealMolalSoln() {
  }

  /**
   *
   */
  ThermoPhase* IdealMolalSoln::duplMyselfAsThermoPhase() const {
    IdealMolalSoln* mtp = new IdealMolalSoln(*this);
    return (ThermoPhase *) mtp;
  }

  //
  // -------- Molar Thermodynamic Properties of the Solution --------------- 
  //
  /*
   * Molar enthalpy of the solution: Units: J/kmol.
   *
   * Returns the amount of enthalpy per mole of solution.
   * For an ideal molal solution,
   * \f[
   * \bar{h}(T, P, X_k) = \sum_k X_k \bar{h}_k(T)  
   * \f]
   * The formula is written in terms of the partial molar enthalpies.
   * \f$ \bar{h}_k(T, p, m_k) \f$.
   * See the partial molar enthalpy function, getPartialMolarEnthalpies(),
   * for details.
   *
   * Units: J/kmol
   */
  doublereal IdealMolalSoln::enthalpy_mole() const {
    getPartialMolarEnthalpies(DATA_PTR(m_tmpV));
    getMoleFractions(DATA_PTR(m_pp));
    double val = mean_X(DATA_PTR(m_tmpV));
    return val;
  }

  /*
   * Molar internal energy of the solution: Units: J/kmol.
   *
   * Returns the amount of internal energy per mole of solution.
   * For an ideal molal solution,
   * \f[
   * \bar{u}(T, P, X_k) = \sum_k X_k \bar{u}_k(T)  
   * \f]
   * The formula is written in terms of the partial molar internal energy.
   * \f$ \bar{u}_k(T, p, m_k) \f$.
   */
  doublereal IdealMolalSoln::intEnergy_mole() const {
    getPartialMolarEnthalpies(DATA_PTR(m_tmpV));
    return mean_X(DATA_PTR(m_tmpV));
  }

  /*
   * Molar entropy of the solution: Units J/kmol/K.
   *
   * Returns the amount of entropy per mole of solution.
   * For an ideal molal solution,
   * \f[
   * \bar{s}(T, P, X_k) = \sum_k X_k \bar{s}_k(T)  
   * \f]
   * The formula is written in terms of the partial molar entropies.
   * \f$ \bar{s}_k(T, p, m_k) \f$.
   * See the partial molar entropies function, getPartialMolarEntropies(),
   * for details.
   *
   * Units: J/kmol/K.
   */
  doublereal IdealMolalSoln::entropy_mole() const {
    getPartialMolarEntropies(DATA_PTR(m_tmpV));
    return mean_X(DATA_PTR(m_tmpV));
  }

  /*
   * Molar Gibbs function for the solution: Units J/kmol. 
   *
   * Returns the gibbs free energy of the solution per mole
   * of the solution.
   *
   * \f[
   * \bar{g}(T, P, X_k) = \sum_k X_k \mu_k(T)  
   * \f]
   *
   * Units: J/kmol
   */
  doublereal IdealMolalSoln::gibbs_mole() const {
    getChemPotentials(DATA_PTR(m_tmpV));
    return mean_X(DATA_PTR(m_tmpV));
  }

  /*
   * Molar heat capacity at constant pressure: Units: J/kmol/K. 
   *  * \f[
   * \bar{c}_p(T, P, X_k) = \sum_k X_k \bar{c}_{p,k}(T)  
   * \f]
   *
   * Units: J/kmol/K
   */
  doublereal IdealMolalSoln::cp_mole() const {
    getPartialMolarCp(DATA_PTR(m_tmpV));
    double val = mean_X(DATA_PTR(m_tmpV));
    return val;
  }

  /*
   * Molar heat capacity at constant volume: Units: J/kmol/K. 
   * NOT IMPLEMENTED.
   * Units: J/kmol/K
   */
  doublereal IdealMolalSoln::cv_mole() const {
    return err("not implemented");
  }

  //
  // ------- Mechanical Equation of State Properties ------------------------
  //



  /*
   * Set the pressure at constant temperature. Units: Pa.
   * This method sets a constant within the object.
   * The mass density is not a function of pressure.
   */
  void IdealMolalSoln::setPressure(doublereal p) {
    setState_TP(temperature(), p);
  }

  void IdealMolalSoln::calcDensity() {
    double *vbar = &m_pp[0];
    getPartialMolarVolumes(vbar);
    double *x = &m_tmpV[0];
    getMoleFractions(x);
    doublereal vtotal = 0.0;
    for (int i = 0; i < m_kk; i++) {
      vtotal += vbar[i] * x[i];
    }
    doublereal dd = meanMolecularWeight() / vtotal;
    State::setDensity(dd);
  }

  /*
   * The isothermal compressibility. Units: 1/Pa.
   * The isothermal compressibility is defined as
   * \f[
   * \kappa_T = -\frac{1}{v}\left(\frac{\partial v}{\partial P}\right)_T
   * \f]
   *
   *  It's equal to zero for this model, since the molar volume
   *  doesn't change with pressure or temperature.
   */
  doublereal IdealMolalSoln::isothermalCompressibility() const {
    return 0.0;
  }

  /*
   * The thermal expansion coefficient. Units: 1/K.
   * The thermal expansion coefficient is defined as
   *
   * \f[
   * \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P
   * \f]
   *
   *  It's equal to zero for this model, since the molar volume
   *  doesn't change with pressure or temperature.
   */
  doublereal IdealMolalSoln::thermalExpansionCoeff() const {
    return 0.0;
  }
    
  /*
   * Overwritten setDensity() function is necessary because the
   * density is not an indendent variable.
   *
   * This function will now throw an error condition
   *
   * @internal May have to adjust the strategy here to make
   * the eos for these materials slightly compressible, in order
   * to create a condition where the density is a function of
   * the pressure.
   *
   * This function will now throw an error condition.
   *
   *  NOTE: This is an overwritten function from the State.h
   *        class
   */
  void IdealMolalSoln::setDensity(const doublereal rho) {
    double dens = density();
    if (rho != dens) {
      throw CanteraError("Idea;MolalSoln::setDensity",
                         "Density is not an independent variable");
    }
  }

  /*
   * Overwritten setMolarDensity() function is necessary because the
   * density is not an indendent variable.
   *
   * This function will now throw an error condition.
   *
   *  NOTE: This is a virtual function, overwritten function from the State.h
   *        class
   */
  void IdealMolalSoln::setMolarDensity(const doublereal conc) {
    double concI = State::molarDensity();
    if (conc != concI) {
      throw CanteraError("IdealMolalSoln::setMolarDensity",
                         "molarDensity/denisty is not an independent variable");
    }
  }

  void IdealMolalSoln::setState_TP(doublereal temp, doublereal pres) {
    State::setTemperature(temp);
    m_Pcurrent = pres;
    updateStandardStateThermo();
    //m_densWaterSS = m_waterSS->density();
    calcDensity();
  }

  //
  // ------- Activities and Activity Concentrations
  //

  /*
   * This method returns an array of activity concentrations \f$ C^a_k\f$.
   * \f$ C^a_k\f$ are defined such that 
   * \f$ a_k = C^a_k / C^s_k, \f$ where \f$ C^s_k \f$ 
   * is a standard concentration
   * defined below.  These activity concentrations are used
   * by kinetics manager classes to compute the forward and
   * reverse rates of elementary reactions. 
   *
   * @param c Array of activity concentrations. The 
   *           units depend upon the implementation of the
   *           reaction rate expressions within the phase.
   */
  void IdealMolalSoln::getActivityConcentrations(doublereal* c) const {
    if (m_formGC != 1) {
      double c_solvent = standardConcentration();
      getActivities(c);
      for (int k = 0; k < m_kk; k++) {
        c[k] *= c_solvent;
      }
    } else {
      getActivities(c);
      for (int k = 0; k < m_kk; k++) {
        double c0 = standardConcentration(k);
        c[k] *= c0;
      }
    }
  }

  /*
   * The standard concentration \f$ C^s_k \f$ used to normalize
   * the activity concentration. In many cases, this quantity
   * will be the same for all species in a phase - for example,
   * for an ideal gas \f$ C^s_k = P/\hat R T \f$. For this
   * reason, this method returns a single value, instead of an
   * array.  However, for phases in which the standard
   * concentration is species-specific (e.g. surface species of
   * different sizes), this method may be called with an
   * optional parameter indicating the species.
   *
   */
  doublereal IdealMolalSoln::standardConcentration(int k) const {
    double c0 = 1.0, mvSolvent;
    switch (m_formGC) {
    case 0:
      break;
    case 1:
      c0 = 1.0 /m_speciesMolarVolume[m_indexSolvent];
      break;
    case 2:
      mvSolvent = m_speciesMolarVolume[m_indexSolvent];
      c0 = 1.0 / mvSolvent;
      break;
    }
    return c0;
  }
    
  /*
   * Returns the natural logarithm of the standard 
   * concentration of the kth species
   */
  doublereal IdealMolalSoln::logStandardConc(int k) const {
    double c0 = standardConcentration(k);
    return log(c0);
  }
    
  /*
   * Returns the units of the standard and general concentrations
   * Note they have the same units, as their divisor is 
   * defined to be equal to the activity of the kth species
   * in the solution, which is unitless.
   *
   * This routine is used in print out applications where the
   * units are needed. Usually, MKS units are assumed throughout
   * the program and in the XML input files. 
   *
   * On return uA contains the powers of the units (MKS assumed)
   * of the standard concentrations and generalized concentrations
   * for the kth species.
   *
   *  uA[0] = kmol units - default  = 1
   *  uA[1] = m    units - default  = -nDim(), the number of spatial
   *                                dimensions in the Phase class.
   *  uA[2] = kg   units - default  = 0;
   *  uA[3] = Pa(pressure) units - default = 0;
   *  uA[4] = Temperature units - default = 0;
   *  uA[5] = time units - default = 0
   */
  void IdealMolalSoln::getUnitsStandardConc(double *uA, int k, int sizeUA) const {
    int eos = eosType();
    if (eos == 0) {
      for (int i = 0; i < sizeUA; i++) {
        uA[i] = 0.0;
      }
    } else {
      for (int i = 0; i < sizeUA; i++) {
        if (i == 0) uA[0] = 1.0;
        if (i == 1) uA[1] = -nDim();
        if (i == 2) uA[2] = 0.0;
        if (i == 3) uA[3] = 0.0;
        if (i == 4) uA[4] = 0.0;
        if (i == 5) uA[5] = 0.0;
      }
    }
  }

  /*
   * Get the array of non-dimensional molality-based
   * activities at the current solution temperature, 
   * pressure, and solution concentration.
   *
   *  The max against xmolSolventMIN is to limit the activity
   *  coefficient to be finite as the solvent mf goes to zero.
   */
  void IdealMolalSoln::getActivities(doublereal* ac) const {  
    _updateStandardStateThermo();
    /*
     * Update the molality array, m_molalities()
     *   This requires an update due to mole fractions
     */
    if (IMS_typeCutoff_ == 0) {
      calcMolalities();
      for (int k = 0; k < m_kk; k++) {
        ac[k] = m_molalities[k];
      }
      double xmolSolvent = moleFraction(m_indexSolvent);
      xmolSolvent = fmaxx(m_xmolSolventMIN, xmolSolvent);
      ac[m_indexSolvent] = 
        exp((xmolSolvent - 1.0)/xmolSolvent);
    } else {

      s_updateIMS_lnMolalityActCoeff();
      /*
       * Now calculate the array of activities.
       */
      for (int k = 1; k < m_kk; k++) {
          ac[k] = m_molalities[k] * exp(IMS_lnActCoeffMolal_[k]);
      }
      double xmolSolvent = moleFraction(m_indexSolvent);
      ac[m_indexSolvent] =
      exp(IMS_lnActCoeffMolal_[m_indexSolvent]) * xmolSolvent;
     
    }
  }

  /*
   * Get the array of non-dimensional Molality based
   * activity coefficients at
   * the current solution temperature, pressure, and
   * solution concentration.
   * See Denbigh
   * (note solvent activity coefficient is on the molar scale).
   *
   *  The max against xmolSolventMIN is to limit the activity
   *  coefficient to be finite as the solvent mf goes to zero.
   */
  void IdealMolalSoln::
  getMolalityActivityCoefficients(doublereal* acMolality) const {
    if (IMS_typeCutoff_ == 0) {
      for (int k = 0; k < m_kk; k++) {
        acMolality[k] = 1.0;
      }
      double xmolSolvent = moleFraction(m_indexSolvent);
      xmolSolvent = fmaxx(m_xmolSolventMIN, xmolSolvent);
      acMolality[m_indexSolvent] = 
        exp((xmolSolvent - 1.0)/xmolSolvent) / xmolSolvent;
    } else {
      s_updateIMS_lnMolalityActCoeff();
      std::copy(IMS_lnActCoeffMolal_.begin(), IMS_lnActCoeffMolal_.end(), acMolality);
      for (int k = 0; k < m_kk; k++) {
        acMolality[k] = exp(acMolality[k]);
      }
    }
  }

  //
  // ------ Partial Molar Properties of the Solution -----------------
  //

  /*
   * Get the species chemical potentials: Units: J/kmol.
   *
   * This function returns a vector of chemical potentials of the 
   * species in solution.
   *
   * \f[
   *    \mu_k = \mu^{o}_k(T,P) + R T \ln(\frac{m_k}{m^\Delta})
   * \f]
   * \f[
   *    \mu_w = \mu^{o}_w(T,P) +
   *            R T ((X_w - 1.0) / X_w)
   * \f]
   *
   * \f$ w \f$ refers to the solvent species. 
   * \f$ X_w \f$ is the mole fraction of the solvent.
   * \f$ m_k \f$ is the molality of the kth solute.
   * \f$ m^\Delta is 1 gmol solute per kg solvent. \f$
   *
   * Units: J/kmol.
   */
  void IdealMolalSoln::getChemPotentials(doublereal* mu) const{
    double xx;
    //const double xxSmall = 1.0E-150; 

    // Assertion is made for speed
    AssertThrow(m_indexSolvent == 0, "solvent not the first species");
  
    /*
     * First get the standard chemical potentials
     *  -> this requires updates of standard state as a function
     *     of T and P
     * These are defined at unit molality.
     */
    getStandardChemPotentials(mu);
    /*
     * Update the molality array, m_molalities()
     *   This requires an update due to mole fractions
     */
    calcMolalities();
    /*
     * get the solvent mole fraction
     */
    double xmolSolvent = moleFraction(m_indexSolvent);
    /*
     *   
     */
    doublereal RT = GasConstant * temperature();

    if (IMS_typeCutoff_ == 0 || xmolSolvent > 3.* IMS_X_o_cutoff_/2.0) {
  
      for (int k = 1; k < m_kk; k++) {
        xx = fmaxx(m_molalities[k], xxSmall);
        mu[k] += RT * log(xx);
      }
      /*
       * Do the solvent 
       *  -> see my notes
       */
   
      xx = fmaxx(xmolSolvent, xxSmall);
      mu[m_indexSolvent] += 
        (RT * (xmolSolvent - 1.0) / xx);
    } else {
      /*
       * Update the activity coefficients
       * This also updates the internal molality array.
       */
      s_updateIMS_lnMolalityActCoeff();


      for (int k = 1; k < m_kk; k++) {
        xx = MAX(m_molalities[k], xxSmall);
        mu[k] += RT * (log(xx) + IMS_lnActCoeffMolal_[k]);
      }
      xx = MAX(xmolSolvent, xxSmall);
      mu[m_indexSolvent] +=
        RT * (log(xx) + IMS_lnActCoeffMolal_[m_indexSolvent]);
    }

  }

  /*
   * Returns an array of partial molar enthalpies for the species
   * in the mixture: Units (J/kmol).
   *
   */
  void IdealMolalSoln::getPartialMolarEnthalpies(doublereal* hbar) const {
    getEnthalpy_RT(hbar);
    doublereal RT = _RT();
    for (int k = 0; k < m_kk; k++) {
      hbar[k] *= RT;
    }
  }

  /*
   * Returns an array of partial molar entropies of the species in the
   * solution: Units: J/kmol.
   *
   * Maxwell's equations provide an insight in how to calculate this
   * (p.215 Smith and Van Ness)
   * \f[
   *      \frac{d(\mu_k)}{dT} = -\bar{s}_i
   * \f]
   * For this phase, the partial molar entropies are equal to the
   * standard state species entropies plus the ideal molal solution contribution.
   * 
   * \f[
   *   \bar{s}_k(T,P) =  s^0_k(T) - R log( m_k )
   * \f]
   * \f[
   *   \bar{s}_w(T,P) =  s^0_w(T) - R ((X_w - 1.0) / X_w)
   * \f]
   *
   * The subscript, w, refers to the solvent species. \f$ X_w \f$ is
   * the mole fraction of solvent.
   * The reference-state pure-species entropies,\f$ s^0_k(T) \f$,
   * at the reference pressure, \f$ P_{ref} \f$, are computed by the
   * species thermodynamic
   * property manager. They are polynomial functions of temperature.
   * @see SpeciesThermo
   */
  void IdealMolalSoln::getPartialMolarEntropies(doublereal* sbar) const {
    getEntropy_R(sbar);
    doublereal R = GasConstant;
    doublereal mm;
    calcMolalities();
    if (IMS_typeCutoff_ == 0) {
      for (int k = 0; k < m_kk; k++) {
        if (k != m_indexSolvent) {
          mm = fmaxx(SmallNumber, m_molalities[k]);
          sbar[k] -= R * log(mm);
        }
      }
      double xmolSolvent = moleFraction(m_indexSolvent);
      sbar[m_indexSolvent] -= (R * (xmolSolvent - 1.0) / xmolSolvent);
    } else {
      /*
       * Update the activity coefficients, This also update the
       * internally stored molalities.
       */
      s_updateIMS_lnMolalityActCoeff();
      /*
       * First we will add in the obvious dependence on the T
       * term out front of the log activity term
       */
      doublereal mm;
      for (int k = 0; k < m_kk; k++) {
        if (k != m_indexSolvent) {
          mm = fmaxx(SmallNumber, m_molalities[k]);
          sbar[k] -= R * (log(mm) + IMS_lnActCoeffMolal_[k]);
        }
      }
      double xmolSolvent = moleFraction(m_indexSolvent);
      mm = fmaxx(SmallNumber, xmolSolvent);
      sbar[m_indexSolvent] -= R *(log(mm) + IMS_lnActCoeffMolal_[m_indexSolvent]);

    }
  }

  /*
   * Returns an array of partial molar volumes of the species
   * in the solution: Units: m^3 kmol-1.
   *
   * For this solution, the partial molar volumes are equal to the
   * constant species molar volumes.
   *
   * Units: m^3 kmol-1.
   */
  void IdealMolalSoln::getPartialMolarVolumes(doublereal* vbar) const {
    getStandardVolumes(vbar);
  }

  /*
   * Partial molar heat capacity of the solution: Units: J/kmol/K.
   *
   *   The kth partial molar heat capacity is equal to 
   *   the temperature derivative of the partial molar
   *   enthalpy of the kth species in the solution at constant
   *   P and composition (p. 220 Smith and Van Ness).
   *  \f[
   *  \bar{Cp}_k(T,P) =  {Cp}^0_k(T) 
   * \f]
   *
   *   For this solution, this is equal to the reference state
   *   heat capacities.
   *
   *  Units: J/kmol/K
   */
  void IdealMolalSoln::getPartialMolarCp(doublereal* cpbar) const {
    /*
     * Get the nondimensional gibbs standard state of the
     * species at the T and P of the solution.
     */
    getCp_R(cpbar);
        
    for (int k = 0; k < m_kk; k++) {
      cpbar[k] *= GasConstant;
    }
  }
        
  /*
   * -------- Properties of the Standard State of the Species
   *           in the Solution ------------------
   */



  /*
   * ------ Thermodynamic Values for the Species Reference States ---
   */

  // -> This is handled by VPStandardStatesTP

  /*
   *  -------------- Utilities -------------------------------
   */

  /*
   *  Initialization routine for an IdealMolalSoln phase.
   *
   * This is a virtual routine. This routine will call initThermo()
   * for the parent class as well.
   */
  void IdealMolalSoln::initThermo() {
    initLengths();
    MolalityVPSSTP::initThermo();
  }

  /*
   * Initialization of an IdealMolalSoln phase using an
   * xml file
   *
   * This routine is a precursor to constructPhaseFile(XML_Node*)
   * routine, which does most of the work.
   *
   * @param inputFile XML file containing the description of the
   *        phase
   *
   * @param id  Optional parameter identifying the name of the
   *            phase. If none is given, the first XML
   *            phase element will be used.
   */
  void IdealMolalSoln::constructPhaseFile(std::string inputFile, 
                                          std::string id) {

    if (inputFile.size() == 0) {
      throw CanteraError("IdealMolalSoln::constructPhaseFile",
                         "input file is null");
    }
    std::string path = findInputFile(inputFile);
    std::ifstream fin(path.c_str());
    if (!fin) {
      throw CanteraError("IdealMolalSoln::constructPhaseFile",
                         "could not open "
                         +path+" for reading.");
    }
    /*
     * The phase object automatically constructs an XML object.
     * Use this object to store information.
     */
    XML_Node &phaseNode_XML = xml();
    XML_Node *fxml = new XML_Node();
    fxml->build(fin);
    XML_Node *fxml_phase = findXMLPhase(fxml, id);
    if (!fxml_phase) {
      throw CanteraError("IdealMolalSoln::constructPhaseFile",
                         "ERROR: Can not find phase named " +
                         id + " in file named " + inputFile);
    }
    fxml_phase->copy(&phaseNode_XML);
    constructPhaseXML(*fxml_phase, id);
    delete fxml;
  }
        
  /*
   *   Import and initialize an IdealMolalSoln phase 
   *   specification in an XML tree into the current object.
   *   Here we read an XML description of the phase.
   *   We import descriptions of the elements that make up the
   *   species in a phase.
   *   We import information about the species, including their
   *   reference state thermodynamic polynomials. We then freeze
   *   the state of the species.
   *
   *   Then, we read the species molar volumes from the xml 
   *   tree to finish the initialization.
   *
   * @param phaseNode This object must be the phase node of a
   *             complete XML tree
   *             description of the phase, including all of the
   *             species data. In other words while "phase" must
   *             point to an XML phase object, it must have
   *             sibling nodes "speciesData" that describe
   *             the species in the phase.
   * @param id   ID of the phase. If nonnull, a check is done
   *             to see if phaseNode is pointing to the phase
   *             with the correct id. 
   */
  void IdealMolalSoln::constructPhaseXML(XML_Node& phaseNode, 
                                         std::string id) {
    if (id.size() > 0) {
      std::string idp = phaseNode.id();
      if (idp != id) {
        throw CanteraError("IdealMolalSoln::constructPhaseXML", 
                           "phasenode and Id are incompatible");
      }
    }

    /*
     * Find the Thermo XML node 
     */
    if (!phaseNode.hasChild("thermo")) {
      throw CanteraError("IdealMolalSoln::constructPhaseXML",
                         "no thermo XML node");
    }
 
    /*
     * Call the Cantera importPhase() function. This will import
     * all of the species into the phase. Then, it will call
     * initThermoXML() below.
     */
    bool m_ok = importPhase(phaseNode, this);
    if (!m_ok) {
      throw CanteraError("IdealMolalSoln::constructPhaseXML","importPhase failed "); 
    }
  }

  /*
   *   Import and initialize an IdealMolalSoln phase 
   *   specification in an XML tree into the current object.
   *
   *   This routine is called from importPhase() to finish
   *   up the initialization of the thermo object. It reads in the
   *   species molar volumes.
   *
   * @param phaseNode This object must be the phase node of a
   *             complete XML tree
   *             description of the phase, including all of the
   *             species data. In other words while "phase" must
   *             point to an XML phase object, it must have
   *             sibling nodes "speciesData" that describe
   *             the species in the phase.
   * @param id   ID of the phase. If nonnull, a check is done
   *             to see if phaseNode is pointing to the phase
   *             with the correct id.
   */
  void IdealMolalSoln::initThermoXML(XML_Node& phaseNode, std::string id) {

    /*
     * Initialize the whole thermo object, using a virtual function.
     */
    initThermo();

    if (id.size() > 0) {
      std::string idp = phaseNode.id();
      if (idp != id) {
        throw CanteraError("IdealMolalSoln::initThermo", 
                           "phasenode and Id are incompatible");
      }
    }

    /*
     * Find the Thermo XML node 
     */
    if (!phaseNode.hasChild("thermo")) {
      throw CanteraError("IdealMolalSoln::initThermo",
                         "no thermo XML node");
    }
    XML_Node& thermoNode = phaseNode.child("thermo");

    /*
     * Possible change the form of the standard concentrations
     */
    if (thermoNode.hasChild("standardConc")) {
      XML_Node& scNode = thermoNode.child("standardConc");
      m_formGC = 2;
      std::string formString = scNode.attrib("model");
      if (formString != "") {
        if (formString == "unity") {
          m_formGC = 0;
        } else if (formString == "molar_volume") {
          m_formGC = 1;
        } else if (formString == "solvent_volume") {
          m_formGC = 2;
        } else {
          throw CanteraError("IdealMolalSoln::initThermo",
                             "Unknown standardConc model: " + formString);
        }
      }
    }

    /*
     * Get the Name of the Solvent:
     *      <solvent> solventName </solvent>
     */
    std::string solventName = "";
    if (thermoNode.hasChild("solvent")) {
      XML_Node& scNode = thermoNode.child("solvent");
      std::vector<std::string> nameSolventa;
      getStringArray(scNode, nameSolventa);
      int nsp = static_cast<int>(nameSolventa.size());
      if (nsp != 1) {
        throw CanteraError("IdealMolalSoln::initThermoXML",
                           "badly formed solvent XML node");
      }
      solventName = nameSolventa[0];
    }

    if (thermoNode.hasChild("activityCoefficients")) {
      XML_Node& acNode = thermoNode.child("activityCoefficients");
      std::string modelString = acNode.attrib("model");
      IMS_typeCutoff_ = 0;
      if (modelString != "IdealMolalSoln") {
        throw CanteraError("IdealMolalSoln::initThermoXML",
                           "unknown ActivityCoefficient model: " + modelString);
      }
      if (acNode.hasChild("idealMolalSolnCutoff")) {
        XML_Node& ccNode = acNode.child("idealMolalSolnCutoff");
        modelString = ccNode.attrib("model");
        if (modelString != "") {
          if (modelString == "polyExp") {
            IMS_typeCutoff_ = 2;
          } else if (modelString == "poly") {
            IMS_typeCutoff_ = 1;
          } else {
            throw CanteraError("IdealMolalSoln::initThermoXML",
                               "Unknown  idealMolalSolnCutoff form: " + modelString);
          }

          if (ccNode.hasChild("gamma_o_limit")) {
            IMS_gamma_o_min_ = getFloat(ccNode, "gamma_o_limit");
          }       
          if (ccNode.hasChild("gamma_k_limit")) {
            IMS_gamma_k_min_ = getFloat(ccNode, "gamma_k_limit");
          }
          if (ccNode.hasChild("X_o_cutoff")) {
            IMS_X_o_cutoff_ = getFloat(ccNode, "X_o_cutoff");
          }
          if (ccNode.hasChild("c_0_param")) {
            IMS_cCut_ = getFloat(ccNode, "c_0_param");
          }
          if (ccNode.hasChild("slope_f_limit")) {
            IMS_slopefCut_ = getFloat(ccNode, "slope_f_limit");
          }
          if (ccNode.hasChild("slope_g_limit")) {
            IMS_slopegCut_ = getFloat(ccNode, "slope_g_limit");
          }

        }
      }
    }


    /*
     * Reconcile the solvent name and index.
     */
    for (int k = 0; k < m_kk; k++) {
      std::string sname = speciesName(k);
      if (solventName == sname) {
        m_indexSolvent = k;
        break;
      }
    }
    if (m_indexSolvent == -1) {
      std::cout << "IdealMolalSoln::initThermo: Solvent Name not found" 
           << std::endl;
      throw CanteraError("IdealMolalSoln::initThermo",
                         "Solvent name not found");
    }
    if (m_indexSolvent != 0) {
      throw CanteraError("IdealMolalSoln::initThermo",
                         "Solvent " + solventName +
                         " should be first species");
    }


    /*
     * Now go get the molar volumes
     */
    XML_Node& speciesList = phaseNode.child("speciesArray");
    XML_Node* speciesDB =
      get_XML_NameID("speciesData", speciesList["datasrc"],
                     &phaseNode.root());
    const std::vector<std::string> &sss = speciesNames();

    for (int k = 0; k < m_kk; k++) {
      XML_Node* s =  speciesDB->findByAttr("name", sss[k]);
      XML_Node *ss = s->findByName("standardState");
      m_speciesMolarVolume[k] = getFloat(*ss, "molarVolume", "toSI");
    }

    IMS_typeCutoff_ = 2;
    if (IMS_typeCutoff_ == 2) {
      calcIMSCutoffParams_();
    }

    MolalityVPSSTP::initThermoXML(phaseNode, id);


    setMoleFSolventMin(1.0E-5);
    /*
     * Set the state
     */
    if (phaseNode.hasChild("state")) {
      XML_Node& stateNode = phaseNode.child("state");
      setStateFromXML(stateNode);
    }
    
  }

  /*
   * @internal
   * Set equation of state parameters. The number and meaning of
   * these depends on the subclass. 
   * @param n number of parameters
   * @param c array of \i n coefficients
   * 
   */
  void IdealMolalSoln::setParameters(int n, doublereal* const c) {
  }

  void IdealMolalSoln::getParameters(int &n, doublereal * const c) const {
  }

  /*
   * Set equation of state parameter values from XML
   * entries. This method is called by function importPhase in
   * file importCTML.cpp when processing a phase definition in
   * an input file. It should be overloaded in subclasses to set
   * any parameters that are specific to that particular phase
   * model.
   *
   * @param eosdata An XML_Node object corresponding to
   * the "thermo" entry for this phase in the input file.
   *
   * HKM -> Right now, the parameters are set elsewhere (initThermo)
   *        It just didn't seem to fit.
   */
  void IdealMolalSoln::setParametersFromXML(const XML_Node& eosdata) {
  }

  /*
   * ----------- Critical State Properties --------------------------
   */

  /*
   * ------------ Private and Restricted Functions ------------------
   */

  /*
   * Bail out of functions with an error exit if they are not
   * implemented.
   */
  doublereal IdealMolalSoln::err(std::string msg) const {
    throw CanteraError("IdealMolalSoln",
                       "Unfinished func called: " + msg );
    return 0.0;
  }

  

  // This function will be called to update the internally storred
  // natural logarithm of the molality activity coefficients
  /*
   * Normally they are all one. However, sometimes they are not,
   * due to stability schemes
   *
   *    gamma_k_molar =  gamma_k_molal / Xmol_solvent
   *
   *    gamma_o_molar = gamma_o_molal
   */
  void  IdealMolalSoln::s_updateIMS_lnMolalityActCoeff() const {
    int k;
    double tmp;
    /*
     * Calculate the molalities. Currently, the molalities
     * may not be current with respect to the contents of the
     * State objects' data.
     */
    calcMolalities();

    double xmolSolvent = moleFraction(m_indexSolvent);
    double xx = MAX(m_xmolSolventMIN, xmolSolvent);

    if (IMS_typeCutoff_ == 0) {
      for (k = 1; k < m_kk; k++) {
        IMS_lnActCoeffMolal_[k]= 0.0;
      }
      IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
      return;
    } else if (IMS_typeCutoff_ == 1) {
      if (xmolSolvent > 3.0 * IMS_X_o_cutoff_/2.0 ) {
        for (k = 1; k < m_kk; k++) {
          IMS_lnActCoeffMolal_[k]= 0.0;
        }
        IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
        return;
      } else if  (xmolSolvent < IMS_X_o_cutoff_/2.0) {  
        tmp = log(xx * IMS_gamma_k_min_);
        for (k = 1; k < m_kk; k++) {
          IMS_lnActCoeffMolal_[k]= tmp;
        }
        IMS_lnActCoeffMolal_[m_indexSolvent] = log(IMS_gamma_o_min_);
        return;
      } else {
      
        /*
         * If we are in the middle region, calculate the connecting polynomials
         */
        double xminus  = xmolSolvent - IMS_X_o_cutoff_/2.0;
        double xminus2 = xminus * xminus;
        double xminus3 = xminus2 * xminus;
        double x_o_cut2 = IMS_X_o_cutoff_ * IMS_X_o_cutoff_;
        double x_o_cut3 =  x_o_cut2 * IMS_X_o_cutoff_;
    
        double h2 = 3.5 * xminus2 /  IMS_X_o_cutoff_ - 2.0 * xminus3 / x_o_cut2;    
        double h2_prime = 7.0 * xminus /  IMS_X_o_cutoff_ - 6.0 * xminus2 /  x_o_cut2;

        double h1 =   (1.0 - 3.0 * xminus2 /  x_o_cut2 + 2.0 *  xminus3/ x_o_cut3);
        double h1_prime = (- 6.0 * xminus /  x_o_cut2 + 6.0 *  xminus2/ x_o_cut3);

        double h1_g = h1 / IMS_gamma_o_min_;
        double h1_g_prime  = h1_prime / IMS_gamma_o_min_;

        double alpha = 1.0 / ( exp(1.0) * IMS_gamma_k_min_);
        double h1_f = h1 * alpha;
        double h1_f_prime  = h1_prime * alpha;

        double f = h2 + h1_f;
        double f_prime = h2_prime + h1_f_prime;

        double g = h2 + h1_g;
        double g_prime = h2_prime + h1_g_prime;

        tmp = (xmolSolvent/ g * g_prime + (1.0-xmolSolvent) / f * f_prime);
        double lngammak = -1.0 - log(f) + tmp * xmolSolvent;
        double lngammao =-log(g) - tmp * (1.0-xmolSolvent);
      
        tmp = log(xmolSolvent) + lngammak;
        for (k = 1; k < m_kk; k++) {
          IMS_lnActCoeffMolal_[k]= tmp;
        }
        IMS_lnActCoeffMolal_[m_indexSolvent] = lngammao;
      }
    } 

    // Exponentials - trial 2
    else if (IMS_typeCutoff_ == 2) {
      if (xmolSolvent > IMS_X_o_cutoff_) {
        for (k = 1; k < m_kk; k++) {
          IMS_lnActCoeffMolal_[k]= 0.0;
        }
        IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
        return;
      } else {
     
        double xoverc = xmolSolvent/IMS_cCut_;
        double eterm = std::exp(-xoverc);
       
        double fptmp = IMS_bfCut_  - IMS_afCut_ / IMS_cCut_ - IMS_bfCut_*xoverc
          + 2.0*IMS_dfCut_*xmolSolvent - IMS_dfCut_*xmolSolvent*xoverc;
        double f_prime = 1.0 + eterm*fptmp;     
        double f = xmolSolvent + IMS_efCut_ + eterm * (IMS_afCut_ + xmolSolvent * (IMS_bfCut_ + IMS_dfCut_*xmolSolvent));

        double gptmp = IMS_bgCut_  - IMS_agCut_ / IMS_cCut_ - IMS_bgCut_*xoverc
          + 2.0*IMS_dgCut_*xmolSolvent - IMS_dgCut_*xmolSolvent*xoverc;
        double g_prime = 1.0 + eterm*gptmp;
        double g = xmolSolvent + IMS_egCut_ + eterm * (IMS_agCut_ + xmolSolvent * (IMS_bgCut_ + IMS_dgCut_*xmolSolvent));

        tmp = (xmolSolvent / g * g_prime + (1.0 - xmolSolvent) / f * f_prime);
        double lngammak = -1.0 - log(f) + tmp * xmolSolvent;
        double lngammao =-log(g) - tmp * (1.0-xmolSolvent);
      
        tmp = log(xx) + lngammak;
        for (k = 1; k < m_kk; k++) {
          IMS_lnActCoeffMolal_[k]= tmp;
        }
        IMS_lnActCoeffMolal_[m_indexSolvent] = lngammao;
      }
    }
    return;
  }

  /*
   * This internal function adjusts the lengths of arrays.
   *
   * This function is not virtual nor is it inherited
   */
  void IdealMolalSoln::initLengths() {
    m_kk = nSpecies();
    /*
     * Obtain the limits of the temperature from the species
     * thermo handler's limits.
     */
    int leng = m_kk;
    m_expg0_RT.resize(leng);
    m_pe.resize(leng, 0.0);
    m_pp.resize(leng);
    m_speciesMolarVolume.resize(leng);
    m_tmpV.resize(leng);
    IMS_lnActCoeffMolal_.resize(leng);
  }


  void  IdealMolalSoln::calcIMSCutoffParams_() {
    IMS_afCut_ = 1.0 / (std::exp(1.0) *  IMS_gamma_k_min_);
    IMS_efCut_ = 0.0;
    bool converged = false;
    double oldV = 0.0;
    int its;
    for (its = 0; its < 100 && !converged; its++) {
      oldV = IMS_efCut_;
      IMS_afCut_ = 1.0 / (std::exp(1.0) * IMS_gamma_k_min_)  - IMS_efCut_;
      IMS_bfCut_ = IMS_afCut_ / IMS_cCut_ + IMS_slopefCut_ - 1.0;
      IMS_dfCut_ = ((- IMS_afCut_/IMS_cCut_ + IMS_bfCut_ - IMS_bfCut_*IMS_X_o_cutoff_/IMS_cCut_) 
                    /
                (IMS_X_o_cutoff_*IMS_X_o_cutoff_/IMS_cCut_ - 2.0 * IMS_X_o_cutoff_));
      double tmp = IMS_afCut_ + IMS_X_o_cutoff_*( IMS_bfCut_ + IMS_dfCut_ * IMS_X_o_cutoff_);
      double eterm = std::exp(-IMS_X_o_cutoff_/IMS_cCut_);
      IMS_efCut_ = - eterm * (tmp);
      if (fabs(IMS_efCut_ - oldV) < 1.0E-14) {
        converged = true;
      } 
    }
    if (!converged) {
      throw CanteraError(" IdealMolalSoln::calcCutoffParams_()",
                         " failed to converge on the f polynomial");
    }
    converged = false;
    double f_0 = IMS_afCut_ + IMS_efCut_;
    double f_prime_0 = 1.0 - IMS_afCut_ / IMS_cCut_ + IMS_bfCut_;
    IMS_egCut_ = 0.0;
    for (its = 0; its < 100 && !converged; its++) {
      oldV = IMS_egCut_;
      double lng_0 = -log(IMS_gamma_o_min_) -  f_prime_0 / f_0;
      IMS_agCut_ = exp(lng_0) - IMS_egCut_;
      IMS_bgCut_ = IMS_agCut_ / IMS_cCut_ + IMS_slopegCut_ - 1.0;
      IMS_dgCut_ = ((- IMS_agCut_/IMS_cCut_ + IMS_bgCut_ - IMS_bgCut_*IMS_X_o_cutoff_/IMS_cCut_) 
                    / 
                    (IMS_X_o_cutoff_*IMS_X_o_cutoff_/IMS_cCut_ - 2.0 * IMS_X_o_cutoff_));
      double tmp = IMS_agCut_ + IMS_X_o_cutoff_*( IMS_bgCut_ + IMS_dgCut_ *IMS_X_o_cutoff_);
      double eterm = std::exp(-IMS_X_o_cutoff_/IMS_cCut_);
      IMS_egCut_ = - eterm * (tmp);
      if (fabs(IMS_egCut_ - oldV) < 1.0E-14) {
        converged = true;
      } 
    }
    if (!converged) {
      throw CanteraError(" IdealMolalSoln::calcCutoffParams_()",
                         " failed to converge on the g polynomial");
    }
  }
 
}