IonsFromNeutralVPSSTP.cpp

Go to the documentation of this file.

/**
 *  @file IonsFromNeutralVPSSTP.cpp
 *   Definitions for the object which treats ionic liquids as made of ions as species
 *   even though the thermodynamics is obtained from the neutral molecule representation.
 *  (see \ref thermoprops
 *   and class \link Cantera::IonsFromNeutralVPSSTP IonsFromNeutralVPSSTP\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 expressions
 * for the excess gibbs free energy expressed as a function of
 * the mole fractions.
 */
/*
 * Copyright (2009) Sandia Corporation. Under the terms of
 * Contract DE-AC04-94AL85000 with Sandia Corporation, the
 * U.S. Government retains certain rights in this software.
 */
#include "cantera/thermo/IonsFromNeutralVPSSTP.h"
#include "cantera/thermo/ThermoFactory.h"
#include "cantera/thermo/PDSS_IonsFromNeutral.h"
#include "cantera/thermo/mix_defs.h"
#include "cantera/base/stringUtils.h"

#include <cmath>
#include <iomanip>
#include <fstream>

using namespace std;

namespace Cantera
{

static  const double xxSmall = 1.0E-150;
//====================================================================================================================
/*
 * Default constructor.
 *
 */
IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP() :
    GibbsExcessVPSSTP(),
    ionSolnType_(cIonSolnType_SINGLEANION),
    numNeutralMoleculeSpecies_(0),
    indexSpecialSpecies_(npos),
    indexSecondSpecialSpecies_(npos),
    numCationSpecies_(0),
    numAnionSpecies_(0),
    numPassThroughSpecies_(0),
    neutralMoleculePhase_(0),
    IOwnNThermoPhase_(true),
    moleFractionsTmp_(0),
    muNeutralMolecule_(0),
    lnActCoeff_NeutralMolecule_(0)
{
}

//====================================================================================================================
// Construct and initialize an IonsFromNeutralVPSSTP object
// directly from an ASCII input file
/*
 * Working constructors
 *
 *  The two constructors below are the normal way
 *  the phase initializes itself. They are shells that call
 *  the routine initThermo(), with a reference to the
 *  XML database to get the info for the phase.
 *
 * @param inputFile Name of the input file containing the phase XML data
 *                  to set up the object
 * @param id        ID of the phase in the input file. Defaults to the
 *                  empty string.
 * @param neutralPhase   The object takes a neutralPhase ThermoPhase
 *                       object as input. It can either take a pointer
 *                       to an existing object in the parameter list,
 *                       in which case it does not own the object, or
 *                       it can construct a neutral Phase as a slave
 *                       object, in which case, it does own the slave
 *                       object, for purposes of who gets to destroy
 *                       the object.
 *                       If this parameter is zero, then a slave
 *                       neutral phase object is created and used.
 */
IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP(std::string inputFile, std::string id,
        ThermoPhase* neutralPhase) :
    GibbsExcessVPSSTP(),
    ionSolnType_(cIonSolnType_SINGLEANION),
    numNeutralMoleculeSpecies_(0),
    indexSpecialSpecies_(npos),
    indexSecondSpecialSpecies_(npos),
    numCationSpecies_(0),
    numAnionSpecies_(0),
    numPassThroughSpecies_(0),
    neutralMoleculePhase_(neutralPhase),
    IOwnNThermoPhase_(true),
    moleFractionsTmp_(0),
    muNeutralMolecule_(0),
    lnActCoeff_NeutralMolecule_(0)
{
    if (neutralPhase) {
        IOwnNThermoPhase_ = false;
    }
    constructPhaseFile(inputFile, id);
}
//====================================================================================================================
IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP(XML_Node& phaseRoot, std::string id,
        ThermoPhase* neutralPhase) :
    GibbsExcessVPSSTP(),
    ionSolnType_(cIonSolnType_SINGLEANION),
    numNeutralMoleculeSpecies_(0),
    indexSpecialSpecies_(npos),
    indexSecondSpecialSpecies_(npos),
    numCationSpecies_(0),
    numAnionSpecies_(0),
    numPassThroughSpecies_(0),
    neutralMoleculePhase_(neutralPhase),
    IOwnNThermoPhase_(true),
    moleFractionsTmp_(0),
    muNeutralMolecule_(0),

    lnActCoeff_NeutralMolecule_(0)
{
    if (neutralPhase) {
        IOwnNThermoPhase_ = false;
    }
    constructPhaseXML(phaseRoot, id);
}

//====================================================================================================================

/*
 * Copy Constructor:
 *
 *  Note this stuff will not work until the underlying phase
 *  has a working copy constructor
 */
IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP(const IonsFromNeutralVPSSTP& b) :
    GibbsExcessVPSSTP(),
    ionSolnType_(cIonSolnType_SINGLEANION),
    numNeutralMoleculeSpecies_(0),
    indexSpecialSpecies_(npos),
    indexSecondSpecialSpecies_(npos),
    numCationSpecies_(0),
    numAnionSpecies_(0),
    numPassThroughSpecies_(0),
    neutralMoleculePhase_(0),
    IOwnNThermoPhase_(true),
    moleFractionsTmp_(0),
    muNeutralMolecule_(0),

    lnActCoeff_NeutralMolecule_(0)
{
    IonsFromNeutralVPSSTP::operator=(b);
}
//====================================================================================================================
/*
 * operator=()
 *
 *  Note this stuff will not work until the underlying phase
 *  has a working assignment operator
 */
IonsFromNeutralVPSSTP& IonsFromNeutralVPSSTP::
operator=(const IonsFromNeutralVPSSTP& b)
{
    if (&b == this) {
        return *this;
    }

    /*
     *  If we own the underlying neutral molecule phase, then we do a deep
     *  copy. If not, we do a shallow copy. We get a valid pointer for
     *  neutralMoleculePhase_ first, because we need it to assign the pointers
     *  within the PDSS_IonsFromNeutral object. which is done in the
     *  GibbsExcessVPSSTP::operator=(b) step.
     */
    if (IOwnNThermoPhase_) {
        if (b.neutralMoleculePhase_) {
            if (neutralMoleculePhase_) {
                delete neutralMoleculePhase_;
            }
            neutralMoleculePhase_   = (b.neutralMoleculePhase_)->duplMyselfAsThermoPhase();
        } else {
            neutralMoleculePhase_   = 0;
        }
    } else {
        neutralMoleculePhase_     = b.neutralMoleculePhase_;
    }

    GibbsExcessVPSSTP::operator=(b);

    ionSolnType_                = b.ionSolnType_;
    numNeutralMoleculeSpecies_  = b.numNeutralMoleculeSpecies_;
    indexSpecialSpecies_        = b.indexSpecialSpecies_;
    indexSecondSpecialSpecies_  = b.indexSecondSpecialSpecies_;
    fm_neutralMolec_ions_       = b.fm_neutralMolec_ions_;
    fm_invert_ionForNeutral     = b.fm_invert_ionForNeutral;
    NeutralMolecMoleFractions_  = b.NeutralMolecMoleFractions_;
    cationList_                 = b.cationList_;
    numCationSpecies_           = b.numCationSpecies_;
    anionList_                  = b.anionList_;
    numAnionSpecies_            = b.numAnionSpecies_;
    passThroughList_            = b.passThroughList_;
    numPassThroughSpecies_      = b.numPassThroughSpecies_;

    IOwnNThermoPhase_           = b.IOwnNThermoPhase_;
    moleFractionsTmp_           = b.moleFractionsTmp_;
    muNeutralMolecule_          = b.muNeutralMolecule_;
    //  gammaNeutralMolecule_       = b.gammaNeutralMolecule_;
    lnActCoeff_NeutralMolecule_ = b.lnActCoeff_NeutralMolecule_;
    dlnActCoeffdT_NeutralMolecule_ = b.dlnActCoeffdT_NeutralMolecule_;
    dlnActCoeffdlnX_diag_NeutralMolecule_ = b.dlnActCoeffdlnX_diag_NeutralMolecule_;
    dlnActCoeffdlnN_diag_NeutralMolecule_ = b.dlnActCoeffdlnN_diag_NeutralMolecule_;
    dlnActCoeffdlnN_NeutralMolecule_ = b.dlnActCoeffdlnN_NeutralMolecule_;

    return *this;
}

/*
 *
 * ~IonsFromNeutralVPSSTP():   (virtual)
 *
 * Destructor: does nothing:
 *
 */
IonsFromNeutralVPSSTP::~IonsFromNeutralVPSSTP()
{
    if (IOwnNThermoPhase_) {
        delete neutralMoleculePhase_;
        neutralMoleculePhase_ = 0;
    }
}

/*
 * This routine duplicates the current object and returns
 * a pointer to ThermoPhase.
 */
ThermoPhase*
IonsFromNeutralVPSSTP::duplMyselfAsThermoPhase() const
{
    IonsFromNeutralVPSSTP* mtp = new IonsFromNeutralVPSSTP(*this);
    return (ThermoPhase*) mtp;
}

/*
 *  -------------- 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 IonsFromNeutralVPSSTP class also returns
 * zero, as it is a non-complete class.
 */
int IonsFromNeutralVPSSTP::eosType() const
{
    return cIonsFromNeutral;
}



/*
 * ------------ Molar Thermodynamic Properties ----------------------
 */
/*
 * Molar enthalpy of the solution. Units: J/kmol.
 */
doublereal IonsFromNeutralVPSSTP::enthalpy_mole() const
{
    getPartialMolarEnthalpies(DATA_PTR(m_pp));
    return mean_X(DATA_PTR(m_pp));
}

/**
 * Molar internal energy of the solution. Units: J/kmol.
 *
 * This is calculated from the soln enthalpy and then
 * subtracting pV.
 */
doublereal IonsFromNeutralVPSSTP::intEnergy_mole() const
{
    double hh = enthalpy_mole();
    double pres = pressure();
    double molarV = 1.0/molarDensity();
    double uu = hh - pres * molarV;
    return uu;
}

/**
 *  Molar soln entropy at constant pressure. Units: J/kmol/K.
 *
 *  This is calculated from the partial molar entropies.
 */
doublereal IonsFromNeutralVPSSTP::entropy_mole() const
{
    getPartialMolarEntropies(DATA_PTR(m_pp));
    return mean_X(DATA_PTR(m_pp));
}

/// Molar Gibbs function. Units: J/kmol.
doublereal IonsFromNeutralVPSSTP::gibbs_mole() const
{
    getChemPotentials(DATA_PTR(m_pp));
    return mean_X(DATA_PTR(m_pp));
}
/** Molar heat capacity at constant pressure. Units: J/kmol/K.
  *
  * Returns the solution heat capacition at constant pressure.
  * This is calculated from the partial molar heat capacities.
  */
doublereal IonsFromNeutralVPSSTP::cp_mole() const
{
    getPartialMolarCp(DATA_PTR(m_pp));
    double val = mean_X(DATA_PTR(m_pp));
    return val;
}

/// Molar heat capacity at constant volume. Units: J/kmol/K.
doublereal IonsFromNeutralVPSSTP::cv_mole() const
{
    // Need to revisit this, as it is wrong
    getPartialMolarCp(DATA_PTR(m_pp));
    return mean_X(DATA_PTR(m_pp));
    //err("not implemented");
    //return 0.0;
}
//===========================================================================================================
/*
 * - Activities, Standard States, Activity Concentrations -----------
 */
//===========================================================================================================
void IonsFromNeutralVPSSTP::getDissociationCoeffs(vector_fp& coeffs,
        vector_fp& charges, std::vector<size_t>& neutMolIndex) const
{
    coeffs = fm_neutralMolec_ions_;
    charges = m_speciesCharge;
    neutMolIndex = fm_invert_ionForNeutral;
}
//===========================================================================================================
// Get the array of non-dimensional molar-based activity coefficients at
// the current solution temperature, pressure, and solution concentration.
/*
 * @param ac Output vector of activity coefficients. Length: m_kk.
 */
void IonsFromNeutralVPSSTP::getActivityCoefficients(doublereal* ac) const
{

    // This stuff has moved to the setState routines
    //   calcNeutralMoleculeMoleFractions();
    //   neutralMoleculePhase_->setState_TPX(temperature(), pressure(), DATA_PTR(NeutralMolecMoleFractions_));
    //   neutralMoleculePhase_->getStandardChemPotentials(DATA_PTR(muNeutralMolecule_));

    /*
     * Update the activity coefficients
     */
    s_update_lnActCoeff();

    /*
     * take the exp of the internally stored coefficients.
     */
    for (size_t k = 0; k < m_kk; k++) {
        ac[k] = exp(lnActCoeff_Scaled_[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 at the current temperature, pressure
 * and mole fraction of the solution.
 *
 * @param mu  Output vector of species chemical
 *            potentials. Length: m_kk. Units: J/kmol
 */
void
IonsFromNeutralVPSSTP::getChemPotentials(doublereal* mu) const
{
    size_t icat, jNeut;
    doublereal xx, fact2;
    /*
     *  Transfer the mole fractions to the slave neutral molecule
     *  phase
     *   Note we may move this in the future.
     */
    //calcNeutralMoleculeMoleFractions();
    //neutralMoleculePhase_->setState_TPX(temperature(), pressure(), DATA_PTR(NeutralMolecMoleFractions_));

    /*
     * Get the standard chemical potentials of netural molecules
     */
    neutralMoleculePhase_->getStandardChemPotentials(DATA_PTR(muNeutralMolecule_));

    doublereal RT_ = GasConstant * temperature();

    switch (ionSolnType_) {
    case cIonSolnType_PASSTHROUGH:
        neutralMoleculePhase_->getChemPotentials(mu);
        break;
    case cIonSolnType_SINGLEANION:
        //  neutralMoleculePhase_->getActivityCoefficients(DATA_PTR(gammaNeutralMolecule_));
        neutralMoleculePhase_->getLnActivityCoefficients(DATA_PTR(lnActCoeff_NeutralMolecule_));

        fact2 = 2.0 * RT_ * log(2.0);

        // Do the cation list
        for (size_t k = 0; k < cationList_.size(); k++) {
            //! Get the id for the next cation
            icat = cationList_[k];
            jNeut = fm_invert_ionForNeutral[icat];
            xx = std::max(SmallNumber, moleFractions_[icat]);
            mu[icat] = muNeutralMolecule_[jNeut] + fact2 + RT_ * (lnActCoeff_NeutralMolecule_[jNeut] + log(xx));
        }

        // Do the anion list
        icat = anionList_[0];
        jNeut = fm_invert_ionForNeutral[icat];
        xx = std::max(SmallNumber, moleFractions_[icat]);
        mu[icat] = RT_ * log(xx);

        // Do the list of neutral molecules
        for (size_t k = 0; k < numPassThroughSpecies_; k++) {
            icat = passThroughList_[k];
            jNeut = fm_invert_ionForNeutral[icat];
            xx = std::max(SmallNumber, moleFractions_[icat]);
            mu[icat] =  muNeutralMolecule_[jNeut]  + RT_ * (lnActCoeff_NeutralMolecule_[jNeut] + log(xx));
        }
        break;

    case cIonSolnType_SINGLECATION:
        throw CanteraError("eosType", "Unknown type");
        break;
    case cIonSolnType_MULTICATIONANION:
        throw CanteraError("eosType", "Unknown type");
        break;
    default:
        throw CanteraError("eosType", "Unknown type");
        break;
    }
}


// Returns an array of partial molar enthalpies for the species
// in the mixture.
/*
 * Units (J/kmol)
 *
 * For this phase, the partial molar enthalpies are equal to the
 * standard state enthalpies modified by the derivative of the
 * molality-based activity coefficient wrt temperature
 *
 *  \f[
 * \bar h_k(T,P) = h^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT}
 * \f]
 *
 */
void IonsFromNeutralVPSSTP::getPartialMolarEnthalpies(doublereal* hbar) const
{
    /*
      * Get the nondimensional standard state enthalpies
      */
    getEnthalpy_RT(hbar);
    /*
     * dimensionalize it.
     */
    double T = temperature();
    double RT = GasConstant * T;
    for (size_t k = 0; k < m_kk; k++) {
        hbar[k] *= RT;
    }
    /*
     * Update the activity coefficients, This also update the
     * internally stored molalities.
     */
    s_update_lnActCoeff();
    s_update_dlnActCoeffdT();
    double RTT = RT * T;
    for (size_t k = 0; k < m_kk; k++) {
        hbar[k] -= RTT * dlnActCoeffdT_Scaled_[k];
    }
}

// Returns an array of partial molar entropies for the species
// in the mixture.
/*
 * Units (J/kmol)
 *
 * For this phase, the partial molar enthalpies are equal to the
 * standard state enthalpies modified by the derivative of the
 * activity coefficient wrt temperature
 *
 *  \f[
 * \bar s_k(T,P) = s^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT}
 * \f]
 *
 */
void IonsFromNeutralVPSSTP::getPartialMolarEntropies(doublereal* sbar) const
{
    double xx;
    /*
     * Get the nondimensional standard state entropies
     */
    getEntropy_R(sbar);
    double T = temperature();
    /*
     * Update the activity coefficients, This also update the
     * internally stored molalities.
     */
    s_update_lnActCoeff();
    s_update_dlnActCoeffdT();

    for (size_t k = 0; k < m_kk; k++) {
        xx = std::max(moleFractions_[k], xxSmall);
        sbar[k] += - lnActCoeff_Scaled_[k] -log(xx) - T * dlnActCoeffdT_Scaled_[k];
    }
    /*
     * dimensionalize it.
     */
    for (size_t k = 0; k < m_kk; k++) {
        sbar[k] *= GasConstant;
    }
}


// Get the array of log concentration-like derivatives of the
// log activity coefficients
/*
 * This function is a virtual method.  For ideal mixtures
 * (unity activity coefficients), this can return zero.
 * Implementations should take the derivative of the
 * logarithm of the activity coefficient with respect to the
 * logarithm of the concentration-like variable (i.e. mole fraction,
 * molality, etc.) that represents the standard state.
 * This quantity is to be used in conjunction with derivatives of
 * that concentration-like variable when the derivative of the chemical
 * potential is taken.
 *
 *  units = dimensionless
 *
 * @param dlnActCoeffdlnX    Output vector of log(mole fraction)
 *                 derivatives of the log Activity Coefficients.
 *                 length = m_kk
 */
void IonsFromNeutralVPSSTP::getdlnActCoeffdlnX_diag(doublereal* dlnActCoeffdlnX_diag) const
{
    s_update_lnActCoeff();
    s_update_dlnActCoeff_dlnX_diag();

    for (size_t k = 0; k < m_kk; k++) {
        dlnActCoeffdlnX_diag[k] = dlnActCoeffdlnX_diag_[k];
    }
}
//====================================================================================================================
// Get the array of log concentration-like derivatives of the
// log activity coefficients
/*
 * This function is a virtual method.  For ideal mixtures
 * (unity activity coefficients), this can return zero.
 * Implementations should take the derivative of the
 * logarithm of the activity coefficient with respect to the
 * logarithm of the concentration-like variable (i.e. moles)
 * that represents the standard state.
 * This quantity is to be used in conjunction with derivatives of
 * that concentration-like variable when the derivative of the chemical
 * potential is taken.
 *
 *  units = dimensionless
 *
 * @param dlnActCoeffdlnN_diag    Output vector of log(mole fraction)
 *                 derivatives of the log Activity Coefficients.
 *                 length = m_kk
 */
void IonsFromNeutralVPSSTP::getdlnActCoeffdlnN_diag(doublereal* dlnActCoeffdlnN_diag) const
{
    s_update_lnActCoeff();
    s_update_dlnActCoeff_dlnN_diag();

    for (size_t k = 0; k < m_kk; k++) {
        dlnActCoeffdlnN_diag[k] = dlnActCoeffdlnN_diag_[k];
    }
}
//====================================================================================================================
void IonsFromNeutralVPSSTP::getdlnActCoeffdlnN(const size_t ld, doublereal* dlnActCoeffdlnN)
{
    s_update_lnActCoeff();
    s_update_dlnActCoeff_dlnN();
    double* data =  & dlnActCoeffdlnN_(0,0);
    for (size_t k = 0; k < m_kk; k++) {
        for (size_t m = 0; m < m_kk; m++) {
            dlnActCoeffdlnN[ld * k + m] = data[m_kk * k + m];
        }
    }
}
//====================================================================================================================
void IonsFromNeutralVPSSTP::setTemperature(const doublereal temp)
{
    double p = pressure();
    IonsFromNeutralVPSSTP::setState_TP(temp, p);
}
//====================================================================================================================
void IonsFromNeutralVPSSTP::setPressure(doublereal p)
{
    double t = temperature();
    IonsFromNeutralVPSSTP::setState_TP(t, p);
}
//====================================================================================================================
// Set the temperature (K) and pressure (Pa)
/*
 * Setting the pressure may involve the solution of a nonlinear equation.
 *
 * @param t    Temperature (K)
 * @param p    Pressure (Pa)
 */
void IonsFromNeutralVPSSTP::setState_TP(doublereal t, doublereal p)
{
    /*
     *  This is a two phase process. First, we calculate the standard states
     *  within the neutral molecule phase.
     */
    neutralMoleculePhase_->setState_TP(t, p);
    VPStandardStateTP::setState_TP(t,p);

    /*
     * Calculate the partial molar volumes, and then the density of the fluid
     */

    //calcDensity();
    double dd = neutralMoleculePhase_->density();
    Phase::setDensity(dd);
}

// Calculate ion mole fractions from neutral molecule
// mole fractions.
/*
 *  @param mf Dump the mole fractions into this vector.
 */
void IonsFromNeutralVPSSTP::calcIonMoleFractions(doublereal* const mf) const
{
    doublereal fmij;
    /*
     * Download the neutral mole fraction vector into the
     * vector, NeutralMolecMoleFractions_[]
     */
    neutralMoleculePhase_->getMoleFractions(DATA_PTR(NeutralMolecMoleFractions_));

    // Zero the mole fractions
    for (size_t k = 0; k < m_kk; k++) {
        mf[k] = 0.0;
    }

    /*
     *  Use the formula matrix to calculate the relative mole numbers.
     */
    for (size_t jNeut = 0; jNeut <  numNeutralMoleculeSpecies_; jNeut++) {
        for (size_t k = 0; k < m_kk; k++) {
            fmij =  fm_neutralMolec_ions_[k + jNeut * m_kk];
            mf[k] += fmij * NeutralMolecMoleFractions_[jNeut];
        }
    }

    /*
     * Normalize the new mole fractions
     */
    doublereal sum = 0.0;
    for (size_t k = 0; k < m_kk; k++) {
        sum += mf[k];
    }
    for (size_t k = 0; k < m_kk; k++) {
        mf[k] /= sum;
    }

}
//====================================================================================================================
// Calculate neutral molecule mole fractions
/*
 *  This routine calculates the neutral molecule mole
 *  fraction given the vector of ion mole fractions,
 *  i.e., the mole fractions from this ThermoPhase.
 *  Note, this routine basically assumes that there
 *  is charge neutrality. If there isn't, then it wouldn't
 *  make much sense.
 *
 *  for the case of  cIonSolnType_SINGLEANION, some slough
 *  in the charge neutrality is allowed. The cation number
 *  is followed, while the difference in charge neutrality
 *  is dumped into the anion mole number to fix the imbalance.
 */
void IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions() const
{
    size_t icat, jNeut;
    doublereal fmij;
    doublereal sum = 0.0;

    //! Zero the vector we are trying to find.
    for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
        NeutralMolecMoleFractions_[k] = 0.0;
    }
#ifdef DEBUG_MODE
    sum = -1.0;
    for (int k = 0; k < m_kk; k++) {
        sum += moleFractions_[k];
    }
    if (fabs(sum) > 1.0E-11)  {
        throw CanteraError("IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions",
                           "molefracts don't sum to one: " + fp2str(sum));
    }
#endif

    // bool fmSimple = true;

    switch (ionSolnType_) {

    case cIonSolnType_PASSTHROUGH:
        for (size_t k = 0; k < m_kk; k++) {
            NeutralMolecMoleFractions_[k] = moleFractions_[k];
        }
        break;

    case cIonSolnType_SINGLEANION:
        for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
            NeutralMolecMoleFractions_[k] = 0.0;
        }

        for (size_t k = 0; k < cationList_.size(); k++) {
            //! Get the id for the next cation
            icat = cationList_[k];
            jNeut = fm_invert_ionForNeutral[icat];
            if (jNeut != npos) {
                fmij =  fm_neutralMolec_ions_[icat + jNeut * m_kk];
                AssertTrace(fmij != 0.0);
                NeutralMolecMoleFractions_[jNeut] += moleFractions_[icat] / fmij;
            }
        }

        for (size_t k = 0; k <  numPassThroughSpecies_; k++) {
            icat = passThroughList_[k];
            jNeut = fm_invert_ionForNeutral[icat];
            fmij = fm_neutralMolec_ions_[ icat + jNeut * m_kk];
            NeutralMolecMoleFractions_[jNeut] += moleFractions_[icat] / fmij;
        }

#ifdef DEBUG_MODE
        for (size_t k = 0; k < m_kk; k++) {
            moleFractionsTmp_[k] = moleFractions_[k];
        }
        for (jNeut = 0; jNeut <  numNeutralMoleculeSpecies_; jNeut++) {
            for (size_t k = 0; k < m_kk; k++) {
                fmij =  fm_neutralMolec_ions_[k + jNeut * m_kk];
                moleFractionsTmp_[k] -= fmij * NeutralMolecMoleFractions_[jNeut];
            }
        }
        for (size_t k = 0; k < m_kk; k++) {
            if (fabs(moleFractionsTmp_[k]) > 1.0E-13) {
                //! Check to see if we have in fact found the inverse.
                if (anionList_[0] != k) {
                    throw CanteraError("", "neutral molecule calc error");
                } else {
                    //! For the single anion case, we will allow some slippage
                    if (fabs(moleFractionsTmp_[k]) > 1.0E-5) {
                        throw CanteraError("", "neutral molecule calc error - anion");
                    }
                }
            }
        }
#endif

        // Normalize the Neutral Molecule mole fractions
        sum = 0.0;
        for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
            sum += NeutralMolecMoleFractions_[k];
        }
        for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
            NeutralMolecMoleFractions_[k] /= sum;
        }

        break;

    case  cIonSolnType_SINGLECATION:

        throw CanteraError("eosType", "Unknown type");

        break;

    case  cIonSolnType_MULTICATIONANION:

        throw CanteraError("eosType", "Unknown type");
        break;

    default:

        throw CanteraError("eosType", "Unknown type");
        break;

    }
}
//====================================================================================================================
// Calculate neutral molecule mole fractions
/*
 *  This routine calculates the neutral molecule mole
 *  fraction given the vector of ion mole fractions,
 *  i.e., the mole fractions from this ThermoPhase.
 *  Note, this routine basically assumes that there
 *  is charge neutrality. If there isn't, then it wouldn't
 *  make much sense.
 *
 *  for the case of  cIonSolnType_SINGLEANION, some slough
 *  in the charge neutrality is allowed. The cation number
 *  is followed, while the difference in charge neutrality
 *  is dumped into the anion mole number to fix the imbalance.
 */
void IonsFromNeutralVPSSTP::getNeutralMoleculeMoleGrads(const doublereal* const dx, doublereal* const dy) const
{
    doublereal fmij;
    vector_fp y;
    y.resize(numNeutralMoleculeSpecies_,0.0);
    doublereal sumy, sumdy;

    //check sum dx = 0

    //! Zero the vector we are trying to find.
    for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
        dy[k] = 0.0;
    }


    // bool fmSimple = true;

    switch (ionSolnType_) {

    case cIonSolnType_PASSTHROUGH:
        for (size_t k = 0; k < m_kk; k++) {
            dy[k] = dx[k];
        }
        break;

    case cIonSolnType_SINGLEANION:
        for (size_t k = 0; k < cationList_.size(); k++) {
            //! Get the id for the next cation
            size_t icat = cationList_[k];
            size_t jNeut = fm_invert_ionForNeutral[icat];
            if (jNeut != npos) {
                fmij =  fm_neutralMolec_ions_[icat + jNeut * m_kk];
                AssertTrace(fmij != 0.0);
                dy[jNeut] += dx[icat] / fmij;
                y[jNeut] += moleFractions_[icat] / fmij;
            }
        }

        for (size_t k = 0; k <  numPassThroughSpecies_; k++) {
            size_t icat = passThroughList_[k];
            size_t jNeut = fm_invert_ionForNeutral[icat];
            fmij = fm_neutralMolec_ions_[ icat + jNeut * m_kk];
            dy[jNeut] += dx[icat] / fmij;
            y[jNeut] += moleFractions_[icat] / fmij;
        }
#ifdef DEBUG_MODE_NOT
        //check dy sum to zero
        for (size_t k = 0; k < m_kk; k++) {
            moleFractionsTmp_[k] = dx[k];
        }
        for (jNeut = 0; jNeut <  numNeutralMoleculeSpecies_; jNeut++) {
            for (size_t k = 0; k < m_kk; k++) {
                fmij =  fm_neutralMolec_ions_[k + jNeut * m_kk];
                moleFractionsTmp_[k] -= fmij * dy[jNeut];
            }
        }
        for (size_t k = 0; k < m_kk; k++) {
            if (fabs(moleFractionsTmp_[k]) > 1.0E-13) {
                //! Check to see if we have in fact found the inverse.
                if (anionList_[0] != k) {
                    throw CanteraError("", "neutral molecule calc error");
                } else {
                    //! For the single anion case, we will allow some slippage
                    if (fabs(moleFractionsTmp_[k]) > 1.0E-5) {
                        throw CanteraError("", "neutral molecule calc error - anion");
                    }
                }
            }
        }
#endif
        // Normalize the Neutral Molecule mole fractions
        sumy = 0.0;
        sumdy = 0.0;
        for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
            sumy += y[k];
            sumdy += dy[k];
        }
        for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
            dy[k] = dy[k]/sumy - y[k]*sumdy/sumy/sumy;
        }

        break;

    case  cIonSolnType_SINGLECATION:

        throw CanteraError("eosType", "Unknown type");

        break;

    case  cIonSolnType_MULTICATIONANION:

        throw CanteraError("eosType", "Unknown type");
        break;

    default:

        throw CanteraError("eosType", "Unknown type");
        break;

    }
}


void IonsFromNeutralVPSSTP::setMassFractions(const doublereal* const y)
{
    GibbsExcessVPSSTP::setMassFractions(y);
    calcNeutralMoleculeMoleFractions();
    neutralMoleculePhase_->setMoleFractions(DATA_PTR(NeutralMolecMoleFractions_));
}

void IonsFromNeutralVPSSTP::setMassFractions_NoNorm(const doublereal* const y)
{
    GibbsExcessVPSSTP::setMassFractions_NoNorm(y);
    calcNeutralMoleculeMoleFractions();
    neutralMoleculePhase_->setMoleFractions(DATA_PTR(NeutralMolecMoleFractions_));
}

void IonsFromNeutralVPSSTP::setMoleFractions(const doublereal* const x)
{
    GibbsExcessVPSSTP::setMoleFractions(x);
    calcNeutralMoleculeMoleFractions();
    neutralMoleculePhase_->setMoleFractions(DATA_PTR(NeutralMolecMoleFractions_));
}

void IonsFromNeutralVPSSTP::setMoleFractions_NoNorm(const doublereal* const x)
{
    GibbsExcessVPSSTP::setMoleFractions_NoNorm(x);
    calcNeutralMoleculeMoleFractions();
    neutralMoleculePhase_->setMoleFractions_NoNorm(DATA_PTR(NeutralMolecMoleFractions_));
}


void IonsFromNeutralVPSSTP::setConcentrations(const doublereal* const c)
{
    GibbsExcessVPSSTP::setConcentrations(c);
    calcNeutralMoleculeMoleFractions();
    neutralMoleculePhase_->setMoleFractions(DATA_PTR(NeutralMolecMoleFractions_));
}

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


doublereal IonsFromNeutralVPSSTP::err(std::string msg) const
{
    throw CanteraError("IonsFromNeutralVPSSTP","Base class method "
                       +msg+" called. Equation of state type: "+int2str(eosType()));
    return 0;
}
/*
 *   Import, construct, and initialize a phase
 *   specification from an XML tree into the current object.
 *
 * This routine is a precursor to constructPhaseXML(XML_Node*)
 * routine, which does most of the work.
 *
 * @param infile 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 IonsFromNeutralVPSSTP::constructPhaseFile(std::string inputFile, std::string id)
{

    if (inputFile.size() == 0) {
        throw CanteraError("MargulesVPSSTP:constructPhaseFile",
                           "input file is null");
    }
    string path = findInputFile(inputFile);
    std::ifstream fin(path.c_str());
    if (!fin) {
        throw CanteraError("MargulesVPSSTP: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("MargulesVPSSTP:constructPhaseFile",
                           "ERROR: Can not find phase named " +
                           id + " in file named " + inputFile);
    }
    fxml_phase->copy(&phaseNode_XML);
    constructPhaseXML(*fxml_phase, id);
    delete fxml;
}

/*
 *   Import, construct, and initialize a HMWSoln phase
 *   specification from an XML tree into the current object.
 *
 *   Most of the work is carried out by the cantera base
 *   routine, importPhase(). That routine imports all of the
 *   species and element data, including the standard states
 *   of the species.
 *
 *   Then, In this routine, we read the information
 *   particular to the specification of the activity
 *   coefficient model for the Pitzer parameterization.
 *
 *   We also read information about the molar volumes of the
 *   standard states if present in the XML file.
 *
 * @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 IonsFromNeutralVPSSTP::constructPhaseXML(XML_Node& phaseNode, std::string id)
{
    string stemp;
    if (id.size() > 0) {
        string idp = phaseNode.id();
        if (idp != id) {
            throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML",
                               "phasenode and Id are incompatible");
        }
    }

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



    /*
     * Make sure that the thermo model is IonsFromNeutralMolecule
     */
    stemp = thermoNode.attrib("model");
    string formString = lowercase(stemp);
    if (formString != "ionsfromneutralmolecule") {
        throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML",
                           "model name isn't IonsFromNeutralMolecule: " + formString);
    }

    /*
     * Find the Neutral Molecule Phase
     */
    if (!thermoNode.hasChild("neutralMoleculePhase")) {
        throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML",
                           "no neutralMoleculePhase XML node");
    }
    XML_Node& neutralMoleculeNode = thermoNode.child("neutralMoleculePhase");

    string nsource = neutralMoleculeNode["datasrc"];
    XML_Node* neut_ptr = get_XML_Node(nsource, 0);
    if (!neut_ptr) {
        throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML",
                           "neut_ptr = 0");
    }

    /*
     *  Create the neutralMolecule ThermoPhase if we haven't already
     */
    if (!neutralMoleculePhase_) {
        neutralMoleculePhase_  = newPhase(*neut_ptr);
    }

    /*
     * Call the Cantera importPhase() function. This will import
     * all of the species into the phase. This will also handle
     * all of the solvent and solute standard states
     */
    bool m_ok = importPhase(phaseNode, this);
    if (!m_ok) {
        throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML",
                           "importPhase failed ");
    }

}

//====================================================================================================================
/*
 * @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
 */
void IonsFromNeutralVPSSTP::initThermo()
{
    initLengths();
    GibbsExcessVPSSTP::initThermo();
}
//====================================================================================================================

//   Initialize lengths of local variables after all species have
//   been identified.
void  IonsFromNeutralVPSSTP::initLengths()
{
    m_kk = nSpecies();
    numNeutralMoleculeSpecies_ =  neutralMoleculePhase_->nSpecies();
    moleFractions_.resize(m_kk);
    fm_neutralMolec_ions_.resize(numNeutralMoleculeSpecies_ * m_kk);
    fm_invert_ionForNeutral.resize(m_kk);
    NeutralMolecMoleFractions_.resize(numNeutralMoleculeSpecies_);
    cationList_.resize(m_kk);
    anionList_.resize(m_kk);
    passThroughList_.resize(m_kk);
    moleFractionsTmp_.resize(m_kk);
    muNeutralMolecule_.resize(numNeutralMoleculeSpecies_);
    lnActCoeff_NeutralMolecule_.resize(numNeutralMoleculeSpecies_);
    dlnActCoeffdT_NeutralMolecule_.resize(numNeutralMoleculeSpecies_);
    dlnActCoeffdlnX_diag_NeutralMolecule_.resize(numNeutralMoleculeSpecies_);
    dlnActCoeffdlnN_diag_NeutralMolecule_.resize(numNeutralMoleculeSpecies_);
    dlnActCoeffdlnN_NeutralMolecule_.resize(numNeutralMoleculeSpecies_, numNeutralMoleculeSpecies_, 0.0);
}
//====================================================================================================================
//!  Return the factor overlap
/*!
 *     @param elnamesVN
 *     @param elemVectorN
 *     @param nElementsN
 *     @param elnamesVI
 *     @param elemVectorI
 *     @param nElementsI
 *
 */
static double factorOverlap(const std::vector<std::string>&  elnamesVN ,
                            const std::vector<double>& elemVectorN,
                            const size_t nElementsN,
                            const std::vector<std::string>&  elnamesVI ,
                            const std::vector<double>& elemVectorI,
                            const size_t nElementsI)
{
    double fMax = 1.0E100;
    for (size_t mi = 0; mi < nElementsI; mi++) {
        if (elnamesVI[mi] != "E") {
            if (elemVectorI[mi] > 1.0E-13) {
                double eiNum = elemVectorI[mi];
                for (size_t mn = 0; mn < nElementsN; mn++) {
                    if (elnamesVI[mi] == elnamesVN[mn]) {
                        if (elemVectorN[mn] <= 1.0E-13) {
                            return 0.0;
                        }
                        fMax = std::min(fMax, elemVectorN[mn]/eiNum);
                    }
                }
            }
        }
    }
    return fMax;
}
//====================================================================================================================
/*
 * initThermoXML()                (virtual from ThermoPhase)
 *   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 IonsFromNeutralVPSSTP::initThermoXML(XML_Node& phaseNode, std::string id)
{
    size_t k;
    /*
     *   variables that need to be populated
     *
     *    cationList_
     *      numCationSpecies_;
     */

    numCationSpecies_ = 0;
    cationList_.clear();
    for (k = 0; k < m_kk; k++) {
        if (charge(k) > 0) {
            cationList_.push_back(k);
            numCationSpecies_++;
        }
    }

    numAnionSpecies_ = 0;
    anionList_.clear();
    for (k = 0; k < m_kk; k++) {
        if (charge(k) < 0) {
            anionList_.push_back(k);
            numAnionSpecies_++;
        }
    }

    numPassThroughSpecies_= 0;
    passThroughList_.clear();
    for (k = 0; k < m_kk; k++) {
        if (charge(k) == 0) {
            passThroughList_.push_back(k);
            numPassThroughSpecies_++;
        }
    }

    PDSS_IonsFromNeutral* speciesSS = 0;
    indexSpecialSpecies_ = npos;
    for (k = 0; k < m_kk; k++) {
        speciesSS = dynamic_cast<PDSS_IonsFromNeutral*>(providePDSS(k));
        if (!speciesSS) {
            throw CanteraError("initThermoXML", "Dynamic cast failed");
        }
        if (speciesSS->specialSpecies_ == 1) {
            indexSpecialSpecies_ = k;
        }
        if (speciesSS->specialSpecies_ == 2) {
            indexSecondSpecialSpecies_ = k;
        }
    }


    size_t nElementsN =  neutralMoleculePhase_->nElements();
    const std::vector<std::string>&  elnamesVN = neutralMoleculePhase_->elementNames();
    std::vector<double> elemVectorN(nElementsN);
    std::vector<double> elemVectorN_orig(nElementsN);

    size_t nElementsI =  nElements();
    const std::vector<std::string>&  elnamesVI = elementNames();
    std::vector<double> elemVectorI(nElementsI);

    vector<doublereal> fm_tmp(m_kk);
    for (size_t k = 0; k <  m_kk; k++) {
        fm_invert_ionForNeutral[k] = npos;
    }
    /*    for (int jNeut = 0; jNeut <  numNeutralMoleculeSpecies_; jNeut++) {
      fm_invert_ionForNeutral[jNeut] = -1;
      }*/
    for (size_t jNeut = 0; jNeut <  numNeutralMoleculeSpecies_; jNeut++) {
        for (size_t m = 0; m < nElementsN; m++) {
            elemVectorN[m] = neutralMoleculePhase_->nAtoms(jNeut, m);
        }
        elemVectorN_orig = elemVectorN;
        fm_tmp.assign(m_kk, 0.0);

        for (size_t m = 0; m < nElementsI; m++) {
            elemVectorI[m] = nAtoms(indexSpecialSpecies_, m);
        }
        double fac = factorOverlap(elnamesVN, elemVectorN, nElementsN,
                                   elnamesVI ,elemVectorI, nElementsI);
        if (fac > 0.0) {
            for (size_t m = 0; m < nElementsN; m++) {
                std::string mName = elnamesVN[m];
                for (size_t mi = 0; mi < nElementsI; mi++) {
                    std::string eName = elnamesVI[mi];
                    if (mName == eName) {
                        elemVectorN[m] -= fac * elemVectorI[mi];
                    }

                }
            }
        }
        fm_neutralMolec_ions_[indexSpecialSpecies_  + jNeut * m_kk ] += fac;


        for (k = 0; k < m_kk; k++) {
            for (size_t m = 0; m < nElementsI; m++) {
                elemVectorI[m] = nAtoms(k, m);
            }
            double fac = factorOverlap(elnamesVN, elemVectorN, nElementsN,
                                       elnamesVI ,elemVectorI, nElementsI);
            if (fac > 0.0) {
                for (size_t m = 0; m < nElementsN; m++) {
                    std::string mName = elnamesVN[m];
                    for (size_t mi = 0; mi < nElementsI; mi++) {
                        std::string eName = elnamesVI[mi];
                        if (mName == eName) {
                            elemVectorN[m] -= fac * elemVectorI[mi];
                        }

                    }
                }
                bool notTaken = true;
                for (size_t iNeut = 0; iNeut < jNeut; iNeut++) {
                    if (fm_invert_ionForNeutral[k] == iNeut) {
                        notTaken = false;
                    }
                }
                if (notTaken) {
                    fm_invert_ionForNeutral[k] = jNeut;
                } else {
                    throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML",
                                       "Simple formula matrix generation failed, one cation is shared between two salts");
                }
            }
            fm_neutralMolec_ions_[k  + jNeut * m_kk] += fac;
        }

        // Ok check the work
        for (size_t m = 0; m < nElementsN; m++) {
            if (fabs(elemVectorN[m]) > 1.0E-13) {
                throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML",
                                   "Simple formula matrix generation failed");
            }
        }


    }
    /*
     * This includes the setStateFromXML calls
     */
    GibbsExcessVPSSTP::initThermoXML(phaseNode, id);

    /*
     * There is one extra step here. We assure ourselves that we
     * have charge conservation.
     */
}
//====================================================================================================================
// Update the activity coefficients
/*
 * This function will be called to update the internally stored
 * natural logarithm of the activity coefficients
 *
 */
void IonsFromNeutralVPSSTP::s_update_lnActCoeff() const
{
    size_t icat, jNeut;
    doublereal fmij;
    /*
     * Get the activity coefficiens of the neutral molecules
     */
    neutralMoleculePhase_->getLnActivityCoefficients(DATA_PTR(lnActCoeff_NeutralMolecule_));

    switch (ionSolnType_) {
    case cIonSolnType_PASSTHROUGH:
        break;
    case cIonSolnType_SINGLEANION:

        // Do the cation list
        for (size_t k = 0; k < cationList_.size(); k++) {
            //! Get the id for the next cation
            icat = cationList_[k];
            jNeut = fm_invert_ionForNeutral[icat];
            fmij =  fm_neutralMolec_ions_[icat + jNeut * m_kk];
            lnActCoeff_Scaled_[icat] = lnActCoeff_NeutralMolecule_[jNeut] / fmij;
        }

        // Do the anion list
        icat = anionList_[0];
        jNeut = fm_invert_ionForNeutral[icat];
        lnActCoeff_Scaled_[icat]= 0.0;

        // Do the list of neutral molecules
        for (size_t k = 0; k <  numPassThroughSpecies_; k++) {
            icat = passThroughList_[k];
            jNeut = fm_invert_ionForNeutral[icat];
            lnActCoeff_Scaled_[icat] = lnActCoeff_NeutralMolecule_[jNeut];
        }
        break;

    case cIonSolnType_SINGLECATION:
        throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff", "Unimplemented type");
        break;
    case cIonSolnType_MULTICATIONANION:
        throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff", "Unimplemented type");
        break;
    default:
        throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff", "Unimplemented type");
        break;
    }

}
//====================================================================================================================
// Get the change in activity coefficients w.r.t. change in state (temp, mole fraction, etc.) along
// a line in parameter space or along a line in physical space
/*
 *
 * @param dTds           Input of temperature change along the path
 * @param dXds           Input vector of changes in mole fraction along the path. length = m_kk
 *                       Along the path length it must be the case that the mole fractions sum to one.
 * @param dlnActCoeffds  Output vector of the directional derivatives of the
 *                       log Activity Coefficients along the path. length = m_kk
 */
void IonsFromNeutralVPSSTP::getdlnActCoeffds(const doublereal dTds, const doublereal* const dXds,
        doublereal* dlnActCoeffds) const
{
    size_t icat, jNeut;
    doublereal fmij;
    /*
     * Get the activity coefficients of the neutral molecules
     */
    GibbsExcessVPSSTP* geThermo = dynamic_cast<GibbsExcessVPSSTP*>(neutralMoleculePhase_);
    if (!geThermo) {
        for (size_t k = 0; k < m_kk; k++) {
            dlnActCoeffds[k] = dXds[k] / moleFractions_[k];
        }
        return;
    }

    size_t numNeutMolSpec = geThermo->nSpecies();
    vector_fp dlnActCoeff_NeutralMolecule(numNeutMolSpec);
    vector_fp dX_NeutralMolecule(numNeutMolSpec);


    getNeutralMoleculeMoleGrads(DATA_PTR(dXds),DATA_PTR(dX_NeutralMolecule));

    // All mole fractions returned to normal

    geThermo->getdlnActCoeffds(dTds, DATA_PTR(dX_NeutralMolecule), DATA_PTR(dlnActCoeff_NeutralMolecule));

    switch (ionSolnType_) {
    case cIonSolnType_PASSTHROUGH:
        break;
    case cIonSolnType_SINGLEANION:

        // Do the cation list
        for (size_t k = 0; k < cationList_.size(); k++) {
            //! Get the id for the next cation
            icat = cationList_[k];
            jNeut = fm_invert_ionForNeutral[icat];
            fmij =  fm_neutralMolec_ions_[icat + jNeut * m_kk];
            dlnActCoeffds[icat] = dlnActCoeff_NeutralMolecule[jNeut]/fmij;
        }

        // Do the anion list
        icat = anionList_[0];
        jNeut = fm_invert_ionForNeutral[icat];
        dlnActCoeffds[icat]= 0.0;

        // Do the list of neutral molecules
        for (size_t k = 0; k < numPassThroughSpecies_; k++) {
            icat = passThroughList_[k];
            jNeut = fm_invert_ionForNeutral[icat];
            dlnActCoeffds[icat] = dlnActCoeff_NeutralMolecule[jNeut];
        }
        break;

    case cIonSolnType_SINGLECATION:
        throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeffds", "Unimplemented type");
        break;
    case cIonSolnType_MULTICATIONANION:
        throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeffds", "Unimplemented type");
        break;
    default:
        throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeffds", "Unimplemented type");
        break;
    }

}
//====================================================================================================================
// Update the temperature derivative of the ln activity coefficients
/*
 * This function will be called to update the internally stored
 * temperature derivative of the natural logarithm of the activity coefficients
 */
void IonsFromNeutralVPSSTP::s_update_dlnActCoeffdT() const
{
    size_t icat, jNeut;
    doublereal fmij;
    /*
     * Get the activity coefficients of the neutral molecules
     */
    GibbsExcessVPSSTP* geThermo = dynamic_cast<GibbsExcessVPSSTP*>(neutralMoleculePhase_);
    if (!geThermo) {
        dlnActCoeffdT_Scaled_.assign(m_kk, 0.0);
        return;
    }

    geThermo->getdlnActCoeffdT(DATA_PTR(dlnActCoeffdT_NeutralMolecule_));

    switch (ionSolnType_) {
    case cIonSolnType_PASSTHROUGH:
        break;
    case cIonSolnType_SINGLEANION:

        // Do the cation list
        for (size_t k = 0; k < cationList_.size(); k++) {
            //! Get the id for the next cation
            icat = cationList_[k];
            jNeut = fm_invert_ionForNeutral[icat];
            fmij =  fm_neutralMolec_ions_[icat + jNeut * m_kk];
            dlnActCoeffdT_Scaled_[icat] = dlnActCoeffdT_NeutralMolecule_[jNeut]/fmij;
        }

        // Do the anion list
        icat = anionList_[0];
        jNeut = fm_invert_ionForNeutral[icat];
        dlnActCoeffdT_Scaled_[icat]= 0.0;

        // Do the list of neutral molecules
        for (size_t k = 0; k <  numPassThroughSpecies_; k++) {
            icat = passThroughList_[k];
            jNeut = fm_invert_ionForNeutral[icat];
            dlnActCoeffdT_Scaled_[icat] = dlnActCoeffdT_NeutralMolecule_[jNeut];
        }
        break;

    case cIonSolnType_SINGLECATION:
        throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeffdT", "Unimplemented type");
        break;
    case cIonSolnType_MULTICATIONANION:
        throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeffdT", "Unimplemented type");
        break;
    default:
        throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeffdT", "Unimplemented type");
        break;
    }

}
//====================================================================================================================
/*
 * This function will be called to update the internally stored
 * temperature derivative of the natural logarithm of the activity coefficients
 */
void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnX_diag() const
{
    size_t icat, jNeut;
    doublereal fmij;
    /*
     * Get the activity coefficients of the neutral molecules
     */
    GibbsExcessVPSSTP* geThermo = dynamic_cast<GibbsExcessVPSSTP*>(neutralMoleculePhase_);
    if (!geThermo) {
        dlnActCoeffdlnX_diag_.assign(m_kk, 0.0);
        return;
    }

    geThermo->getdlnActCoeffdlnX_diag(DATA_PTR(dlnActCoeffdlnX_diag_NeutralMolecule_));

    switch (ionSolnType_) {
    case cIonSolnType_PASSTHROUGH:
        break;
    case cIonSolnType_SINGLEANION:

        // Do the cation list
        for (size_t k = 0; k < cationList_.size(); k++) {
            //! Get the id for the next cation
            icat = cationList_[k];
            jNeut = fm_invert_ionForNeutral[icat];
            fmij =  fm_neutralMolec_ions_[icat + jNeut * m_kk];
            dlnActCoeffdlnX_diag_[icat] = dlnActCoeffdlnX_diag_NeutralMolecule_[jNeut]/fmij;
        }

        // Do the anion list
        icat = anionList_[0];
        jNeut = fm_invert_ionForNeutral[icat];
        dlnActCoeffdlnX_diag_[icat]= 0.0;

        // Do the list of neutral molecules
        for (size_t k = 0; k <  numPassThroughSpecies_; k++) {
            icat = passThroughList_[k];
            jNeut = fm_invert_ionForNeutral[icat];
            dlnActCoeffdlnX_diag_[icat] = dlnActCoeffdlnX_diag_NeutralMolecule_[jNeut];
        }
        break;

    case cIonSolnType_SINGLECATION:
        throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnX_diag()", "Unimplemented type");
        break;
    case cIonSolnType_MULTICATIONANION:
        throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnX_diag()", "Unimplemented type");
        break;
    default:
        throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnX_diag()", "Unimplemented type");
        break;
    }

}
//====================================================================================================================
/*
 * This function will be called to update the internally stored
 * temperature derivative of the natural logarithm of the activity coefficients
 */
void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnN_diag() const
{
    size_t icat, jNeut;
    doublereal fmij;
    /*
     * Get the activity coefficients of the neutral molecules
     */
    GibbsExcessVPSSTP* geThermo = dynamic_cast<GibbsExcessVPSSTP*>(neutralMoleculePhase_);
    if (!geThermo) {
        dlnActCoeffdlnN_diag_.assign(m_kk, 0.0);
        return;
    }

    geThermo->getdlnActCoeffdlnN_diag(DATA_PTR(dlnActCoeffdlnN_diag_NeutralMolecule_));

    switch (ionSolnType_) {
    case cIonSolnType_PASSTHROUGH:
        break;
    case cIonSolnType_SINGLEANION:

        // Do the cation list
        for (size_t k = 0; k < cationList_.size(); k++) {
            //! Get the id for the next cation
            icat = cationList_[k];
            jNeut = fm_invert_ionForNeutral[icat];
            fmij =  fm_neutralMolec_ions_[icat + jNeut * m_kk];
            dlnActCoeffdlnN_diag_[icat] = dlnActCoeffdlnN_diag_NeutralMolecule_[jNeut]/fmij;
        }

        // Do the anion list
        icat = anionList_[0];
        jNeut = fm_invert_ionForNeutral[icat];
        dlnActCoeffdlnN_diag_[icat]= 0.0;

        // Do the list of neutral molecules
        for (size_t k = 0; k < numPassThroughSpecies_; k++) {
            icat = passThroughList_[k];
            jNeut = fm_invert_ionForNeutral[icat];
            dlnActCoeffdlnN_diag_[icat] = dlnActCoeffdlnN_diag_NeutralMolecule_[jNeut];
        }
        break;

    case cIonSolnType_SINGLECATION:
        throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnN_diag()", "Unimplemented type");
        break;
    case cIonSolnType_MULTICATIONANION:
        throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnN_diag()", "Unimplemented type");
        break;
    default:
        throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnN_diag()", "Unimplemented type");
        break;
    }

}
//====================================================================================================================
// Update the derivative of the log of the activity coefficients
//  wrt log(number of moles) - diagonal components
/*
 * This function will be called to update the internally stored
 * derivative of the natural logarithm of the activity coefficients
 * wrt logarithm of the number of moles of given species.
 */
void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnN() const
{
    size_t kcat, kNeut, mcat, mNeut;
    doublereal fmij, mfmij;
    dlnActCoeffdlnN_.zero();
    /*
     * Get the activity coefficients of the neutral molecules
     */
    GibbsExcessVPSSTP* geThermo = dynamic_cast<GibbsExcessVPSSTP*>(neutralMoleculePhase_);
    if (!geThermo) {
        throw CanteraError("IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnN()", "dynamic cast failed");
    }
    size_t nsp_ge = geThermo->nSpecies();
    geThermo->getdlnActCoeffdlnN(nsp_ge, &(dlnActCoeffdlnN_NeutralMolecule_(0,0)));

    switch (ionSolnType_) {
    case cIonSolnType_PASSTHROUGH:
        break;
    case cIonSolnType_SINGLEANION:

        // Do the cation list
        for (size_t k = 0; k < cationList_.size(); k++) {
            for (size_t m = 0; m < cationList_.size(); m++) {
                kcat = cationList_[k];

                kNeut = fm_invert_ionForNeutral[kcat];
                fmij =  fm_neutralMolec_ions_[kcat + kNeut * m_kk];
                dlnActCoeffdlnN_diag_[kcat] = dlnActCoeffdlnN_diag_NeutralMolecule_[kNeut]/fmij;

                mcat = cationList_[m];
                mNeut = fm_invert_ionForNeutral[mcat];
                mfmij =  fm_neutralMolec_ions_[mcat + mNeut * m_kk];

                dlnActCoeffdlnN_(kcat,mcat) = dlnActCoeffdlnN_NeutralMolecule_(kNeut,mNeut) * mfmij / fmij;

            }
            for (size_t m = 0; m < numPassThroughSpecies_; m++) {
                mcat = passThroughList_[m];
                mNeut = fm_invert_ionForNeutral[mcat];
                dlnActCoeffdlnN_(kcat, mcat) =  dlnActCoeffdlnN_NeutralMolecule_(kNeut, mNeut) / fmij;
            }
        }

        // Do the anion list -> anion activity coefficient is one
        kcat = anionList_[0];
        kNeut = fm_invert_ionForNeutral[kcat];
        for (size_t k = 0; k < m_kk; k++) {
            dlnActCoeffdlnN_(kcat, k) = 0.0;
            dlnActCoeffdlnN_(k, kcat) = 0.0;
        }

        // Do the list of neutral molecules
        for (size_t k = 0; k <  numPassThroughSpecies_; k++) {
            kcat = passThroughList_[k];
            kNeut = fm_invert_ionForNeutral[kcat];
            dlnActCoeffdlnN_diag_[kcat] = dlnActCoeffdlnN_diag_NeutralMolecule_[kNeut];

            for (size_t m = 0; m < m_kk; m++) {
                mcat = passThroughList_[m];
                mNeut = fm_invert_ionForNeutral[mcat];
                dlnActCoeffdlnN_(kcat, mcat) =  dlnActCoeffdlnN_NeutralMolecule_(kNeut, mNeut);
            }


            for (size_t m = 0; m < cationList_.size(); m++) {
                mcat = cationList_[m];
                mNeut = fm_invert_ionForNeutral[mcat];
                mfmij =  fm_neutralMolec_ions_[mcat + mNeut * m_kk];
                dlnActCoeffdlnN_(kcat, mcat) =  dlnActCoeffdlnN_NeutralMolecule_(kNeut,mNeut);
            }

        }
        break;

    case cIonSolnType_SINGLECATION:
        throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnN", "Unimplemented type");
        break;
    case cIonSolnType_MULTICATIONANION:
        throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnN", "Unimplemented type");
        break;
    default:
        throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnN", "Unimplemented type");
        break;
    }
}
//====================================================================================================================
}
//======================================================================================================================