MolalityVPSSTP.h

Go to the documentation of this file.

/**
 *  @file MolalityVPSSTP.h
 *   Header for intermediate ThermoPhase object for phases which
 *   employ molality based activity coefficient formulations
 *  (see \ref thermoprops
 * and class \link Cantera::MolalityVPSSTP MolalityVPSSTP\endlink).
 *
 * Header file for a derived class of ThermoPhase that handles
 * variable pressure standard state methods for calculating
 * thermodynamic properties that are further based upon activities
 * based on the molality scale.  These include most of the methods for
 * calculating liquid electrolyte thermodynamics.
 */
/*
 * Copyright (2006) Sandia Corporation. Under the terms of
 * Contract DE-AC04-94AL85000 with Sandia Corporation, the
 * U.S. Government retains certain rights in this software.
 */
#ifndef CT_MOLALITYVPSSTP_H
#define CT_MOLALITYVPSSTP_H

#include "VPStandardStateTP.h"

namespace Cantera
{

/**
 * @ingroup thermoprops
 */

/*!
 * MolalityVPSSTP is a derived class of ThermoPhase that handles
 * variable pressure standard state methods for calculating
 * thermodynamic properties that are further based on
 * molality-scaled activities.
 * This category incorporates most of the methods
 * for calculating liquid electrolyte thermodynamics that have been
 * developed since the 1970's.
 *
 * This class adds additional functions onto the %ThermoPhase interface
 * that handle molality based standard states. The %ThermoPhase
 * class includes a member function, ThermoPhase::activityConvention()
 * that indicates which convention the activities are based on. The
 * default is to assume activities are based on the molar convention.
 * However, classes which derive from the MolalityVPSSTP class return
 * <b>cAC_CONVENTION_MOLALITY</b> from this member function.
 *
 * The molality of a solute, \f$ m_i \f$, is defined as
 *
 * \f[
 *     m_i = \frac{n_i}{\tilde{M}_o n_o}
 * \f]
 * where
 * \f[
 *     \tilde{M}_o = \frac{M_o}{1000}
 * \f]
 *
 * where \f$ M_o \f$ is the molecular weight of the solvent. The molality
 * has units of gmol kg<SUP>-1</SUP>. For the solute, the molality may be
 * considered as the amount of gmol's of solute per kg of solvent, a natural
 * experimental quantity.
 *
 * The formulas for calculating mole fractions if given the molalities of
 * the solutes is stated below. First calculate \f$ L^{sum} \f$, an intermediate
 * quantity.
 *
 *   \f[
 *       L^{sum} = \frac{1}{\tilde{M}_o X_o} = \frac{1}{\tilde{M}_o} + \sum_{i\ne o} m_i
 *  \f]
 *  Then,
 *  \f[
 *       X_o =   \frac{1}{\tilde{M}_o L^{sum}}
 *  \f]
 *  \f[
 *       X_i =   \frac{m_i}{L^{sum}}
 *  \f]
 * where \f$ X_o \f$ is the mole fraction of solvent, and \f$ X_o \f$ is the
 * mole fraction of solute <I>i</I>. Thus, the molality scale and the mole fraction
 * scale offer a one-to-one mapping between each other, except in the limit
 * of a zero solvent mole fraction.
 *
 * The standard states for thermodynamic objects that derive from <b>MolalityVPSSTP</b>
 * are on the unit molality basis. Chemical potentials
 * of the solutes,  \f$ \mu_k \f$, and the solvent, \f$ \mu_o \f$, which are based
 * on the molality form, have the following general format:
 *
 * \f[
 *    \mu_k = \mu^{\triangle}_k(T,P) + R T ln(\gamma_k^{\triangle} \frac{m_k}{m^\triangle})
 * \f]
 * \f[
 *    \mu_o = \mu^o_o(T,P) + RT ln(a_o)
 * \f]
 *
 * where \f$ \gamma_k^{\triangle} \f$ is the molality based activity coefficient for species
 * \f$k\f$.
 *
 * The chemical potential of the solvent is thus expressed in a different format
 * than the chemical potential of the solutes. Additionally, the activity of the
 * solvent, \f$ a_o \f$, is further reexpressed in terms of an osmotic coefficient,
 * \f$ \phi \f$.
 *   \f[
 *       \phi = \frac{- ln(a_o)}{\tilde{M}_o \sum_{i \ne o} m_i}
 *   \f]
 *
 *  MolalityVPSSTP::osmoticCoefficient() returns the value of \f$ \phi \f$.
 *  Note there  are a few of definitions of the osmotic coefficient floating
 *  around. We use the one defined in
 *  (Activity Coefficients in Electrolyte Solutions, K. S. Pitzer
 *   CRC Press, Boca Raton, 1991, p. 85, Eqn. 28). This definition is most clearly
 *  related to theoretical calculation.
 *
 * The molar-based activity coefficients \f$ \gamma_k \f$ may be calculated
 * from the molality-based
 * activity coefficients, \f$ \gamma_k^\triangle \f$ by the following
 * formula.
 * \f[
 *     \gamma_k = \frac{\gamma_k^\triangle}{X_o}
 * \f]
 * For purposes of establishing a convention, the molar activity coefficient of the
 * solvent is set equal to the molality-based activity coefficient of the
 * solvent:
 * \f[
 *     \gamma_o = \gamma_o^\triangle
 * \f]
 *
 *  The molality-based and molarity-based standard states may be related to one
 *  another by the following formula.
 *
 * \f[
 *   \mu_k^\triangle(T,P) = \mu_k^o(T,P) + R T \ln(\tilde{M}_o m^\triangle)
 * \f]
 *
 *  An important convention is followed in all routines that derive from <b>%MolalityVPSSTP</b>.
 *  Standard state thermodynamic functions and reference state thermodynamic functions
 *  return the molality-based quantities. Also all functions which return
 *  activities return the molality-based activities. The reason for this convention
 *  has been discussed in supporting memos. However, it's important because the
 *  term in the equation above is non-trivial. For example it's equal
 *  to 2.38 kcal gmol<SUP>-1</SUP> for water at 298 K.
 *
 *
 * In order to prevent a singularity, this class includes the concept of a minimum
 * value for the solvent mole fraction. All calculations involving the formulation
 * of activity coefficients and other non-ideal solution behavior adhere to
 * this concept of a minimal value for the solvent mole fraction. This makes sense
 * because these solution behavior were all designed and measured far away from
 * the zero solvent singularity condition and are not applicable in that limit.
 *
 *
 * This objects add a layer that supports molality. It inherits from VPStandardStateTP.
 *
 *  All objects that derive from this are assumed to have molality based standard states.
 *
 *  Molality based activity coefficients are scaled according to the current
 *  pH scale. See the Eq3/6 manual for details.
 *
 *  Activity coefficients for species k may be altered between scales s1 to s2
 *  using the following formula
 *
 *   \f[
 *       ln(\gamma_k^{s2}) = ln(\gamma_k^{s1})
 *          + \frac{z_k}{z_j} \left(  ln(\gamma_j^{s2}) - ln(\gamma_j^{s1}) \right)
 *   \f]
 *
 *  where j is any one species. For the NBS scale, j is equal to the Cl- species
 *  and
 *
 *  \f[
 *       ln(\gamma_{Cl-}^{s2}) = \frac{-A_{\phi} \sqrt{I}}{1.0 + 1.5 \sqrt{I}}
 *  \f]
 *
 *  The Pitzer scale doesn't actually change anything. The pitzer scale is defined
 *  as the raw unscaled activity coefficients produced by the underlying objects.
 *
 *  <H3> SetState Strategy  </H3>
 *
 *   The MolalityVPSSTP object does not have a setState strategy concerning the
 *   molalities. It does not keep track of whether the molalities have changed.
 *   It's strictly an interfacial layer that writes the current mole fractions to the
 *   State object. When molalities are needed it recalculates the molalities from
 *   the State object's mole fraction vector.
 *
 *
 * @todo Make two solvent minimum fractions. One would be for calculation of the non-ideal
 *       factors. The other one would be for purposes of stoichiometry evaluation. the
 *       stoichiometry evaluation one would be a 1E-13 limit. Anything less would create
 *       problems with roundoff error.
 *
 */
class MolalityVPSSTP : public VPStandardStateTP
{

public:

    /// Constructors
    /*!
     * This doesn't do much more than initialize constants with
     * default values for water at 25C. Water molecular weight
     * comes from the default elements.xml file. It actually
     * differs slightly from the IAPWS95 value of 18.015268. However,
     * density conservation and therefore element conservation
     * is the more important principle to follow.
     */
    MolalityVPSSTP();

    //! Copy constructor
    /*!
     *  Note this stuff will not work until the underlying phase
     *  has a working copy constructor
     *
     * @param b class to be copied
     */
    MolalityVPSSTP(const MolalityVPSSTP& b);

    /// Assignment operator
    /*!
     *  Note this stuff will not work until the underlying phase
     *  has a working assignment operator
     *
     * @param b class to be copied.
     */
    MolalityVPSSTP& operator=(const MolalityVPSSTP& b);

    /// Destructor.
    virtual ~MolalityVPSSTP();

    //! Duplication routine for objects which inherit from  ThermoPhase.
    /*!
     *  This virtual routine can be used to duplicate thermophase objects
     *  inherited from ThermoPhase even if the application only has
     *  a pointer to ThermoPhase to work with.
     */
    virtual ThermoPhase* duplMyselfAsThermoPhase() const;

    /**
     *
     * @name  Utilities
     * @{
     */


    //! Equation of state type flag.
    /*!
     * The ThermoPhase base class returns
     * zero. Subclasses should define this to return a unique
     * non-zero value. Known constants defined for this purpose are
     * listed in mix_defs.h. The MolalityVPSSTP class also returns
     * zero, as it is a non-complete class.
     */
    virtual int eosType() const;

    //! Set the pH scale, which determines the scale for single-ion activity
    //! coefficients.
    /*!
     *  Single ion activity coefficients are not unique in terms of the
     *  representing actual measurable quantities.
     *
     * @param pHscaleType  Integer representing the pHscale
     */
    void setpHScale(const int pHscaleType);

    //! Reports the pH scale, which determines the scale for single-ion activity
    //! coefficients.
    /*!
     *  Single ion activity coefficients are not unique in terms of the
     *  representing actual measurable quantities.
     *
     * @return Return the pHscale type
     */
    int pHScale() const;

    /**
     * @}
     * @name  Molar Thermodynamic Properties
     * @{
     */


    /**
     * @}
     * @name Utilities for Solvent ID and Molality
     * @{
     */

    /**
     * This routine sets the index number of the solvent for
     * the phase.
     *
     *  Note, having a solvent
     *   is a precursor to many things having to do with molality.
     *
     * @param k the solvent index number
     */
    void setSolvent(size_t k);

    /**
     * Sets the minimum mole fraction in the molality formulation.
     * Note the molality formulation is singular in the limit that
     * the solvent mole fraction goes to zero. Numerically, how
     * this limit is treated and resolved is an ongoing issue within
     * Cantera.
     *
     * @param xmolSolventMIN  Input double containing the minimum mole fraction
     */
    void setMoleFSolventMin(doublereal xmolSolventMIN);

    //! Returns the solvent index.
    size_t solventIndex() const;

    /**
     * Returns the minimum mole fraction in the molality
     * formulation.
     */
    doublereal moleFSolventMin() const;

    //! Calculates the molality of all species and stores the result internally.
    /*!
     *   We calculate the vector of molalities of the species
     *   in the phase and store the result internally:
     *   \f[
     *     m_i = \frac{X_i}{1000 * M_o * X_{o,p}}
     *   \f]
     *    where
     *    - \f$ M_o \f$ is the molecular weight of the solvent
     *    - \f$ X_o \f$ is the mole fraction of the solvent
     *    - \f$ X_i \f$ is the mole fraction of the solute.
     *    - \f$ X_{o,p} = max (X_{o}^{min}, X_o) \f$
     *    - \f$ X_{o}^{min} \f$ = minimum mole fraction of solvent allowed
     *              in the denominator.
     */
    void calcMolalities() const;

    //!  This function will return the molalities of the species.
    /*!
     *   We calculate the vector of molalities of the species
     *   in the phase
     * \f[
     *     m_i = \frac{X_i}{1000 * M_o * X_{o,p}}
     * \f]
     *    where
     *    - \f$ M_o \f$ is the molecular weight of the solvent
     *    - \f$ X_o \f$ is the mole fraction of the solvent
     *    - \f$ X_i \f$ is the mole fraction of the solute.
     *    - \f$ X_{o,p} = \max (X_{o}^{min}, X_o) \f$
     *    - \f$ X_{o}^{min} \f$ = minimum mole fraction of solvent allowed
     *              in the denominator.
     *
     * @param molal       Output vector of molalities. Length: m_kk.
     */
    void getMolalities(doublereal* const molal) const;

    //! Set the molalities of the solutes in a phase
    /*!
     * Note, the entry for the solvent is not used.
     *   We are supplied with the molalities of all of the
     *   solute species. We then calculate the mole fractions of all
     *   species and update the %ThermoPhase object.
     *  \f[
     *     m_i = \frac{X_i}{M_o/1000 * X_{o,p}}
     *  \f]
     *    where
     *    -  \f$M_o\f$ is the molecular weight of the solvent
     *    -  \f$X_o\f$ is the mole fraction of the solvent
     *    -  \f$X_i\f$ is the mole fraction of the solute.
     *    -  \f$X_{o,p} = \max(X_o^{min}, X_o)\f$
     *    -  \f$X_o^{min}\f$ = minimum mole fraction of solvent allowed
     *                     in the denominator.
     *
     * The formulas for calculating mole fractions are
     *  \f[
     *       L^{sum} = \frac{1}{\tilde{M}_o X_o} = \frac{1}{\tilde{M}_o} + \sum_{i\ne o} m_i
     *  \f]
     *  Then,
     *  \f[
     *       X_o =   \frac{1}{\tilde{M}_o L^{sum}}
     *  \f]
     *  \f[
     *       X_i =   \frac{m_i}{L^{sum}}
     *  \f]
     *  It is currently an error if the solvent mole fraction is attempted to be set
     *  to a value lower than  \f$X_o^{min}\f$.
     *
     * @param molal   Input vector of molalities. Length: m_kk.
     */
    void setMolalities(const doublereal* const molal);

    //! Set the molalities of a phase
    /*!
     * Set the molalities of the solutes in a phase. Note, the entry for the
     * solvent is not used.
     *
     * @param xMap  Composition Map containing the molalities.
     */
    void setMolalitiesByName(compositionMap& xMap);

    //! Set the molalities of a phase
    /*!
     * Set the molalities of the solutes in a phase. Note, the entry for the
     * solvent is not used.
     *
     * @param name  String containing the information for a composition map.
     */
    void setMolalitiesByName(const std::string& name);

    /**
     * @}
     * @name Mechanical Properties
     * @{
     */

    /**
     * @}
     * @name Potential Energy
     *
     * Species may have an additional potential energy due to the
     * presence of external gravitation or electric fields. These
     * methods allow specifying a potential energy for individual
     * species.
     * @{
     */

    /**
     * @}
     * @name Activities, Standard States, and Activity Concentrations
     *
     * The activity \f$a_k\f$ of a species in solution is
     * related to the chemical potential by \f[ \mu_k = \mu_k^0(T)
     * + \hat R T \log a_k. \f] The quantity \f$\mu_k^0(T,P)\f$ is
     * the chemical potential at unit activity, which depends only
     * on temperature and pressure.
     * @{
     */

    /**
     * This method returns the activity convention.
     * Currently, there are two activity conventions
     *  Molar-based activities
     *       Unit activity of species at either a hypothetical pure
     *       solution of the species or at a hypothetical
     *       pure ideal solution at infinite dilution
     *   cAC_CONVENTION_MOLAR 0
     *      - default
     *
     *  Molality based activities
     *       (unit activity of solutes at a hypothetical 1 molal
     *        solution referenced to infinite dilution at all
     *        pressures and temperatures).
     *   cAC_CONVENTION_MOLALITY 1
     *
     *  We set the convention to molality here.
     */
    int activityConvention() const;

    /**
     * This method returns an array of generalized concentrations
     * \f$ C_k\f$ that are defined such that
     * \f$ a_k = C_k / C^0_k, \f$ where \f$ C^0_k \f$
     * is a standard concentration
     * defined below.  These generalized concentrations are used
     * by kinetics manager classes to compute the forward and
     * reverse rates of elementary reactions.
     *
     * @param c Array of generalized concentrations. The
     *           units depend upon the implementation of the
     *           reaction rate expressions within the phase.
     */
    virtual void getActivityConcentrations(doublereal* c) const;

    /**
     * The standard concentration \f$ C^0_k \f$ used to normalize
     * the generalized concentration. In many cases, this quantity
     * will be the same for all species in a phase - for example,
     * for an ideal gas \f$ C^0_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.
     *
     * @param k species index. Defaults to zero.
     */
    virtual doublereal standardConcentration(size_t k=0) const;

    /**
     * Returns the natural logarithm of the standard
     * concentration of the kth species
     *
     * @param k  species index
     */
    virtual doublereal logStandardConc(size_t k=0) const;

    /**
     * Returns the units of the standard and generalized
     * concentrations Note they have the same units, as their
     * ratio 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.
     *
     * @param uA Output vector containing the units
     *  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
     * @param k species index. Defaults to 0.
     * @param sizeUA output int containing the size of the vector.
     *        Currently, this is equal to 6.
     */
    virtual void getUnitsStandardConc(double* uA, int k = 0,
                                      int sizeUA = 6) const;


    //! Get the array of non-dimensional activities (molality
    //! based for this class and classes that derive from it) at
    //! the current solution temperature, pressure, and solution concentration.
    /*!
     * All standard state properties for molality-based phases are
     * evaluated consistent with the molality scale. Therefore, this function
     * must return molality-based activities.
     *
     * \f[
     *  a_i^\triangle = \gamma_k^{\triangle} \frac{m_k}{m^\triangle}
     * \f]
     *
     * This function must be implemented in derived classes.
     *
     * @param ac     Output vector of molality-based activities. Length: m_kk.
     */
    virtual void getActivities(doublereal* ac) const;

    //! Get the array of non-dimensional activity coefficients at
    //! the current solution temperature, pressure, and solution concentration.
    /*!
     * These are mole-fraction based activity coefficients. In this
     * object, their calculation is based on translating the values
     * of the molality-based activity coefficients.
     *  See Denbigh p. 278 for a thorough discussion.
     *
     * The molar-based activity coefficients \f$ \gamma_k \f$ may be calculated from the
     * molality-based
     * activity coefficients, \f$ \gamma_k^\triangle \f$ by the following
     * formula.
     * \f[
     *     \gamma_k = \frac{\gamma_k^\triangle}{X_o}
     * \f]
     *
     * For purposes of establishing a convention, the molar activity coefficient of the
     * solvent is set equal to the molality-based activity coefficient of the
     * solvent:
     *
     * \f[
     *     \gamma_o = \gamma_o^\triangle
     * \f]
     *
     * Derived classes don't need to overload this function. This function is
     * handled at this level.
     *
     * @param ac  Output vector containing the mole-fraction based activity coefficients.
     *            length: m_kk.
     */
    void getActivityCoefficients(doublereal* ac) const;

    //!  Get the array of non-dimensional molality based
    //!  activity coefficients at the current solution temperature,
    //!  pressure, and  solution concentration.
    /*!
     *  See Denbigh p. 278 for a thorough discussion. This class must be overwritten in
     *  classes which derive from %MolalityVPSSTP. This function takes over from the
     *  molar-based activity coefficient calculation, getActivityCoefficients(), in
     *  derived classes.
     *
     *  These molality based activity coefficients are scaled according to the current
     *  pH scale. See the Eq3/6 manual for details.
     *
     *  Activity coefficients for species k may be altered between scales s1 to s2
     *  using the following formula
     *
     *   \f[
     *       ln(\gamma_k^{s2}) = ln(\gamma_k^{s1})
     *          + \frac{z_k}{z_j} \left(  ln(\gamma_j^{s2}) - ln(\gamma_j^{s1}) \right)
     *   \f]
     *
     *  where j is any one species. For the NBS scale, j is equal to the Cl- species
     *  and
     *
     *  \f[
     *       ln(\gamma_{Cl-}^{s2}) = \frac{-A_{\phi} \sqrt{I}}{1.0 + 1.5 \sqrt{I}}
     *  \f]
     *
     * @param acMolality Output vector containing the molality based activity coefficients.
     *                   length: m_kk.
     */
    virtual void getMolalityActivityCoefficients(doublereal* acMolality) const;



    //! Calculate the osmotic coefficient
    /*!
     *   \f[
     *       \phi = \frac{- ln(a_o)}{\tilde{M}_o \sum_{i \ne o} m_i}
     *   \f]
     *
     *  Note there  are a few of definitions of the osmotic coefficient floating
     *  around. We use the one defined in
     *  (Activity Coefficients in Electrolyte Solutions, K. S. Pitzer
     *   CRC Press, Boca Raton, 1991, p. 85, Eqn. 28). This definition is most clearly
     *  related to theoretical calculation.
     *
     *  units = dimensionless
     */
    virtual double osmoticCoefficient() const;

    //@}
    /// @name  Partial Molar Properties of the Solution
    //@{


    /**
     * Get the species electrochemical potentials.
     * These are partial molar quantities.
     * This method adds a term \f$ Fz_k \phi_k \f$ to the
     * to each chemical potential.
     *
     * Units: J/kmol
     *
     * @param mu     output vector containing the species electrochemical potentials.
     *               Length: m_kk.
     */
    void getElectrochemPotentials(doublereal* mu) const;


    //@}
    /// @name  Properties of the Standard State of the Species in the Solution
    //@{



    //@}
    /// @name Thermodynamic Values for the Species Reference States
    //@{


    ///////////////////////////////////////////////////////
    //
    //  The methods below are not virtual, and should not
    //  be overloaded.
    //
    //////////////////////////////////////////////////////

    /**
     * @name Specific Properties
     * @{
     */


    /**
     * @name Setting the State
     *
     * These methods set all or part of the thermodynamic
     * state.
     * @{
     */

    //@}

    /**
     * @name Chemical Equilibrium
     * Routines that implement the Chemical equilibrium capability
     * for a single phase, based on the element-potential method.
     * @{
     */

    /**
     * This method is used by the ChemEquil element-potential
     * based equilibrium solver.
     * It sets the state such that the chemical potentials of the
     * species within the current phase satisfy
     * \f[ \frac{\mu_k}{\hat R T} = \sum_m A_{k,m}
     * \left(\frac{\lambda_m} {\hat R T}\right) \f] where
     * \f$ \lambda_m \f$ is the element potential of element m. The
     * temperature is unchanged.  Any phase (ideal or not) that
     * implements this method can be equilibrated by ChemEquil.
     *
     * @param lambda_RT Input vector containing the dimensionless
     *                  element potentials.
     */
    virtual void setToEquilState(const doublereal* lambda_RT);


    //@}


    //! 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.
     *
     * The MolalityVPSSTP object defines a new method for setting
     * the concentrations of a phase. The new method is defined by a
     * block called "soluteMolalities". If this block
     * is found, the concentrations within that phase are
     * set to the "name":"molalities pairs found within that
     * XML block. The solvent concentration is then set
     * to everything else.
     *
     * The function first calls the overloaded function ,
     * VPStandardStateTP::setStateFromXML(), to pick up the parent class
     * behavior.
     *
     * usage: Overloaded functions should call this function
     *        before carrying out their own behavior.
     *
     * @param state An XML_Node object corresponding to
     *              the "state" entry for this phase in the input file.
     */
    virtual void setStateFromXML(const XML_Node& state);

    /// The following methods are used in the process of constructing
    /// the phase and setting its parameters from a specification in an
    /// input file. They are not normally used in application programs.
    /// To see how they are used, see files importCTML.cpp and
    /// ThermoFactory.cpp.


    /*!
     * @internal Initialize. This method is provided to allow
     * subclasses to perform any initialization required after all
     * species have been added. For example, it might be used to
     * resize internal work arrays that must have an entry for
     * each species.  The base class implementation does nothing,
     * and subclasses that do not require initialization do not
     * need to overload this method.  When importing a CTML phase
     * description, this method is called just prior to returning
     * from function importPhase.
     *
     * @see importCTML.cpp
     */
    virtual void initThermo();


    /**
     *   Import and initialize a ThermoPhase object
     *
     * @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 initThermoXML(XML_Node& phaseNode, std::string id);


    //! Set the temperature (K), pressure (Pa), and molalities
    //!(gmol kg-1) of the solutes
    /*!
     * @param t           Temperature (K)
     * @param p           Pressure (Pa)
     * @param molalities  Input vector of molalities of the solutes.
     *                    Length: m_kk.
     */
    void setState_TPM(doublereal t, doublereal p,
                      const doublereal* const molalities);

    //! Set the temperature (K), pressure (Pa), and molalities.
    /*!
     * @param t           Temperature (K)
     * @param p           Pressure (Pa)
     * @param m           compositionMap containing the molalities
     */
    void setState_TPM(doublereal t, doublereal p, compositionMap& m);

    //! Set the temperature (K), pressure (Pa), and molalities.
    /*!
     * @param t           Temperature (K)
     * @param p           Pressure (Pa)
     * @param m           String which gets translated into a composition
     *                    map for the molalities of the solutes.
     */
    void setState_TPM(doublereal t, doublereal p, const std::string& m);

    //! Get the array of derivatives of the log activity coefficients with respect to the log of the species mole numbers
    /*!
     * Implementations should take the derivative of the logarithm of the activity coefficient with respect to a
     * species log mole number (with all other species mole numbers held constant). The default treatment in the
     * %ThermoPhase object is to set this vector to zero.
     *
     *  units = 1 / kmol
     *
     *  dlnActCoeffdlnN[ ld * k  + m]  will contain the derivative of log act_coeff for the <I>m</I><SUP>th</SUP>
     *                               species with respect to the number of moles of the <I>k</I><SUP>th</SUP> species.
     *
     * \f[
     *        \frac{d \ln(\gamma_m) }{d \ln( n_k ) }\Bigg|_{n_i}
     * \f]
     *
     * @param ld               Number of rows in the matrix
     * @param dlnActCoeffdlnN    Output vector of derivatives of the
     *                           log Activity Coefficients. length = m_kk * m_kk
     */
    virtual void getdlnActCoeffdlnN(const size_t ld, doublereal* const dlnActCoeffdlnN) {
        getdlnActCoeffdlnN_numderiv(ld, dlnActCoeffdlnN);
    }

    //! returns a summary of the state of the phase as a string
    /*!
     * @param show_thermo If true, extra information is printed out
     *                    about the thermodynamic state of the system.
     */
    virtual std::string report(bool show_thermo = true) const;

    //! returns a summary of the state of the phase to specified
    //! comma separated files
    /*!
     * @param csvFile     ofstream file to print comma separated data for
     *                    the phase
     */
    virtual void reportCSV(std::ofstream& csvFile) const;

protected:

    //! Get the array of unscaled non-dimensional molality based
    //!  activity coefficients at the current solution temperature,
    //!  pressure, and  solution concentration.
    /*!
     *  See Denbigh p. 278 for a thorough discussion. This class must be overwritten in
     *  classes which derive from %MolalityVPSSTP. This function takes over from the
     *  molar-based activity coefficient calculation, getActivityCoefficients(), in
     *  derived classes.
     *
     * @param acMolality Output vector containing the molality based activity coefficients.
     *                   length: m_kk.
     */
    virtual void getUnscaledMolalityActivityCoefficients(doublereal* acMolality) const;

    //! Apply the current phScale to a set of activity Coefficients or activities
    /*!
     *  See the Eq3/6 Manual for a thorough discussion.
     *
     * @param acMolality input/Output vector containing the molality based
     *                   activity coefficients. length: m_kk.
     */
    virtual void applyphScale(doublereal* acMolality) const;

private:
    //! Returns the index of the Cl- species.
    /*!
     *  The Cl- species is special in the sense that its single ion
     *  molality-based activity coefficient is used in the specification
     *  of the pH scale for single ions. Therefore, we need to know
     *  what species index is Cl-. If the species isn't in the species
     *  list then this routine returns -1, and we can't use the NBS
     *  pH scale.
     *
     *  Right now we use a restrictive interpretation. The species
     *  must be named "Cl-". It must consist of exactly one Cl and one E
     *  atom.
     */
    virtual size_t findCLMIndex() const;

    //! Initialize lengths of local variables after all species have
    //! been identified.
    void initLengths();

protected:

    //! Index of the solvent
    /*!
     * Currently the index of the solvent is hard-coded to the value 0
     */
    size_t m_indexSolvent;

    //! Scaling to be used for output of single-ion species activity
    //! coefficients.
    /*!
     *   Index of the species to be used in the single-ion scaling
     *   law. This is the identity of the Cl- species for the PHSCALE_NBS
     *   scaling.
     *   Either PHSCALE_PITZER or PHSCALE_NBS
     */
    int        m_pHScalingType;

    //! Index of the phScale species
    /*!
     *   Index of the species to be used in the single-ion scaling
     *   law. This is the identity of the Cl- species for the PHSCALE_NBS
     *   scaling
     */
    size_t m_indexCLM;

    //! Molecular weight of the Solvent
    doublereal m_weightSolvent;

    /*!
     * In any molality implementation, it makes sense to have
     * a minimum solvent mole fraction requirement, since the
     * implementation becomes singular in the xmolSolvent=0
     * limit. The default is to set it to 0.01.
     * We then modify the molality definition to ensure that
     * molal_solvent = 0 when xmol_solvent = 0.
     */
    doublereal m_xmolSolventMIN;

    //! This is the multiplication factor that goes inside
    //! log expressions involving the molalities of species.
    /*!
     * It's equal to Wt_0 / 1000,
     *     where Wt_0 = weight of solvent (kg/kmol)
     */
    doublereal m_Mnaught;

    //! Current value of the molalities of the species in the phase.
    /*!
     * Note this vector is a mutable quantity.
     * units are (kg/kmol)
     */
    mutable vector_fp  m_molalities;

private:
    //! Error function
    /*!
     *  Print an error string and exit
     *
     * @param msg  Message to be printed
     */
    doublereal err(std::string msg) const;

};


//! Scale to be used for the output of single-ion activity coefficients
//! is that used by Pitzer.
/*!
 *  This is the internal scale used within the code. One property is that
 *  the activity coefficients for the cation and anion of a single salt
 *  will be equal. This scale is the one presumed by the formulation of the
 *  single-ion activity coefficients described in this report.
 *
 *  Activity coefficients for species k may be altered between scales s1 to s2
 *  using the following formula
 *
 *   \f[
 *       ln(\gamma_k^{s2}) = ln(\gamma_k^{s1})
 *          + \frac{z_k}{z_j} \left(  ln(\gamma_j^{s2}) - ln(\gamma_j^{s1}) \right)
 *   \f]
 *
 *  where j is any one species.
 *
 *
 */
const int PHSCALE_PITZER = 0;

//! Scale to be used for evaluation of single-ion activity coefficients
//! is that used by the NBS standard for evaluation of the pH variable.
/*!
 *  This is not the internal scale used within the code.
 *
 *  Activity coefficients for species k may be altered between scales s1 to s2
 *  using the following formula
 *
 *   \f[
 *       ln(\gamma_k^{s2}) = ln(\gamma_k^{s1})
 *          + \frac{z_k}{z_j} \left(  ln(\gamma_j^{s2}) - ln(\gamma_j^{s1}) \right)
 *   \f]
 *
 *  where j is any one species. For the NBS scale, j is equal to the Cl- species
 *  and
 *
 *  \f[
 *       ln(\gamma_{Cl-}^{s2}) = \frac{-A_{\phi} \sqrt{I}}{1.0 + 1.5 \sqrt{I}}
 *  \f]
 *
 *  This is the NBS pH scale, which is used in all conventional pH
 *  measurements. and is based on the Bates-Guggenheim equations.
 *
 */
const int PHSCALE_NBS    = 1;

}

#endif