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/base/stringUtils.h"

#include <fstream>

using namespace std;

namespace Cantera
{

IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP() :
    GibbsExcessVPSSTP(),
    ionSolnType_(cIonSolnType_SINGLEANION),
    numNeutralMoleculeSpecies_(0),
    indexSpecialSpecies_(npos),
    indexSecondSpecialSpecies_(npos),
    numCationSpecies_(0),
    numAnionSpecies_(0),
    numPassThroughSpecies_(0),
    neutralMoleculePhase_(0),
    geThermo(0),
    IOwnNThermoPhase_(true),
    moleFractionsTmp_(0),
    muNeutralMolecule_(0),
    lnActCoeff_NeutralMolecule_(0)
{
}

IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP(const std::string& inputFile,
        const 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_);
    geThermo = dynamic_cast<GibbsExcessVPSSTP*>(neutralMoleculePhase_);
    //y.resize(numNeutralMoleculeSpecies_,0.0);
    //size_t numNeutMolSpec = geThermo->nSpecies();
    //dlnActCoeff_NeutralMolecule.resize(numNeutMolSpec);
    //dX_NeutralMolecule.resize(numNeutMolSpec);
}

IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP(XML_Node& phaseRoot,
        const 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_);
    geThermo = dynamic_cast<GibbsExcessVPSSTP*>(neutralMoleculePhase_);
    y_.resize(numNeutralMoleculeSpecies_,0.0);
    size_t numNeutMolSpec = geThermo->nSpecies();
    dlnActCoeff_NeutralMolecule_.resize(numNeutMolSpec);
    dX_NeutralMolecule_.resize(numNeutMolSpec);
}

IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP(const IonsFromNeutralVPSSTP& b) :
    GibbsExcessVPSSTP(),
    ionSolnType_(cIonSolnType_SINGLEANION),
    numNeutralMoleculeSpecies_(0),
    indexSpecialSpecies_(npos),
    indexSecondSpecialSpecies_(npos),
    numCationSpecies_(0),
    numAnionSpecies_(0),
    numPassThroughSpecies_(0),
    neutralMoleculePhase_(0),
    geThermo(0),
    IOwnNThermoPhase_(true),
    moleFractionsTmp_(0),
    muNeutralMolecule_(0),

    lnActCoeff_NeutralMolecule_(0)
{
    IonsFromNeutralVPSSTP::operator=(b);
}

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_) {
            delete neutralMoleculePhase_;
            neutralMoleculePhase_   = (b.neutralMoleculePhase_)->duplMyselfAsThermoPhase();
        } else {
            neutralMoleculePhase_   = 0;
        }
    } else {
        neutralMoleculePhase_     = b.neutralMoleculePhase_;
    }
    geThermo = dynamic_cast<GibbsExcessVPSSTP*>(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_;

    y_                          = b.y_;
    dlnActCoeff_NeutralMolecule_ = b.dlnActCoeff_NeutralMolecule_;
    dX_NeutralMolecule_         = b.dX_NeutralMolecule_;

    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::~IonsFromNeutralVPSSTP()
{
    if (IOwnNThermoPhase_) {
        delete neutralMoleculePhase_;
        neutralMoleculePhase_ = 0;
    }
}

ThermoPhase*
IonsFromNeutralVPSSTP::duplMyselfAsThermoPhase() const
{
    return new IonsFromNeutralVPSSTP(*this);
}

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;
}

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 ");
    }

}

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

int IonsFromNeutralVPSSTP::eosType() const
{
    return cIonsFromNeutral;
}

/*
 * ------------ Molar Thermodynamic Properties ----------------------
 */

doublereal IonsFromNeutralVPSSTP::enthalpy_mole() const
{
    getPartialMolarEnthalpies(DATA_PTR(m_pp));
    return mean_X(DATA_PTR(m_pp));
}

doublereal IonsFromNeutralVPSSTP::intEnergy_mole() const
{
    double hh = enthalpy_mole();
    double pres = pressure();
    double molarV = 1.0/molarDensity();
    return hh - pres * molarV;
}

doublereal IonsFromNeutralVPSSTP::entropy_mole() const
{
    getPartialMolarEntropies(DATA_PTR(m_pp));
    return mean_X(DATA_PTR(m_pp));
}

doublereal IonsFromNeutralVPSSTP::gibbs_mole() const
{
    getChemPotentials(DATA_PTR(m_pp));
    return mean_X(DATA_PTR(m_pp));
}

doublereal IonsFromNeutralVPSSTP::cp_mole() const
{
    getPartialMolarCp(DATA_PTR(m_pp));
    return mean_X(DATA_PTR(m_pp));
}

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;
}

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 -------------
 */

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;
    }
}

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];
    }
}

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], SmallNumber);
        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;
    }
}

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];
    }
}

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);
}

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);
}

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;
    }

}

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 (size_t 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;

    }
}

void IonsFromNeutralVPSSTP::getNeutralMoleculeMoleGrads(const doublereal* const dx, doublereal* const dy) const
{
    doublereal fmij;
    doublereal sumy, sumdy;

    //check sum dx = 0

    //! Zero the vector we are trying to find.
    for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
        y_[k] = 0.0;
        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);
                const doublereal temp = 1.0/fmij;
                dy[jNeut] += dx[icat] * temp;
                y_[jNeut] += moleFractions_[icat] * temp;
            }
        }

        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];
            const doublereal temp = 1.0/fmij;
            dy[jNeut] += dx[icat] * temp;
            y_[jNeut] += moleFractions_[icat] * temp;
        }
#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];
        }
        sumy = 1.0 / sumy;
        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(const std::string& msg) const
{
    throw CanteraError("IonsFromNeutralVPSSTP","Base class method "
                       +msg+" called. Equation of state type: "+int2str(eosType()));
    return 0;
}

void IonsFromNeutralVPSSTP::initThermo()
{
    initLengths();
    GibbsExcessVPSSTP::initThermo();
}

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);

    y_.resize(numNeutralMoleculeSpecies_, 0.0);
    dlnActCoeff_NeutralMolecule_.resize(numNeutralMoleculeSpecies_, 0.0);
    dX_NeutralMolecule_.resize(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;
}
void IonsFromNeutralVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id_)
{
    string stemp;
    if (id_.size() > 0) {
        string idp = phaseNode.id();
        if (idp != id_) {
            throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML",
                               "phasenode and Id are incompatible");
        }
    }

    /*
     * Find the Thermo XML node
     */
    if (!phaseNode.hasChild("thermo")) {
        throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML",
                           "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::initThermoXML",
                           "model name isn't IonsFromNeutralMolecule: " + formString);
    }

    /*
     * Find the Neutral Molecule Phase
     */
    if (!thermoNode.hasChild("neutralMoleculePhase")) {
        throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML",
                           "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::initThermoXML",
                           "neut_ptr = 0");
    }

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

    /*
     *   variables that need to be populated
     *
     *    cationList_
     *      numCationSpecies_;
     */

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

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

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

    PDSS_IonsFromNeutral* speciesSS = 0;
    indexSpecialSpecies_ = npos;
    for (size_t 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 (size_t k = 0; k < m_kk; k++) {
            for (size_t m = 0; m < nElementsI; m++) {
                elemVectorI[m] = nAtoms(k, m);
            }
            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.
     */
}

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;
    }

}

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
     */
    if (!geThermo) {
        for (size_t k = 0; k < m_kk; k++) {
            dlnActCoeffds[k] = dXds[k] / moleFractions_[k];
        }
        return;
    }

    //    static vector_fp dlnActCoeff_NeutralMolecule(numNeutMolSpec);
    //    static 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;
    }

}

void IonsFromNeutralVPSSTP::s_update_dlnActCoeffdT() const
{
    size_t icat, jNeut;
    doublereal fmij;
    /*
     * Get the activity coefficients of the neutral molecules
     */
    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;
    }

}

void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnX_diag() const
{
    size_t icat, jNeut;
    doublereal fmij;
    /*
     * Get the activity coefficients of the neutral molecules
     */
    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;
    }

}

void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnN_diag() const
{
    size_t icat, jNeut;
    doublereal fmij;
    /*
     * Get the activity coefficients of the neutral molecules
     */
    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;
    }

}

void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnN() const
{
    size_t kcat = 0, kNeut = 0, mcat = 0, mNeut = 0;
    doublereal fmij = 0.0, mfmij;
    dlnActCoeffdlnN_.zero();
    /*
     * Get the activity coefficients of the neutral molecules
     */
    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;
    }
}

}