/**
 * @file vcs_solve_TP.cpp Implementation file that contains the
 *     main algorithm for finding an equilibrium
 */
/*
 * Copyright (2005) Sandia Corporation. Under the terms of
 * Contract DE-AC04-94AL85000 with Sandia Corporation, the
 * U.S. Government retains certain rights in this software.
 */

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cassert>

#include "cantera/equil/vcs_solve.h"
#include "cantera/equil/vcs_internal.h"
#include "cantera/equil/vcs_VolPhase.h"
#include "cantera/equil/vcs_species_thermo.h"
#include "cantera/base/ctexceptions.h"

#include "cantera/base/clockWC.h"
#include "cantera/base/stringUtils.h"

using namespace std;
using namespace Cantera;

#ifndef MAX
#define MAX(x,y) (( (x) > (y) ) ? (x) : (y))
#endif

namespace VCSnonideal
{


/***************************************************************************/
/************ Prototypes for static functions ******************************/


static void print_space(size_t num);



#ifdef DEBUG_MODE
#  ifdef DEBUG_NOT
static void prneav(void);
static int  prnfm(void);
#  endif
#endif
/*****************************************************************************/


#ifdef DEBUG_MODE
void VCS_SOLVE::checkDelta1(double* const dsLocal,
                            double* const delTPhMoles, int kspec)
{
    std::vector<double> dchange(m_numPhases, 0.0);
    for (int k = 0; k < kspec; k++) {
        if (m_speciesUnknownType[k] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
            int iph = m_phaseID[k];
            dchange[iph] += dsLocal[k];
        }
    }
    for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
        double denom = MAX(m_totalMolNum, 1.0E-4);
        if (!vcs_doubleEqual(dchange[iphase]/denom, delTPhMoles[iphase]/denom)) {
            plogf("checkDelta1: we have found a problem\n");
            exit(EXIT_FAILURE);
        }
    }
}
#endif
/*****************************************************************************/

// Main routine that solves for equilibrium at constant T and P
// using a variant of the VCS method
/*
 * This is the main routine that solves for equilibrium at constant T and P
 * using a variant of the VCS method. Nonideal phases can be accommodated
 * as well.
 *
 * Any number of single-species phases and  multi-species phases
 * can be handled by the present version.
 *
 *     Input
 * ------------
 *   @param print_lvl     1 -> Print results to standard output
 *                        0 -> don't report on anything
 *
 *   @param printDetails  1 -> Print intermediate results.
 *
 *   @param maxit         Maximum number of iterations for the algorithm
 *
 *   @return     0 = Equilibrium Achieved
 *               1 = Range space error encountered. The element abundance criteria are
 *                   only partially satisfied. Specifically, the first NC= (number of
 *                   components) conditions are satisfied. However, the full NE
 *                   (number of elements) conditions are not satisfied. The equilibrium
 *                   condition is returned.
 *              -1 = Maximum number of iterations is exceeded. Convergence was not
 *                   found.
 */
int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit)
{
    int retn = VCS_SUCCESS, soldel, solveFail;
    double test, RT;
    size_t j, k, l, l1, kspec, irxn, i;
    bool conv = false, allMinorZeroedSpecies = false, forced, lec;
    size_t iph;
    double dx, xx, par;
    size_t it1 = 0;
    size_t npb, iti, lnospec;
    bool dofast;
    int rangeErrorFound = 0;
    bool giveUpOnElemAbund = false;
    int finalElemAbundAttempts = 0;
    bool uptodate_minors = true;
    bool justDeletedMultiPhase = false;
    bool usedZeroedSpecies; /* return flag from basopt indicating that
                  one of the components had a zero concentration */
    vcs_VolPhase* Vphase;
    double*     sc_irxn = NULL;  /* Stoichiometric coefficients for cur rxn  */
    double* dnPhase_irxn;
    double atomComp;
    size_t iphasePop;
    int forceComponentCalc = 1;
    size_t iphaseDelete;  /* integer that determines which phase is being deleted */
    std::vector<size_t> phasePopPhaseIDs(0);
    size_t doPhaseDeleteIph = npos;
    size_t doPhaseDeleteKspec = npos;

#ifdef DEBUG_MODE
    char ANOTE[128];
    /*
     * Set the debug print lvl to the same as the print lvl.
     */
    m_debug_print_lvl = printDetails;
#endif
    if (printDetails > 0 && print_lvl == 0) {
        print_lvl = 1;
    }
    /*
     *    Initialize and set up all counters
     */
    vcs_counters_init(0);
    Cantera::clockWC ticktock;

    /*
     *  Malloc temporary space for usage in this routine and in
     *  subroutines
     *        sm[ne*ne]
     *        ss[ne]
     *        sa[ne]
     *        aw[m]
     *        wx[ne]
     *        xy[m]
     */


    std::vector<double> sm(m_numElemConstraints*m_numElemConstraints, 0.0);
    std::vector<double> ss(m_numElemConstraints, 0.0);
    std::vector<double> sa(m_numElemConstraints, 0.0);

    std::vector<double> aw(m_numSpeciesTot, 0.0);
    std::vector<double> wx(m_numElemConstraints, 0.0);

    solveFail = false;

#if DEBUG_MODE
    int ll;
#endif
    /* ****************************************************** */
    /* **** Evaluate the elemental composition         ****** */
    /* ****************************************************** */
    vcs_elab();

    /* ******************************************************* */
    /* **** Printout the initial conditions for problem ****** */
    /* ******************************************************* */

    if (print_lvl != 0) {
        plogf("VCS CALCULATION METHOD\n\n ");
        plogf("%s\n", m_title.c_str());
        plogf("\n\n%5d SPECIES\n%5d ELEMENTS\n", m_numSpeciesTot, m_numElemConstraints);
        plogf("%5d COMPONENTS\n", m_numComponents);
        plogf("%5d PHASES\n", m_numPhases);

        plogf(" PRESSURE%22.8g %3s\n", m_pressurePA, "Pa ");
        plogf(" TEMPERATURE%19.3f K\n", m_temperature);
        Vphase = m_VolPhaseList[0];
        if (Vphase->nSpecies() > 0) {
            plogf(" PHASE1 INERTS%17.3f\n", TPhInertMoles[0]);
        }
        if (m_numPhases > 1) {
            plogf(" PHASE2 INERTS%17.3f\n", TPhInertMoles[1]);
        }
        plogf("\n ELEMENTAL ABUNDANCES             CORRECT");
        plogf("          FROM ESTIMATE           Type\n\n");
        for (i = 0; i < m_numElemConstraints; ++i) {
            print_space(26);
            plogf("%-2.2s", (m_elementName[i]).c_str());
            plogf("%20.12E%20.12E     %3d\n", m_elemAbundancesGoal[i], m_elemAbundances[i],
                  m_elType[i]);
        }
        if (m_doEstimateEquil < 0) {
            plogf("\n MODIFIED LINEAR PROGRAMMING ESTIMATE OF EQUILIBRIUM - forced\n");
        } else if (m_doEstimateEquil > 0) {
            plogf("\n MODIFIED LINEAR PROGRAMMING ESTIMATE OF EQUILIBRIUM - where necessary\n");
        }
        if (m_doEstimateEquil == 0) {
            plogf("\n USER ESTIMATE OF EQUILIBRIUM\n");
        }
        if (m_VCS_UnitsFormat == VCS_UNITS_KCALMOL) {
            plogf(" Stan. Chem. Pot. in kcal/mole\n");
        }
        if (m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) {
            plogf(" Stan. Chem. Pot. is MU/RT\n");
        }
        if (m_VCS_UnitsFormat == VCS_UNITS_KJMOL) {
            plogf(" Stan. Chem. Pot. in KJ/mole\n");
        }
        if (m_VCS_UnitsFormat == VCS_UNITS_KELVIN) {
            plogf(" Stan. Chem. Pot. in Kelvin\n");
        }
        if (m_VCS_UnitsFormat == VCS_UNITS_MKS) {
            plogf(" Stan. Chem. Pot. in J/kmol\n");
        }
        plogf("\n SPECIES            FORMULA VECTOR   ");
        print_space(41);
        plogf("   STAN_CHEM_POT  EQUILIBRIUM_EST.  Species_Type\n\n");
        print_space(20);
        for (i = 0; i < m_numElemConstraints; ++i) {
            plogf("%-4.4s    ", m_elementName[i].c_str());
        }
        plogf("   PhaseID\n");
        RT = vcs_nondimMult_TP(m_VCS_UnitsFormat, m_temperature);
        for (i = 0; i < m_numSpeciesTot; ++i) {
            plogf(" %-18.18s", m_speciesName[i].c_str());
            for (j = 0; j < m_numElemConstraints; ++j) {
                plogf("% -7.3g ", m_formulaMatrix[j][i]);
            }
            plogf("  %3d  ", m_phaseID[i]);
            print_space(std::max(55-int(m_numElemConstraints)*8, 0));
            plogf("%12.5E  %12.5E", RT * m_SSfeSpecies[i], m_molNumSpecies_old[i]);
            if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_MOLNUM) {
                plogf("       Mol_Num");
            } else if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
                plogf("       Voltage");
            } else {
                plogf("       Unknown");
            }
            plogendl();
        }
    }

    for (i = 0; i < m_numSpeciesTot; ++i) {
        if (m_molNumSpecies_old[i] < 0.0) {
            plogf("On Input species %-12s has a "
                  "negative MF, setting it small",
                  m_speciesName[i].c_str());
            plogendl();
            iph = m_phaseID[i];
            double tmp = m_tPhaseMoles_old[iph] * VCS_RELDELETE_SPECIES_CUTOFF * 10;
            if (VCS_DELETE_MINORSPECIES_CUTOFF*10. > tmp) {
                tmp = VCS_DELETE_MINORSPECIES_CUTOFF*10.;
            }
            m_molNumSpecies_old[i] = tmp;
        }
    }

    /*
     *  Evaluate the total moles of species in the problem
     */
    vcs_tmoles();

    /*
     * Evaluate all chemical potentials at the old mole numbers at the
     * outset of the calculation.
     */
    vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
    vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc);

    /* *********************************************************** */
    /* **** DETERMINE BASIS SPECIES, EVALUATE STOICHIOMETRY ****** */
    /* *********************************************************** */
    /*
     *   This is an entry point for later in the calculation
     */
L_COMPONENT_CALC:
    ;
    test = -1.0e-10;
    retn = vcs_basopt(false, VCS_DATA_PTR(aw), VCS_DATA_PTR(sa),
                      VCS_DATA_PTR(sm), VCS_DATA_PTR(ss),
                      test, &usedZeroedSpecies);
    if (retn != VCS_SUCCESS) {
        return retn;
    }

    // Update the phase objects with the contents of the soln vector
    vcs_updateVP(VCS_STATECALC_OLD);
    vcs_deltag(0, false, VCS_STATECALC_OLD);
    // Turn off the force componentCalc flag
    forceComponentCalc = 0;

    if (conv) {
        goto L_RETURN_BLOCK;
    }
    it1 = 1;

    /*************************************************************************/
    /************** EVALUATE INITIAL SPECIES STATUS VECTOR *******************/
    /*************************************************************************/
    allMinorZeroedSpecies = vcs_evaluate_speciesType();
    lec = false;

    /*************************************************************************/
    /************** EVALUATE THE ELELEMT ABUNDANCE CHECK    ******************/
    /*************************************************************************/
    if (! vcs_elabcheck(0)) {
#ifdef DEBUG_MODE
        if (m_debug_print_lvl >= 2) {
            plogf("   --- Element Abundance check failed");
            plogendl();
        }
#endif
        vcs_elcorr(VCS_DATA_PTR(sm), VCS_DATA_PTR(wx));
        vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
        vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc);
        // Update the phase objects with the contents of the soln vector
        vcs_updateVP(VCS_STATECALC_OLD);
        vcs_deltag(0, false, VCS_STATECALC_OLD);
    }
#ifdef DEBUG_MODE
    else {
        if (m_debug_print_lvl >= 2) {
            plogf("   --- Element Abundance check passed");
            plogendl();
        }
    }
#endif

    iti = 0;
    goto L_MAINLOOP_ALL_SPECIES;

    /* ********************************************************* */
    /* **** SET INITIAL VALUES FOR ITERATION ******************* */
    /* **** EVALUATE REACTION ADJUSTMENTS    ******************* */
    /* ********************************************************* */
    /*
     *  This is the top of the loop ----------------------------------------
     *  Every 4th iteration ITI = 0. Else, It's equal to a negative number
     */
L_MAINLOOP_MM4_SPECIES:
    ;
    iti = ((it1/4) *4) - it1;
    /*
     *      Entry point when the code wants to force an ITI=0 calculation
     */
L_MAINLOOP_ALL_SPECIES:
    ;
    if (iti == 0) {
        /*
         *          Evaluate the minor non-component species chemical
         *          potentials and delta G for their formation reactions
         *          We have already evaluated the major non-components
         */
        if (!uptodate_minors) {
            vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
            vcs_dfe(VCS_STATECALC_OLD, 1, 0, m_numSpeciesRdc);
            vcs_deltag(1, false, VCS_STATECALC_OLD);
        }
        uptodate_minors = true;
    } else {
        uptodate_minors = false;
    }

    if (printDetails) {
        plogf("\n");
        vcs_print_line("=", 110);
        plogf(" Iteration = %3d, Iterations since last evaluation of "
              "optimal basis = %3d",
              m_VCount->Its, it1 - 1);
        if (iti == 0) {
            plogf(" (all species)\n");
        } else {
            plogf(" (only major species)\n");
        }
    }
    /*
     *  Calculate the total moles in each phase -> old solution
     *   -> Needed for numerical stability when phases disappear.
     *   -> the phase moles tend to drift off without this step.
     */
#ifdef DEBUG_MODE
    check_tmoles();
#endif
    vcs_tmoles();
    /*************************************************************************/
    /************** COPY OLD into NEW and ZERO VECTORS ***********************/
    /*************************************************************************/
    /*
     *  Copy the old solution into the new solution as an initial guess
     */
    vcs_dcopy(VCS_DATA_PTR(m_feSpecies_new), VCS_DATA_PTR(m_feSpecies_old), m_numSpeciesRdc);
    vcs_dcopy(VCS_DATA_PTR(m_actCoeffSpecies_new), VCS_DATA_PTR(m_actCoeffSpecies_old), m_numSpeciesRdc);
    vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_new), VCS_DATA_PTR(m_deltaGRxn_old), m_numRxnRdc);
    vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_Deficient), VCS_DATA_PTR(m_deltaGRxn_old), m_numRxnRdc);
    vcs_dcopy(VCS_DATA_PTR(m_tPhaseMoles_new), VCS_DATA_PTR(m_tPhaseMoles_old), m_numPhases);

    /*
     *     Zero out the entire vector of updates. We sometimes would
     *     query these values below, and we want to be sure that no
     *     information is left from previous iterations.
     */
    vcs_dzero(VCS_DATA_PTR(m_deltaMolNumSpecies), m_numSpeciesTot);

    /*************************************************************************/
    /************** DETERMINE IF DEAD PHASES POP INTO EXISTENCE **************/
    /*************************************************************************/
    /*
     * First step is a major branch in the algorithm.
     * We first determine if a phase pops into existence.
     */
    phasePopPhaseIDs.clear();
    iphasePop = vcs_popPhaseID(phasePopPhaseIDs);
    /*
     *
     */
    soldel = -1;
    if (iphasePop != npos) {
        soldel = vcs_popPhaseRxnStepSizes(iphasePop);
        if (soldel == 3) {
            iphasePop = npos;
#ifdef DEBUG_MODE
            if (m_debug_print_lvl >= 2) {
                plogf("   --- vcs_popPhaseRxnStepSizes() was called but stoich "
                      "prevented phase %d popping\n");
            }
#endif
        }
    }

    /*************************************************************************/
    /* DETERMINE THE REACTION STEP SIZES FOR MAIN STEP AND IF PHASES DIE *****/
    /*************************************************************************/
    /*
     * Don't do this step if there is a phase pop
     */
    iphaseDelete = npos;
    if (iphasePop == npos) {
        /*
         * Figure out the new reaction step sizes
         * for the major species (do minor species in the future too)
         */
        kspec = npos;
        iphaseDelete = vcs_RxnStepSizes(forceComponentCalc, kspec);
    }
#ifdef DEBUG_MODE
    else {
        if (m_debug_print_lvl >= 2) {
            plogf("   --- vcs_RxnStepSizes not called because alternative"
                  "phase creation delta was used instead\n");
        }
    }
#endif
    lec = false;
    doPhaseDeleteIph = npos;
    doPhaseDeleteKspec = npos;
    /*
     *    Zero out the net change in moles of multispecies phases
     */
    vcs_dzero(VCS_DATA_PTR(m_deltaPhaseMoles), m_numPhases);

    /*
     * Check on too many iterations.
     *         If we have too many iterations, Clean up and exit code even though we haven't
     *         converged. -> we have run out of iterations!
     */
    if (m_VCount->Its > maxit) {
        solveFail = -1;
        goto L_RETURN_BLOCK;
    }

    /* ********************************************************************** */
    /* ***************** MAIN LOOP IN CALCULATION *************************** */
    /* ***************** LOOP OVER IRXN TO DETERMINE STEP SIZE ************** */
    /* ********************************************************************** */
    /*
     *   Loop through all of the reactions, irxn, pertaining to the
     *   formation reaction for species kspec in canonical form.
     *
     *   At the end of this loop, we will have a new estimate for the
     *   mole numbers for all species consistent with an extent
     *   of reaction for all noncomponent species formation
     *   reactions. We will have also ensured that all predicted
     *   non-component mole numbers are greater than zero.
     *
     *         Old_Solution               New_Solution             Description
     * -----------------------------------------------------------------------------
     *  m_molNumSpecies_old[kspec]   m_molNumSpecies_new[kspec]   Species Mole Numbers
     *                               m_deltaMolNumSpecies[kspec]  Delta in the Species Mole Numbers
     *
     *
     *
     */
    if (iphaseDelete != npos) {
#ifdef DEBUG_MODE
        if (m_debug_print_lvl >= 2) {
            plogf("   --- Main Loop Treatment -> Circumvented due to Phase Deletion ");
            plogendl();
        }
#endif

        for (k = 0; k < m_numSpeciesTot; k++) {
            m_molNumSpecies_new[k] = m_molNumSpecies_old[k] +  m_deltaMolNumSpecies[k];
            iph = m_phaseID[k];
            m_tPhaseMoles_new[iph] += m_deltaMolNumSpecies[k];
        }
        if (kspec >= m_numComponents) {
            if (m_molNumSpecies_new[k] != 0.0) {
                printf("vcs_solve_tp:: we shouldn't be here!\n");
                exit(EXIT_FAILURE);
            }
            if (m_SSPhase[kspec] == 1) {
                m_speciesStatus[kspec] = VCS_SPECIES_ZEROEDSS;
            } else {
                printf("vcs_solve_tp:: we shouldn't be here!\n");
                exit(EXIT_FAILURE);
            }
            ++m_numRxnMinorZeroed;
            allMinorZeroedSpecies = (m_numRxnMinorZeroed == m_numRxnRdc);
        }
        /*
         *        Set the flags indicating the mole numbers in the vcs_VolPhase
         *        objects are out of date.
         */
        vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW);

        /*
         *         Calculate the new chemical potentials using the tentative
         *         solution values. We only calculate a subset of these, because
         *         we have only updated a subset of the W().
         */
        vcs_dfe(VCS_STATECALC_NEW, 0, 0, m_numSpeciesTot);

        /*
         *         Evaluate DeltaG for all components if ITI=0, and for
         *         major components only if ITI NE 0
         */
        vcs_deltag(0, false, VCS_STATECALC_NEW);
    } else {
#ifdef DEBUG_MODE
        if (m_debug_print_lvl >= 2) {
            plogf("   --- Main Loop Treatment of each non-component species ");
            if (iti == 0) {
                plogf("- Full Calculation:\n");
            } else {
                plogf("- Major Components Calculation:\n");
            }
            plogf("   --- Species     IC    ");
            plogf(" KMoles  Tent_KMoles Rxn_Adj   |    Comment \n");
        }
#endif
        for (irxn = 0; irxn < m_numRxnRdc; irxn++) {
            kspec = m_indexRxnToSpecies[irxn];
            sc_irxn = m_stoichCoeffRxnMatrix[irxn];
            iph = m_phaseID[kspec];
            Vphase = m_VolPhaseList[iph];
#ifdef DEBUG_MODE
            ANOTE[0] = '\0';
#endif
            if (iphasePop != npos) {
                if (iph == iphasePop) {
                    dx =  m_deltaMolNumSpecies[kspec];
                    m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] +  m_deltaMolNumSpecies[kspec];
#ifdef DEBUG_MODE
                    sprintf(ANOTE, "Phase pop");
#endif
                } else {
                    dx = 0.0;
                    m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
                }
            } else {


                if (m_speciesStatus[kspec] == VCS_SPECIES_INTERFACIALVOLTAGE) {
                    /********************************************************************/
                    /************************ VOLTAGE SPECIES ***************************/
                    /********************************************************************/
                    bool soldel_ret;
#ifdef DEBUG_MODE
                    dx = vcs_minor_alt_calc(kspec, irxn, &soldel_ret, ANOTE);
#else
                    dx = vcs_minor_alt_calc(kspec, irxn, &soldel_ret);
#endif
                    soldel = soldel_ret;
                    m_deltaMolNumSpecies[kspec] = dx;
                } else if (m_speciesStatus[kspec] < VCS_SPECIES_MINOR) {
                    /********************************************************************/
                    /********************** ZEROED OUT SPECIES **************************/
                    /********************************************************************/
                    bool resurrect = (m_deltaMolNumSpecies[kspec] > 0.0);
#ifdef DEBUG_MODE
                    if (m_debug_print_lvl >= 3) {
                        plogf("   --- %s currently zeroed (SpStatus=%-2d):",
                              m_speciesName[kspec].c_str(), m_speciesStatus[kspec]);
                        plogf("%3d DG = %11.4E WT = %11.4E W = %11.4E DS = %11.4E\n",
                              irxn, m_deltaGRxn_new[irxn], m_molNumSpecies_new[kspec],
                              m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec]);
                    }
#endif

                    if (m_deltaGRxn_new[irxn] >= 0.0 || !resurrect) {
                        m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
                        m_deltaMolNumSpecies[kspec] = 0.0;
                        resurrect = false;
#ifdef DEBUG_MODE
                        sprintf(ANOTE, "Species stays zeroed: DG = %11.4E", m_deltaGRxn_new[irxn]);
                        if (m_deltaGRxn_new[irxn] < 0.0) {
                            if (m_speciesStatus[kspec] == VCS_SPECIES_STOICHZERO) {
                                sprintf(ANOTE, "Species stays zeroed even though dg neg due to "
                                        "STOICH/PHASEPOP constraint: DG = %11.4E",
                                        m_deltaGRxn_new[irxn]);
                            } else {
                                sprintf(ANOTE, "Species stays zeroed even though dg neg: DG = %11.4E, ds zeroed",
                                        m_deltaGRxn_new[irxn]);
                            }
                        }
#endif
                    } else {
                        for (size_t j = 0; j < m_numElemConstraints; ++j) {
                            int elType = m_elType[j];
                            if (elType == VCS_ELEM_TYPE_ABSPOS) {
                                atomComp = m_formulaMatrix[j][kspec];
                                if (atomComp > 0.0) {
                                    double maxPermissible = m_elemAbundancesGoal[j] / atomComp;
                                    if (maxPermissible < VCS_DELETE_MINORSPECIES_CUTOFF) {
#ifdef DEBUG_MODE
                                        sprintf(ANOTE, "Species stays zeroed even though dG "
                                                "neg, because of %s elemAbund",
                                                m_elementName[j].c_str());
#endif
                                        resurrect = false;
                                        break;
                                    }
                                }
                            }
                        }
                    }
                    /*
                     * Resurrect the species
                     */
                    if (resurrect) {
                        bool phaseResurrected = false;
                        if (Vphase->exists() == VCS_PHASE_EXIST_NO) {
                            //Vphase->setExistence(1);
                            phaseResurrected = true;
                        }

                        if (phaseResurrected) {
#ifdef DEBUG_MODE
                            if (m_debug_print_lvl >= 2) {
                                plogf("   --- Zeroed species changed to major: ");
                                plogf("%-12s\n", m_speciesName[kspec].c_str());
                            }
#endif
                            m_speciesStatus[kspec] = VCS_SPECIES_MAJOR;
                            allMinorZeroedSpecies = false;
                        } else {
#ifdef DEBUG_MODE
                            if (m_debug_print_lvl >= 2) {
                                plogf("   --- Zeroed species changed to minor: ");
                                plogf("%-12s\n", m_speciesName[kspec].c_str());
                            }
#endif
                            m_speciesStatus[kspec] = VCS_SPECIES_MINOR;
                        }
                        if (m_deltaMolNumSpecies[kspec] > 0.0) {
                            dx = m_deltaMolNumSpecies[kspec] * 0.01;
                            m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
                        } else {
                            m_molNumSpecies_new[kspec] = m_totalMolNum * VCS_DELETE_PHASE_CUTOFF * 10.;
                            dx = m_molNumSpecies_new[kspec] - m_molNumSpecies_old[kspec];
                        }
                        m_deltaMolNumSpecies[kspec] = dx;
#ifdef DEBUG_MODE
                        sprintf(ANOTE, "Born:IC=-1 to IC=1:DG=%11.4E", m_deltaGRxn_new[irxn]);
#endif
                    } else {
                        m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
                        m_deltaMolNumSpecies[kspec] = 0.0;
                        dx = 0.0;
                    }
                } else if (m_speciesStatus[kspec] == VCS_SPECIES_MINOR) {
                    /********************************************************************/
                    /***************************** MINOR SPECIES ************************/
                    /********************************************************************/
                    /*
                     *    Unless ITI isn't equal to zero we zero out changes
                     *    to minor species.
                     */
                    if (iti != 0) {
                        m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
                        m_deltaMolNumSpecies[kspec] = 0.0;
                        dx = 0.0;
#ifdef DEBUG_MODE
                        sprintf(ANOTE,"minor species not considered");
                        if (m_debug_print_lvl >= 2) {
                            plogf("   --- ");
                            plogf("%-12s", m_speciesName[kspec].c_str());
                            plogf("%3d%11.4E%11.4E%11.4E | %s",
                                  m_speciesStatus[kspec], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec],
                                  m_deltaMolNumSpecies[kspec], ANOTE);
                            plogendl();
                        }
#endif
                        continue;
                    }
                    /*
                     *        Minor species alternative calculation
                     *       ---------------------------------------
                     *    This is based upon the following approximation:
                     *    The mole fraction changes due to these reactions don't affect
                     *    the mole numbers of the component species. Therefore the
                     *    following approximation is valid for an ideal solution
                     *       0 = DG(I) + log(WT(I)/W(I))
                     *       (DG contains the contribution from FF(I) + log(W(I)/TL) )
                     *    Thus,
                     *        WT(I) = W(I) EXP(-DG(I))
                     *    If soldel is true on return, then we branch to the section
                     *    that deletes a species from the current set of active species.
                     */
                    bool soldel_ret;
#ifdef DEBUG_MODE
                    dx = vcs_minor_alt_calc(kspec, irxn, &soldel_ret, ANOTE);
#else
                    dx = vcs_minor_alt_calc(kspec, irxn, &soldel_ret);
#endif
                    soldel = soldel_ret;
                    m_deltaMolNumSpecies[kspec] = dx;
                    m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;

                    if (soldel) {
                        /*******************************************************************/
                        /*****  DELETE MINOR SPECIES LESS THAN  VCS_DELETE_SPECIES_CUTOFF  */
                        /*****  MOLE NUMBER                                                */
                        /*******************************************************************/
#ifdef DEBUG_MODE
                        if (m_debug_print_lvl >= 2) {
                            plogf("   --- Delete minor species in multispec phase: %-12s",
                                  m_speciesName[kspec].c_str());
                            plogendl();
                        }
#endif
                        m_deltaMolNumSpecies[kspec] = 0.0;
                        /*
                         *       Delete species, kspec. The alternate return is for the case
                         *       where all species become deleted. Then, we need to
                         *       branch to the code where we reevaluate the deletion
                         *       of all species.
                         */
                        lnospec = vcs_delete_species(kspec);
                        if (lnospec) {
                            goto L_RECHECK_DELETED;
                        }
                        /*
                         *       Go back to consider the next species in the list.
                         *       Note, however, that the next species in the list is now
                         *       in slot l. In deleting the previous species L, We have
                         *       exchanged slot MR with slot l, and then have
                         *       decremented MR.
                         *       Therefore, we will decrement the species counter, here.
                         */
                        --irxn;
#ifdef DEBUG_MODE
                        goto L_MAIN_LOOP_END_NO_PRINT;
#else
                        goto L_MAIN_LOOP_END;
#endif
                    }
                } else {
                    /********************************************************************/
                    /*********************** MAJOR SPECIES ******************************/
                    /********************************************************************/
#ifdef DEBUG_MODE
                    sprintf(ANOTE, "Normal Major Calc");
#endif
                    /*
                     * Check for superconvergence of the formation reaction. Do
                     * nothing if it is superconverged. Skip to the end of the
                     * irxn loop if it is superconverged.
                     */
                    if (fabs(m_deltaGRxn_new[irxn]) <= m_tolmaj2) {
                        m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
                        m_deltaMolNumSpecies[kspec] = 0.0;
                        dx = 0.0;
#ifdef DEBUG_MODE
                        sprintf(ANOTE, "major species is converged");
                        if (m_debug_print_lvl >= 2) {
                            plogf("   --- ");
                            plogf("%-12s", m_speciesName[kspec].c_str());
                            plogf("%3d%11.4E%11.4E%11.4E | %s",
                                  m_speciesStatus[kspec], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec],
                                  m_deltaMolNumSpecies[kspec], ANOTE);
                            plogendl();
                        }
#endif
                        continue;
                    }
                    /*
                     *      Set the initial step size, dx, equal to the value produced
                     *      by the routine, vcs_RxnStepSize().
                     *
                     *          Note the multiplication logic is to make sure that
                     *          dg[] didn't change sign due to w[] changing in the
                     *          middle of the iteration. (it can if a single species
                     *          phase goes out of existence).
                     */
                    if ((m_deltaGRxn_new[irxn] * m_deltaMolNumSpecies[kspec]) <= 0.0) {
                        dx = m_deltaMolNumSpecies[kspec];
                    } else {
                        dx = 0.0;
                        m_deltaMolNumSpecies[kspec] = 0.0;
#ifdef DEBUG_MODE
                        sprintf(ANOTE, "dx set to 0, DG flipped sign due to "
                                "changed initial point");
#endif
                    }
                    /*
                     *      Form a tentative value of the new species moles
                     */
                    m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;

                    /*
                     *      Check for non-positive mole fraction of major species.
                     *      If we find one, we branch to a section below. Then,
                     *      depending upon the outcome, we branch to sections below,
                     *      or we restart the entire iteration.
                     */
                    if (m_molNumSpecies_new[kspec] <= 0.0) {
#ifdef DEBUG_MODE
                        sprintf(ANOTE, "initial nonpos kmoles= %11.3E",
                                m_molNumSpecies_new[kspec]);
#endif
                        /* ************************************************* */
                        /* *** NON-POSITIVE MOLES OF MAJOR SPECIES ********* */
                        /* ************************************************* */
                        /*
                         *          We are here when a tentative value of a mole fraction
                         *          created by a tentative value of M_DELTAMOLNUMSPECIES(*) is negative.
                         *          We branch from here depending upon whether this
                         *          species is in a single species phase or in
                         *          a multispecies phase.
                         */
                        if (!(m_SSPhase[kspec])) {
                            /*
                             *   Section for multispecies phases:
                             *     - Cut reaction adjustment for positive kmoles of
                             *       major species in multispecies phases.
                             *       Decrease its concentration by a factor of 10.
                             */
                            dx = -0.9 * m_molNumSpecies_old[kspec];
                            m_deltaMolNumSpecies[kspec] = dx;
                            m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
                        } else {
                            /*
                             *   Section for single species phases:
                             *       Calculate a dx that will wipe out the
                             *       moles in the phase.
                             */
                            dx = -m_molNumSpecies_old[kspec];
                            /*
                             *       Calculate an update that doesn't create a negative mole
                             *       number for a component species. Actually, restrict this
                             *       a little more so that the component values can only be
                             *       reduced by two 99%,
                             */
                            for (j = 0; j < m_numComponents; ++j) {
                                if (sc_irxn[j] != 0.0) {
                                    wx[j] = m_molNumSpecies_old[j] + sc_irxn[j] * dx;
                                    if (wx[j] <= m_molNumSpecies_old[j] * 0.01 - 1.0E-150) {
                                        dx = std::max(dx,  m_molNumSpecies_old[j] * -0.99 / sc_irxn[j]);
                                    }
                                } else {
                                    wx[j] = m_molNumSpecies_old[j];
                                }
                            }
                            m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
                            if (m_molNumSpecies_new[kspec] > 0.0) {
                                m_deltaMolNumSpecies[kspec] = dx;
#ifdef DEBUG_MODE
                                sprintf(ANOTE,
                                        "zeroing SS phase created a neg component species "
                                        "-> reducing step size instead");
#endif
                            } else {
                                /*
                                 *     We are going to zero the single species phase.
                                 *     Set the existence flag
                                 */
                                iph = m_phaseID[kspec];
                                Vphase = m_VolPhaseList[iph];
                                //Vphase->setExistence(0);
#ifdef DEBUG_MODE
                                sprintf(ANOTE, "zeroing out SS phase: ");
#endif
                                /*
                                 *     Change the base mole numbers for the iteration.
                                 *     We need to do this here, because we have decided
                                 *     to eliminate  the phase in this special section
                                 *     outside the main loop.
                                 */
                                m_molNumSpecies_new[kspec] = 0.0;
                                doPhaseDeleteIph = iph;
                                doPhaseDeleteKspec = kspec;

#ifdef DEBUG_MODE
                                if (m_debug_print_lvl >= 2) {
                                    if (m_speciesStatus[kspec] >= 0) {
                                        plogf("   --- SS species changed to zeroedss: ");
                                        plogf("%-12s", m_speciesName[kspec].c_str());
                                        plogendl();
                                    }
                                }
#endif
                                m_speciesStatus[kspec] = VCS_SPECIES_ZEROEDSS;
                                ++m_numRxnMinorZeroed;
                                allMinorZeroedSpecies = (m_numRxnMinorZeroed == m_numRxnRdc);

                                for (size_t kk = 0; kk < m_numSpeciesTot; kk++) {
                                    m_deltaMolNumSpecies[kk] = 0.0;
                                    m_molNumSpecies_new[kk] = m_molNumSpecies_old[kk];
                                }
                                m_deltaMolNumSpecies[kspec] = dx;
                                m_molNumSpecies_new[kspec] = 0.0;

                                for (k = 0; k < m_numComponents; ++k) {
                                    m_deltaMolNumSpecies[k] = 0.0;
                                }
                                for (iph = 0; iph < m_numPhases; iph++) {
                                    m_deltaPhaseMoles[iph] = 0.0;
                                }

                            }
                        }

                    }

#ifdef VCS_LINE_SEARCH
                    /*********************************************************************/
                    /*** LINE SEARCH ALGORITHM FOR MAJOR SPECIES IN NON-IDEAL PHASES *****/
                    /*********************************************************************/
                    /*
                     * Skip the line search if we are birthing a species
                     */
                    if ((dx != 0.0) &&
                            (m_molNumSpecies_old[kspec] > 0.0) &&
                            (doPhaseDeleteIph == -1) &&
                            (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE)) {
                        double dx_old = dx;

#ifdef DEBUG_MODE
                        dx = vcs_line_search(irxn, dx_old, ANOTE);
#else
                        dx = vcs_line_search(irxn, dx_old);
#endif
                        vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW);
                    }
                    m_deltaMolNumSpecies[kspec] = dx;
#endif
                }/* End of Loop on ic[irxn] -> the type of species */
            }
            /***********************************************************************/
            /****** CALCULATE KMOLE NUMBER CHANGE FOR THE COMPONENT BASIS **********/
            /***********************************************************************/
            if (dx != 0.0 && (m_speciesUnknownType[kspec] !=
                              VCS_SPECIES_TYPE_INTERFACIALVOLTAGE)) {
                /*
                 *         Change the amount of the component compounds according
                 *         to the reaction delta that we just computed.
                 *         This should keep the amount of material constant.
                 */
#ifdef DEBUG_MODE
                if (fabs(m_deltaMolNumSpecies[kspec] -dx) >
                        1.0E-14*(fabs(m_deltaMolNumSpecies[kspec]) + fabs(dx) + 1.0E-32)) {
                    plogf(" ds[kspec] = %20.16g dx = %20.16g , kspec = %d\n",
                          m_deltaMolNumSpecies[kspec], dx, kspec);
                    plogf("we have a problem!");
                    plogendl();
                    exit(EXIT_FAILURE);
                }
#endif
                for (k = 0; k < m_numComponents; ++k) {
                    m_deltaMolNumSpecies[k] += sc_irxn[k] * dx;
                }
                /*
                 *         Calculate the tentative change in the total number of
                 *         moles in all of the phases
                 */
                dnPhase_irxn = m_deltaMolNumPhase[irxn];
                for (iph = 0; iph < m_numPhases; iph++) {
                    m_deltaPhaseMoles[iph] += dx * dnPhase_irxn[iph];
                }
            }

#ifdef DEBUG_MODE
            checkDelta1(VCS_DATA_PTR(m_deltaMolNumSpecies),
                        VCS_DATA_PTR(m_deltaPhaseMoles), kspec+1);
#endif
            /*
             *          Branch point for returning -
             */
#ifndef DEBUG_MODE
L_MAIN_LOOP_END:
            ;
#endif
#ifdef DEBUG_MODE
            if (m_debug_print_lvl >= 2) {
                m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + m_deltaMolNumSpecies[kspec];
                plogf("   --- ");
                plogf("%-12.12s", m_speciesName[kspec].c_str());
                plogf("%3d%11.4E%11.4E%11.4E | %s",
                      m_speciesStatus[kspec], m_molNumSpecies_old[kspec],
                      m_molNumSpecies_new[kspec],
                      m_deltaMolNumSpecies[kspec], ANOTE);
                plogendl();
            }
L_MAIN_LOOP_END_NO_PRINT:
            ;
#endif
            if (doPhaseDeleteIph != npos) {
#ifdef DEBUG_MODE
                if (m_debug_print_lvl >= 2) {
                    plogf("   --- ");
                    plogf("%-12.12s Main Loop Special Case deleting phase with species: ",
                          m_speciesName[doPhaseDeleteKspec].c_str());
                    plogendl();
                }
#endif
                break;
            }
        }  /**************** END OF MAIN LOOP OVER FORMATION REACTIONS ************/

#ifdef DEBUG_MODE
        if (m_debug_print_lvl >= 2) {
            for (k = 0; k < m_numComponents; k++) {
                plogf("   --- ");
                plogf("%-12.12s", m_speciesName[k].c_str());
                plogf("  c%11.4E%11.4E%11.4E |\n",
                      m_molNumSpecies_old[k],
                      m_molNumSpecies_old[k]+m_deltaMolNumSpecies[k], m_deltaMolNumSpecies[k]);
            }
            plogf("   ");
            vcs_print_line("-", 80);
            plogf("   --- Finished Main Loop");
            plogendl();
        }
#endif

        /*************************************************************************/
        /*********** LIMIT REDUCTION OF BASIS SPECIES TO 99% *********************/
        /*************************************************************************/
        /*
         *        We have a tentative m_deltaMolNumSpecies[]. Now apply other criteria
         *        to limit its magnitude.
         *
         *
         */
        par = 0.5;
        for (k = 0; k < m_numComponents; ++k) {
            if (m_molNumSpecies_old[k] > 0.0) {
                xx = -m_deltaMolNumSpecies[k] / m_molNumSpecies_old[k];
                if (par < xx) {
                    par = xx;
#ifdef DEBUG_MODE
                    ll = k;
#endif
                }
            } else {
                if (m_deltaMolNumSpecies[k] < 0.0) {
                    /*
                     * If we are here, we then do a step which violates element
                     * conservation.
                     */
                    iph = m_phaseID[k];
                    m_deltaPhaseMoles[iph] -= m_deltaMolNumSpecies[k];
                    m_deltaMolNumSpecies[k] = 0.0;
                }
            }
        }
        par = 1.0 / par;
        if (par <= 1.01 && par > 0.0) {
            /* Reduce the size of the step by the multiplicative factor, par */
            par *= 0.99;
#ifdef DEBUG_MODE
            if (m_debug_print_lvl >= 2) {
                plogf("   --- Reduction in step size due to component ");
                plogf("%s", m_speciesName[ll].c_str());
                plogf(" going negative = %11.3E", par);
                plogendl();
            }
#endif
            for (i = 0; i < m_numSpeciesTot; ++i) {
                m_deltaMolNumSpecies[i] *= par;
            }
            for (iph = 0; iph < m_numPhases; iph++) {
                m_deltaPhaseMoles[iph] *= par;
            }
        } else {
            par = 1.0;
        }
#ifdef DEBUG_MODE
        checkDelta1(VCS_DATA_PTR(m_deltaMolNumSpecies),
                    VCS_DATA_PTR(m_deltaPhaseMoles), m_numSpeciesTot);
#endif

        /*
         *      Now adjust the wt[kspec]'s so that the reflect the decrease in
         *      the overall length of m_deltaMolNumSpecies[kspec] just calculated. At the end
         *      of this section wt[], m_deltaMolNumSpecies[], tPhMoles, and tPhMoles1 should all be
         *      consistent with a new estimate of the state of the system.
         */
        for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
            m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + m_deltaMolNumSpecies[kspec];
            if (m_molNumSpecies_new[kspec] < 0.0 && (m_speciesUnknownType[kspec]
                    != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE)) {
                plogf("vcs_solve_TP: ERROR on step change wt[%d:%s]: %g < 0.0",
                      kspec, m_speciesName[kspec].c_str(), m_molNumSpecies_new[kspec]);
                plogendl();
                exit(EXIT_FAILURE);
            }
        }

        /*
         *        Calculate the tentative total mole numbers for each phase
         */
        for (iph = 0; iph < m_numPhases; iph++) {
            m_tPhaseMoles_new[iph] = m_tPhaseMoles_old[iph] + m_deltaPhaseMoles[iph];
        }

        /*
         *        Set the flags indicating the mole numbers in the vcs_VolPhase
         *        objects are out of date.
         */
        vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW);

        /*
         *         Calculate the new chemical potentials using the tentative
         *         solution values. We only calculate a subset of these, because
         *         we have only updated a subset of the W().
         */
        vcs_dfe(VCS_STATECALC_NEW, 0, 0, m_numSpeciesTot);

        /*
         *         Evaluate DeltaG for all components if ITI=0, and for
         *         major components only if ITI NE 0
         */
        vcs_deltag(0, false, VCS_STATECALC_NEW);

        /* *************************************************************** */
        /* **** CONVERGENCE FORCER SECTION ******************************* */
        /* *************************************************************** */
        if (printDetails) {
            plogf("   --- Total Old       Dimensionless Gibbs Free Energy = %20.13E\n",
                  vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_old),
                                  VCS_DATA_PTR(m_tPhaseMoles_old)));
            plogf("   --- Total tentative Dimensionless Gibbs Free Energy = %20.13E",
                  vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_feSpecies_new),
                                  VCS_DATA_PTR(m_tPhaseMoles_new)));
            plogendl();
        }

        forced = vcs_globStepDamp();

        /*
         *       Print out the changes to the solution that FORCER produced
         */
        if (printDetails && forced) {

            plogf(" -----------------------------------------------------\n");
            plogf("   --- FORCER SUBROUTINE changed the solution:\n");
            plogf("   --- SPECIES Status INIT MOLES TENT_MOLES");
            plogf("  FINAL KMOLES  INIT_DEL_G/RT  TENT_DEL_G/RT  FINAL_DELTA_G/RT\n");
            for (i = 0; i < m_numComponents; ++i) {
                plogf("  --- %-12.12s", m_speciesName[i].c_str());
                plogf("    %14.6E %14.6E %14.6E\n",  m_molNumSpecies_old[i],
                      m_molNumSpecies_old[i] + m_deltaMolNumSpecies[i], m_molNumSpecies_new[i]);
            }
            for (kspec = m_numComponents; kspec < m_numSpeciesRdc; ++kspec) {
                irxn = kspec - m_numComponents;
                plogf("  --- %-12.12s", m_speciesName[kspec].c_str());
                plogf(" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n", m_speciesStatus[kspec],
                      m_molNumSpecies_old[kspec],
                      m_molNumSpecies_old[kspec]+m_deltaMolNumSpecies[kspec],
                      m_molNumSpecies_new[kspec], m_deltaGRxn_old[irxn],
                      m_deltaGRxn_tmp[irxn], m_deltaGRxn_new[irxn]);
            }
            print_space(26);
            plogf("Norms of Delta G():%14.6E%14.6E\n",
                  l2normdg(VCS_DATA_PTR(m_deltaGRxn_old)),
                  l2normdg(VCS_DATA_PTR(m_deltaGRxn_new)));
            plogf("   Total kmoles of gas    = %15.7E\n", m_tPhaseMoles_old[0]);
            if ((m_numPhases > 1) && (!(m_VolPhaseList[1])->m_singleSpecies)) {
                plogf("   Total kmoles of liquid = %15.7E\n", m_tPhaseMoles_old[1]);
            } else {
                plogf("   Total kmoles of liquid = %15.7E\n", 0.0);
            }
            plogf("   Total New Dimensionless Gibbs Free Energy = %20.13E\n",
                  vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_feSpecies_new),
                                  VCS_DATA_PTR(m_tPhaseMoles_new)));
            plogf(" -----------------------------------------------------");
            plogendl();
        }
    }
    /* *************************************************************** */
    /* **** ITERATION SUMMARY PRINTOUT SECTION *********************** */
    /* *************************************************************** */

    if (printDetails) {
        plogf("   ");
        vcs_print_line("-", 103);
        plogf("   --- Summary of the Update ");
        if (iti == 0) {
            plogf(" (all species):");
        } else {
            plogf(" (only major species):");
        }
        if (m_totalMoleScale != 1.0) {
            plogf(" (Total Mole Scale = %g)", m_totalMoleScale);
        }
        plogf("\n");
        plogf("   ---      Species Status Initial_KMoles Final_KMoles Initial_Mu/RT");
        plogf("     Mu/RT     Init_Del_G/RT   Delta_G/RT\n");
        for (i = 0; i < m_numComponents; ++i) {
            plogf("   ---   %-12.12s", m_speciesName[i].c_str());
            plogf("    ");
            plogf("%14.6E%14.6E%14.6E%14.6E\n", m_molNumSpecies_old[i],
                  m_molNumSpecies_new[i], m_feSpecies_old[i], m_feSpecies_new[i]);
        }
        for (i = m_numComponents; i < m_numSpeciesRdc; ++i) {
            l1 = i - m_numComponents;
            plogf("   ---   %-12.12s", m_speciesName[i].c_str());
            plogf(" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n",
                  m_speciesStatus[i], m_molNumSpecies_old[i],
                  m_molNumSpecies_new[i], m_feSpecies_old[i], m_feSpecies_new[i],
                  m_deltaGRxn_old[l1], m_deltaGRxn_new[l1]);
        }
        for (kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) {
            l1 = kspec - m_numComponents;
            plogf("   ---   %-12.12s", m_speciesName[kspec].c_str());
            plogf(" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n",
                  m_speciesStatus[kspec], m_molNumSpecies_old[kspec],
                  m_molNumSpecies_new[kspec], m_feSpecies_old[kspec], m_feSpecies_new[kspec],
                  m_deltaGRxn_old[l1], m_deltaGRxn_new[l1]);
        }
        plogf("   ---");
        print_space(56);
        plogf("Norms of Delta G():%14.6E%14.6E",
              l2normdg(VCS_DATA_PTR(m_deltaGRxn_old)),
              l2normdg(VCS_DATA_PTR(m_deltaGRxn_new)));
        plogendl();

        plogf("   ---           Phase_Name    KMoles(after update)\n");
        plogf("   ---   ");
        vcs_print_line("-", 50);
        for (iph = 0; iph < m_numPhases; iph++) {
            Vphase = m_VolPhaseList[iph];
            plogf("   ---   %18s = %15.7E\n", Vphase->PhaseName.c_str(), m_tPhaseMoles_new[iph]);
        }
        plogf("   ");
        vcs_print_line("-", 103);
        plogf("   --- Total Old Dimensionless Gibbs Free Energy = %20.13E\n",
              vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_old),
                              VCS_DATA_PTR(m_tPhaseMoles_old)));
        plogf("   --- Total New Dimensionless Gibbs Free Energy = %20.13E",
              vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_feSpecies_new),
                              VCS_DATA_PTR(m_tPhaseMoles_new)));
        plogendl();
#ifdef DEBUG_MODE
        if (m_VCount->Its > 550) {
            plogf("   --- Troublesome solve");
            plogendl();
        }
#endif
    }

    /*************************************************************************/
    /******************* RESET VALUES AT END OF ITERATION ********************/
    /******************* UPDATE MOLE NUMBERS *********************************/
    /*************************************************************************/
    /*
     *       If the solution wasn't changed in the forcer routine,
     *       then copy the tentative mole numbers and Phase moles
     *       into the actual mole numbers and phase moles.
     *       We will consider this current step to be completed.
     *
     *   Accept the step. -> the tentative solution now becomes
     *                the real solution. If FORCED is true, then
     *                we have already done this inside the FORCED
     *                loop.
     */
    vcs_updateMolNumVolPhases(VCS_STATECALC_NEW);
    vcs_dcopy(VCS_DATA_PTR(m_tPhaseMoles_old), VCS_DATA_PTR(m_tPhaseMoles_new), m_numPhases);
    vcs_dcopy(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_molNumSpecies_new),
              m_numSpeciesRdc);
    vcs_dcopy(VCS_DATA_PTR(m_actCoeffSpecies_old),
              VCS_DATA_PTR(m_actCoeffSpecies_new), m_numSpeciesRdc);
    vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_old), VCS_DATA_PTR(m_deltaGRxn_new), m_numRxnRdc);
    vcs_dcopy(VCS_DATA_PTR(m_feSpecies_old), VCS_DATA_PTR(m_feSpecies_new), m_numSpeciesRdc);

    vcs_setFlagsVolPhases(true, VCS_STATECALC_OLD);
    /*
     *       Increment the iteration counters
     */
    ++(m_VCount->Its);
    ++it1;
#ifdef DEBUG_MODE
    if (m_debug_print_lvl >= 2) {
        plogf("   --- Increment counter increased, step is accepted: %4d",
              m_VCount->Its);
        plogendl();
    }
#endif
    /*************************************************************************/
    /******************* HANDLE DELETION OF MULTISPECIES PHASES **************/
    /*************************************************************************/
    /*
     *   We delete multiphases, when the total moles in the multiphase
     *   is reduced below a relative threshold.
     *   Set microscopic multispecies phases with total relative
     *   number of moles less than  VCS_DELETE_PHASE_CUTOFF to
     *   absolute zero.
     */
    justDeletedMultiPhase = false;
    for (iph = 0; iph < m_numPhases; iph++) {
        Vphase = m_VolPhaseList[iph];
        if (!(Vphase->m_singleSpecies)) {
            if (m_tPhaseMoles_old[iph] != 0.0 &&
                    m_tPhaseMoles_old[iph]/m_totalMolNum <= VCS_DELETE_PHASE_CUTOFF) {
                soldel = 1;

                for (kspec = 0; kspec < m_numSpeciesRdc; kspec++) {
                    if (m_phaseID[kspec] == iph && m_molNumSpecies_old[kspec] > 0.0) {
                        irxn = kspec - m_numComponents;
                    }
                }
                if (soldel) {
#ifdef DEBUG_MODE
                    if (m_debug_print_lvl >= 1) {
                        plogf("   --- Setting microscopic phase %d to zero", iph);
                        plogendl();
                    }
#endif
                    justDeletedMultiPhase = true;
                    vcs_delete_multiphase(iph);
                }
            }
        }
    }
    /*
     *       If we have deleted a multispecies phase because the
     *       equilibrium moles decreased, then we will update all
     *       the component basis calculation, and therefore all
     *       of the thermo functions just to be safe.
     */
    if (justDeletedMultiPhase) {
        justDeletedMultiPhase = false;
        retn = vcs_basopt(false, VCS_DATA_PTR(aw), VCS_DATA_PTR(sa),
                          VCS_DATA_PTR(sm), VCS_DATA_PTR(ss), test,
                          &usedZeroedSpecies);
        if (retn != VCS_SUCCESS) {
#ifdef DEBUG_MODE
            plogf("   --- BASOPT returned with an error condition\n");
#endif
            exit(EXIT_FAILURE);
        }
        vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
        vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc);

        vcs_deltag(0, true, VCS_STATECALC_OLD);
        iti = 0;
        goto  L_MAINLOOP_ALL_SPECIES ;

    }
    /*************************************************************************/
    /***************** CHECK FOR ELEMENT ABUNDANCE****************************/
    /*************************************************************************/
#ifdef DEBUG_MODE
    if (m_debug_print_lvl >= 2) {
        plogf("   --- Normal element abundance check");
    }
#endif
    vcs_elab();
    if (! vcs_elabcheck(0)) {
#ifdef DEBUG_MODE
        if (m_debug_print_lvl >= 2) {
            plogf(" - failed -> redoing element abundances.");
            plogendl();
        }
#endif
        vcs_elcorr(VCS_DATA_PTR(sm), VCS_DATA_PTR(wx));
        vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
        vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc);
        vcs_deltag(0, true,  VCS_STATECALC_OLD);
        uptodate_minors = true;
    }
#ifdef DEBUG_MODE
    else {
        if (m_debug_print_lvl >= 2) {
            plogf(" - passed");
            plogendl();
        }
    }
#endif
    /*************************************************************************/
    /***************** CHECK FOR OPTIMUM BASIS *******************************/
    /*************************************************************************/
    /*
     *   HKM -> We first evaluate whether the components species are
     *          ordered according to their mole numbers. If they are,
     *          then we can essential do an order(NR) operation instead
     *          of an order(NR*NC) operation to determine whether
     *          a new basis is needed.
     *
     *   HKM -> This section used to be branched to initially if
     *          there was a machine estimate. I took it out to simplify
     *          the code logic.
     */
    dofast = (m_numComponents != 1);
    for (i = 1; i < m_numComponents; ++i) {
        if ((m_molNumSpecies_old[i - 1] * m_spSize[i-1]) < (m_molNumSpecies_old[i] * m_spSize[i])) {
            dofast = false;
            break;
        }
    }
    dofast = false;
    if (dofast) {
        for (i = 0; i < m_numRxnRdc; ++i) {
            l = m_indexRxnToSpecies[i];
            if (m_speciesUnknownType[l] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
                for (j = m_numComponents - 1; j != npos; j--) {
                    bool doSwap = false;
                    if (m_SSPhase[j]) {
                        doSwap = (m_molNumSpecies_old[l] * m_spSize[l]) >
                                 (m_molNumSpecies_old[j] * m_spSize[j] * 1.01);
                        if (!m_SSPhase[i]) {
                            if (doSwap) {
                                doSwap = (m_molNumSpecies_old[l]) > (m_molNumSpecies_old[j] * 1.01);
                            }
                        }
                    } else {
                        if (m_SSPhase[i]) {
                            doSwap = (m_molNumSpecies_old[l] * m_spSize[l]) >
                                     (m_molNumSpecies_old[j] * m_spSize[j] * 1.01);
                            if (!doSwap) {
                                doSwap = (m_molNumSpecies_old[l]) > (m_molNumSpecies_old[j] * 1.01);
                            }
                        } else {
                            doSwap = (m_molNumSpecies_old[l] * m_spSize[l]) >
                                     (m_molNumSpecies_old[j] * m_spSize[j] * 1.01);
                        }
                    }
                    if (doSwap) {
                        if (m_stoichCoeffRxnMatrix[i][j] != 0.0) {
#ifdef DEBUG_MODE
                            if (m_debug_print_lvl >= 2) {
                                plogf("   --- Get a new basis because %s", m_speciesName[l].c_str());
                                plogf(" is better than comp %s", m_speciesName[j].c_str());
                                plogf(" and share nonzero stoic: %-9.1f",
                                      m_stoichCoeffRxnMatrix[i][j]);
                                plogendl();
                            }
#endif
                            goto L_COMPONENT_CALC;
                        }
                    } else {
                        break;
                    }
#ifdef DEBUG_NOT
                    if (m_speciesStatus[l] == VCS_SPECIES_ZEROEDMS) {
                        if (m_molNumSpecies_old[j] == 0.0) {
                            if (m_stoichCoeffRxnMatrix[i][j] != 0.0) {
                                if (dg[i] < 0.0) {
#ifdef DEBUG_MODE
                                    if (m_debug_print_lvl >= 2) {
                                        plogf("   --- Get a new basis because %s", m_speciesName[l].c_str());
                                        plogf(" has dg < 0.0 and comp %s has zero mole num", m_speciesName[j].c_str());
                                        plogf(" and share nonzero stoic: %-9.1f",
                                              m_stoichCoeffRxnMatrix[i][j]);
                                        plogendl();
                                    }
#endif
                                    goto L_COMPONENT_CALC;
                                }
                            }
                        }
                    }
#endif
                }
            }
        }
    } else {
        for (i = 0; i < m_numRxnRdc; ++i) {
            l = m_indexRxnToSpecies[i];
            if (m_speciesUnknownType[l] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
                for (j = 0; j < m_numComponents; ++j) {
                    bool doSwap = false;
                    if (m_SSPhase[j]) {
                        doSwap = (m_molNumSpecies_old[l] * m_spSize[l]) >
                                 (m_molNumSpecies_old[j] * m_spSize[j] * 1.01);
                        if (!m_SSPhase[l]) {
                            if (doSwap) {
                                doSwap = (m_molNumSpecies_old[l]) > (m_molNumSpecies_old[j] * 1.01);
                            }
                        }
                    } else {
                        if (m_SSPhase[l]) {
                            doSwap = (m_molNumSpecies_old[l] * m_spSize[l]) >
                                     (m_molNumSpecies_old[j] * m_spSize[j] * 1.01);
                            if (!doSwap) {
                                doSwap = (m_molNumSpecies_old[l]) > (m_molNumSpecies_old[j] * 1.01);
                            }
                        } else {
                            doSwap = (m_molNumSpecies_old[l] * m_spSize[l]) >
                                     (m_molNumSpecies_old[j] * m_spSize[j] * 1.01);
                        }
                    }
                    if (doSwap) {
                        if (m_stoichCoeffRxnMatrix[i][j] != 0.0) {
#ifdef DEBUG_MODE
                            if (m_debug_print_lvl >= 2) {
                                plogf("   --- Get a new basis because ");
                                plogf("%s", m_speciesName[l].c_str());
                                plogf(" is better than comp ");
                                plogf("%s", m_speciesName[j].c_str());
                                plogf(" and share nonzero stoic: %-9.1f",
                                      m_stoichCoeffRxnMatrix[i][j]);
                                plogendl();
                            }
#endif
                            goto L_COMPONENT_CALC;
                        }
                    }
#ifdef DEBUG_NOT
                    if (m_speciesStatus[l] == VCS_SPECIES_ZEROEDMS) {
                        if (m_molNumSpecies_old[j] == 0.0) {
                            if (m_stoichCoeffRxnMatrix[i][j] != 0.0) {
                                if (dg[i] < 0.0) {
#ifdef DEBUG_MODE
                                    if (m_debug_print_lvl >= 2) {
                                        plogf("   --- Get a new basis because %s", m_speciesName[l].c_str());
                                        plogf(" has dg < 0.0 and comp %s has zero mole num",
                                              m_speciesName[j].c_str());
                                        plogf(" and share nonzero stoic: %-9.1f",
                                              m_stoichCoeffRxnMatrix[i][j]);
                                        plogendl();
                                    }
#endif
                                    goto L_COMPONENT_CALC;
                                }
                            }
                        }
                    }
#endif
                }
            }
        }
    }
#ifdef DEBUG_MODE
    if (m_debug_print_lvl >= 2) {
        plogf("   --- Check for an optimum basis passed");
        plogendl();
    }
#endif
    /*************************************************************************/
    /********************** RE-EVALUATE MAJOR-MINOR VECTOR IF NECESSARY ******/
    /*************************************************************************/
    /*
     *     Skip this section if we haven't done a full calculation.
     *     Go right to the check equilibrium section
     */
    if (iti == 0) {
#ifdef DEBUG_MODE
        if (m_debug_print_lvl >= 2) {
            plogf("   --- Reevaluate major-minor status of noncomponents:\n");
        }
#endif
        m_numRxnMinorZeroed = 0;
        for (irxn = 0; irxn < m_numRxnRdc; irxn++) {
            kspec = m_indexRxnToSpecies[irxn];

            int speciesType = vcs_species_type(kspec);
            if (speciesType < VCS_SPECIES_MINOR) {
#ifdef DEBUG_MODE
                if (m_debug_print_lvl >= 2) {
                    if (m_speciesStatus[kspec] >= VCS_SPECIES_MINOR) {
                        plogf("   ---    major/minor species is now zeroed out: %s\n",
                              m_speciesName[kspec].c_str());
                    }
                }
#endif
                ++m_numRxnMinorZeroed;
            } else if (speciesType == VCS_SPECIES_MINOR) {
#ifdef DEBUG_MODE
                if (m_debug_print_lvl >= 2) {
                    if (m_speciesStatus[kspec] != VCS_SPECIES_MINOR) {
                        if (m_speciesStatus[kspec] == VCS_SPECIES_MAJOR) {
                            plogf("   ---   Noncomponent turned from major to minor: ");
                        } else if (kspec < m_numComponents) {
                            plogf("   ---   Component turned into a minor species: ");
                        } else {
                            plogf("   ---   Zeroed Species turned into a "
                                  "minor species: ");
                        }
                        plogf("%s\n", m_speciesName[kspec].c_str());
                    }
                }
#endif
                ++m_numRxnMinorZeroed;
            } else if (speciesType == VCS_SPECIES_MAJOR) {
                if (m_speciesStatus[kspec] != VCS_SPECIES_MAJOR) {
#ifdef DEBUG_MODE
                    if (m_debug_print_lvl >= 2) {
                        if (m_speciesStatus[kspec] ==  VCS_SPECIES_MINOR) {
                            plogf("   ---   Noncomponent turned from minor to major: ");
                        } else if (kspec < m_numComponents) {
                            plogf("   ---   Component turned into a major: ");
                        } else {
                            plogf("   ---   Noncomponent turned from zeroed to major: ");
                        }
                        plogf("%s\n", m_speciesName[kspec].c_str());
                    }
#endif
                    m_speciesStatus[kspec] = VCS_SPECIES_MAJOR;
                    /*
                     *   For this special case, we must reevaluate thermo functions
                     */
                    if (iti != 0) {
                        vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
                        vcs_dfe(VCS_STATECALC_OLD, 0, kspec, kspec+1);
                        vcs_deltag(0, false, VCS_STATECALC_OLD);
                    }
                }
            }
            m_speciesStatus[kspec] = speciesType;
        }
        /*
         *         This logical variable indicates whether all current
         *         non-component species are minor or nonexistent
         */
        allMinorZeroedSpecies = (m_numRxnMinorZeroed == m_numRxnRdc);
    }
    /*************************************************************************/
    /***************** EQUILIBRIUM CHECK FOR MAJOR SPECIES *******************/
    /*************************************************************************/
L_EQUILIB_CHECK:
    ;
    if (! allMinorZeroedSpecies) {
#ifdef DEBUG_MODE
        if (m_debug_print_lvl >= 2) {
            plogf("   --- Equilibrium check for major species: ");
        }
#endif
        for (irxn = 0; irxn < m_numRxnRdc; ++irxn) {
            kspec = irxn + m_numComponents;
            if (m_speciesStatus[kspec] == VCS_SPECIES_MAJOR && (fabs(m_deltaGRxn_new[irxn]) > m_tolmaj)) {
                if (m_VCount->Its >= maxit) {
                    solveFail = -1;
                    /*
                     *         Clean up and exit code even though we haven't
                     *         converged. -> we have run out of iterations!
                     */
                    goto L_RETURN_BLOCK;
                } else {
#ifdef DEBUG_MODE
                    if (m_debug_print_lvl >= 2) {
                        plogf("%s failed\n", m_speciesName[m_indexRxnToSpecies[irxn]].c_str());
                    }
#endif
                    // Convergence amongst major species has not been achieved
                    /*
                     *   Go back and do another iteration with variable ITI
                     */
                    if (forceComponentCalc) {
                        goto L_COMPONENT_CALC;
                    }
                    goto L_MAINLOOP_MM4_SPECIES;
                }
            }
        }
#ifdef DEBUG_MODE
        if (m_debug_print_lvl >= 2) {
            plogf(" MAJOR SPECIES CONVERGENCE achieved");
            plogendl();
        }
#endif
    }
#ifdef DEBUG_MODE
    else {
        if (m_debug_print_lvl >= 2) {
            plogf(" MAJOR SPECIES CONVERGENCE achieved "
                  "(because there are no major species)");
            plogendl();
        }
    }
#endif
    // Convergence amongst major species has been achieved

    /*************************************************************************/
    /*************** EQUILIBRIUM CHECK FOR MINOR SPECIES *********************/
    /*************************************************************************/
    if (m_numRxnMinorZeroed != 0) {
        /*
         *       Calculate the chemical potential and reaction DeltaG
         *       for minor species, if needed.
         */
        if (iti != 0) {
            vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
            vcs_dfe(VCS_STATECALC_OLD, 1, 0, m_numSpeciesRdc);
            vcs_deltag(1, false, VCS_STATECALC_OLD);
            uptodate_minors = true;
        }
#ifdef DEBUG_MODE
        if (m_debug_print_lvl >= 2) {
            plogf("   --- Equilibrium check for minor species: ");
        }
#endif
        for (irxn = 0; irxn < m_numRxnRdc; ++irxn) {
            kspec = irxn + m_numComponents;
            if (m_speciesStatus[kspec] == VCS_SPECIES_MINOR && (fabs(m_deltaGRxn_new[irxn]) > m_tolmin)) {
                if (m_VCount->Its >= maxit) {
                    solveFail = -1;
                    /*
                     *       Clean up and exit code. -> Even though we have not
                     *        converged, we have run out of iterations !
                     */
                    goto L_RETURN_BLOCK;
                }
#ifdef DEBUG_MODE
                if (m_debug_print_lvl >= 2) {
                    plogf("%s failed\n", m_speciesName[m_indexRxnToSpecies[irxn]].c_str());
                }
#endif
                /*
                 *  Set iti to zero to force a full calculation, and go back
                 *  to the main loop to do another iteration.
                 */
                iti = 0;
                if (forceComponentCalc) {
                    goto L_COMPONENT_CALC;
                }
                goto L_MAINLOOP_ALL_SPECIES;
            }
        }
#ifdef DEBUG_MODE
        if (m_debug_print_lvl >= 2) {
            plogf(" CONVERGENCE achieved\n");
        }
#endif
    }
    /*************************************************************************/
    /*********************** FINAL ELEMENTAL ABUNDANCE CHECK *****************/
    /*************************************************************************/
    /*
     *    Recalculate the element abundance vector again
     */
    vcs_updateVP(VCS_STATECALC_OLD);
    vcs_elab();

    /*        LEC is only true when we are near the end game */
    if (lec) {
        if (!giveUpOnElemAbund) {
#ifdef DEBUG_MODE
            if (m_debug_print_lvl >= 2) {
                plogf("   --- Check the Full Element Abundances: ");
            }
#endif
            /*
             *  Final element abundance check:
             *        If we fail then we need to go back and correct
             *        the element abundances, and then go do a major step
             */
            if (! vcs_elabcheck(1)) {
#ifdef DEBUG_MODE
                if (m_debug_print_lvl >= 2) {
                    if (! vcs_elabcheck(0)) {
                        plogf(" failed\n");
                    } else {
                        plogf(" passed for NC but failed for NE: RANGE ERROR\n");
                    }
                }
#endif
                // delete?
                goto L_ELEM_ABUND_CHECK;
            }
#ifdef DEBUG_MODE
            if (m_debug_print_lvl >= 2) {
                plogf(" passed\n");
            }
#endif
        }
        /*
         *   If we have deleted a species then we need to recheck the
         *   the deleted species, before exiting
         */
        if (m_numSpeciesRdc != m_numSpeciesTot) {
            goto L_RECHECK_DELETED;
        }
        /* - Final checks are passed -> go check out */
        goto L_RETURN_BLOCK;
    }
    lec = true;
    /* *************************************************** */
    /* **** CORRECT ELEMENTAL ABUNDANCES ***************** */
    /* *************************************************** */
L_ELEM_ABUND_CHECK:
    ;
    /*
     *  HKM - Put in an element abundance check. The element abundances
     *        were being corrected even if they were perfectly OK to
     *        start with. This is actually an expensive operation, so
     *        I took it out. Also vcs_dfe() doesn't need to be called if
     *        no changes were made.
     */
    rangeErrorFound = 0;
    if (! vcs_elabcheck(1)) {
        bool ncBefore = vcs_elabcheck(0);
        vcs_elcorr(VCS_DATA_PTR(sm), VCS_DATA_PTR(wx));
        bool ncAfter = vcs_elabcheck(0);
        bool neAfter = vcs_elabcheck(1);
        /*
         *      Go back to evaluate the total moles of gas and liquid.
         */
        vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
        vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc);
        vcs_deltag(0, false, VCS_STATECALC_OLD);
        if (!ncBefore) {
            if (ncAfter) {
                /*
                 * We have breathed new life into the old problem. Now the
                 * element abundances up to NC agree. Go back and
                 * restart the main loop calculation, resetting the
                 * end conditions.
                 */
                lec = false;
                iti = 0;
                goto L_MAINLOOP_ALL_SPECIES;
            } else {
                /*
                 * We are still hosed
                 */
                if (finalElemAbundAttempts >= 3) {
                    giveUpOnElemAbund = true;
                    goto L_EQUILIB_CHECK;
                } else {
                    finalElemAbundAttempts++;
                    lec = false;
                    iti = 0;
                    goto L_MAINLOOP_ALL_SPECIES;
                }
            }
        } else {
            if (ncAfter) {
                if (neAfter) {
                    /*
                     * Recovery of end element abundances
                     * -> go do equilibrium check again and then
                     *    check out.
                     */
                    goto L_EQUILIB_CHECK;
                } else {
                    /*
                     * Probably an unrecoverable range error
                     */
#ifdef DEBUG_MODE
                    if (m_debug_print_lvl >= 2) {
                        plogf(" ---  vcs_solve_tp: RANGE SPACE ERROR ENCOUNTERED\n");
                        plogf(" ---  vcs_solve_tp: - Giving up on NE Element Abundance satisfaction \n");
                        plogf(" ---  vcs_solve_tp: - However, NC Element Abundance criteria is satisfied \n");
                        plogf(" ---  vcs_solve_tp: - Returning the calculated equilibrium condition ");
                        plogendl();
                    }
#endif
                    rangeErrorFound = 1;
                    giveUpOnElemAbund = true;
                    goto L_EQUILIB_CHECK;
                }
            }
        }
    }
    // Calculate delta g's
    vcs_deltag(0, false, VCS_STATECALC_OLD);
    // Go back to equilibrium check as a prep to eventually checking out
    goto L_EQUILIB_CHECK;

    /* *************************************************** */
    /* **** RECHECK DELETED SPECIES ********************** */
    /* *************************************************** */
    /*
     *         We are here for two reasons. One is if we have
     *         achieved convergence, but some species have been eliminated
     *         from the problem because they were in multispecies phases
     *         and their mole fractions drifted less than
     *         VCS_DELETE_SPECIES_CUTOFF .
     *         The other reason why we are here is because all of the
     *         non-component species in the problem have been eliminated
     *         for one reason or another.
     */
L_RECHECK_DELETED:
    ;
    npb = vcs_recheck_deleted();
    /*
     *        If we haven't found any species that needed adding we are done.
     */
    if (npb <= 0) {
        goto L_RETURN_BLOCK_B;
    }
    /*
     *        If we have found something to add, recalculate everything
     *        for minor species and go back to do a full iteration
     */
    vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
    vcs_dfe(VCS_STATECALC_OLD, 1, 0, m_numSpeciesRdc);
    vcs_deltag(0, false, VCS_STATECALC_OLD);
    iti = 0;
    goto L_MAINLOOP_ALL_SPECIES;
    /*************************************************************************/
    /******************** CLEANUP AND RETURN BLOCK ***************************/
    /*************************************************************************/
L_RETURN_BLOCK:
    ;

    npb = vcs_recheck_deleted();
    /*
     *        If we haven't found any species that needed adding we are done.
     */
    if (npb > 0) {
        /*
         *        If we have found something to add, recalculate everything
         *        for minor species and go back to do a full iteration
         */
        vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
        vcs_dfe(VCS_STATECALC_OLD, 1, 0, m_numSpeciesRdc);
        vcs_deltag(0, false, VCS_STATECALC_OLD);
        iti = 0;
        goto L_MAINLOOP_ALL_SPECIES;
    }

L_RETURN_BLOCK_B:
    ;

    /*
     *  Add back deleted species in non-zeroed phases. Estimate their
     *  mole numbers.
     */
    npb = vcs_add_all_deleted();
    if (npb > 0) {
        iti = 0;
#ifdef DEBUG_MODE
        if (m_debug_print_lvl >= 1) {
            plogf("  --- add_all_deleted(): some rxns not converged. RETURNING TO LOOP!");
            plogendl();
        }
#endif
        goto L_MAINLOOP_ALL_SPECIES;
    }

    /*
     * Make sure the volume phase objects hold the same state and
     * information as the vcs object. This also update the Cantera objects
     * with this information.
     */
    vcs_updateVP(VCS_STATECALC_OLD);
    /*
     *    Store the final Delta G values for each non-component species
     *    in the species slot rather than the reaction slot
     */
    // kspec = m_numSpeciesTot;
    // i = m_numRxnTot;
    //for (irxn = 0; irxn < m_numRxnTot; ++irxn) {
    //   --kspec;
    // --i;
    //  m_deltaGRxn_new[kspec] = m_deltaGRxn_new[i];
    //}
    // vcs_dzero(VCS_DATA_PTR(m_deltaGRxn_new), m_numComponents);
    /*
     *       Evaluate the final mole fractions
     *        storing them in wt[]
     */
    vcs_vdzero(m_molNumSpecies_new, m_numSpeciesTot);
    for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
        if (m_SSPhase[kspec]) {
            m_molNumSpecies_new[kspec] = 1.0;
        } else {
            iph = m_phaseID[kspec];
            if (m_tPhaseMoles_old[iph] != 0.0) {
                m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] / m_tPhaseMoles_old[iph];
            } else {
                /*
                 * For MultiSpecies phases that are zeroed out,
                 * return the mole fraction vector from the VolPhase object.
                 * This contains the mole fraction that would be true if
                 * the phase just pops into existence.
                 */
                i = m_speciesLocalPhaseIndex[kspec];
                Vphase = m_VolPhaseList[iph];
                m_molNumSpecies_new[kspec] = Vphase->molefraction(i);
            }
        }
    }
    // Return an error code if a Range Space Error is thought to have occurred.
    if (rangeErrorFound) {
        solveFail = 1;
    }
    /*
     *           Free temporary storage used in this routine
     *           and increment counters
     */
    /*
     *           Calculate counters
     */
    double tsecond = ticktock.secondsWC();
    m_VCount->Time_vcs_TP = tsecond;
    m_VCount->T_Time_vcs_TP += m_VCount->Time_vcs_TP;
    (m_VCount->T_Calls_vcs_TP)++;
    m_VCount->T_Its += m_VCount->Its;
    m_VCount->T_Basis_Opts += m_VCount->Basis_Opts;
    m_VCount->T_Time_basopt += m_VCount->Time_basopt;
    /*
     *          Return a Flag indicating whether convergence occurred
     */
    return solveFail;
}
/*********************************************************************************/

double VCS_SOLVE::vcs_minor_alt_calc(size_t kspec, size_t irxn, bool* do_delete
#ifdef DEBUG_MODE
                                     , char* ANOTE
#endif
                                    ) const
{
    double dx = 0.0, a;
    double w_kspec = m_molNumSpecies_old[kspec];
    double molNum_kspec_new;
    double wTrial, tmp;
    double dg_irxn = m_deltaGRxn_old[irxn];
    doublereal  s;
    size_t iph = m_phaseID[kspec];

    *do_delete = false;
    if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
        if (w_kspec <= 0.0) {
            w_kspec = VCS_DELETE_MINORSPECIES_CUTOFF;
        }
        if (dg_irxn < -200.) {
            dg_irxn = -200.;
        }
#ifdef DEBUG_MODE
        sprintf(ANOTE,"minor species alternative calc");
#endif
        if (dg_irxn >= 23.0) {
            molNum_kspec_new = w_kspec * 1.0e-10;
            if (w_kspec < VCS_DELETE_MINORSPECIES_CUTOFF) {
                goto L_ZERO_SPECIES;
            }
            dx = molNum_kspec_new - w_kspec;
            return dx;
        } else {
            if (fabs(dg_irxn) <= m_tolmin2) {
                molNum_kspec_new = w_kspec;
                return 0.0;
            }
        }

        /*
         * get the diagonal of the activity coefficient jacobian
         */
        s = m_np_dLnActCoeffdMolNum[kspec][kspec] / (m_tPhaseMoles_old[iph]);
        // s *= (m_tPhaseMoles_old[iph]);
        /*
         *   We fit it to a power law approximation of the activity coefficient
         *
         *      gamma = gamma_0 * ( x / x0)**a
         *
         *   where a is forced to be a little bit greater than -1.
         *   We do this so that the resulting expression is always nonnegative
         *
         *   We then solve the resulting calculation:
         *
         *      gamma * x  = gamma_0 * x0 exp (-deltaG/RT);
         *
         *
         */
        a = w_kspec * s;
        if (a < (-1.0 + 1.0E-8)) {
            a = -1.0 + 1.0E-8;
        } else  if (a > 100.0) {
            a = 100.0;
        }
        tmp = -dg_irxn / (1.0 + a);
        if (tmp < -200.) {
            tmp = -200.;
        } else if (tmp > 200.) {
            tmp = 200.;
        }
        wTrial =  w_kspec * exp(tmp);
        // wTrial = w_kspec * exp(-dg_irxn);

        molNum_kspec_new = wTrial;

        if (wTrial > 100. * w_kspec) {
            double molNumMax = 0.0001 * m_tPhaseMoles_old[iph];
            if (molNumMax < 100. * w_kspec) {
                molNumMax =  100. * w_kspec;
            }
            if (wTrial > molNumMax) {
                molNum_kspec_new = molNumMax;
            } else {
                molNum_kspec_new = wTrial;
            }

        } else if (1.0E10 * wTrial < w_kspec) {
            molNum_kspec_new= 1.0E-10 * w_kspec;
        } else {
            molNum_kspec_new = wTrial;
        }


        if ((molNum_kspec_new) < VCS_DELETE_MINORSPECIES_CUTOFF) {
            goto L_ZERO_SPECIES;
        }
        dx = molNum_kspec_new - w_kspec;
        return dx;
        /*
         *
         *  Alternate return based for cases where we need to delete the species
         *  from the current list of active species, because its concentration
         *  has gotten too small.
         */
L_ZERO_SPECIES:
        ;
        *do_delete = true;
        dx = - w_kspec;
        return dx;
    } else {
        /*
         * Voltage calculation
         *   Need to check the sign -> This is good for electrons
         */
        dx = m_deltaGRxn_old[irxn]/ m_Faraday_dim;
#ifdef DEBUG_MODE
        sprintf(ANOTE,"voltage species alternative calc");
#endif
    }
    return dx;
}
/*****************************************************************************/

//  Change the concentration of a species by delta moles.
/*
 *  Make sure to conserve elements and keep track of the total kmoles in all phases.
 *
 *
 *  @param kspec The species index
 *  @param delta_ptr   pointer to the delta for the species. This may change during
 *                     the calculation
 *
 *  @return
 *      1: succeeded without change of dx
 *      0: Had to adjust dx, perhaps to zero, in order to do the delta.
 */
int VCS_SOLVE::delta_species(const size_t kspec, double* const delta_ptr)
{
    size_t irxn = kspec - m_numComponents;
    int retn = 1;
    double tmp;
    double delta = *delta_ptr;
#ifdef DEBUG_MODE
    if (kspec < m_numComponents) {
        plogf("  --- delete_species() ERROR: called for a component %d", kspec);
        plogendl();
        exit(EXIT_FAILURE);
    }
#endif
    if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
        /*
         * Attempt the given dx. If it doesn't work, try to see if a smaller
         * one would work,
         */
        double dx = delta;
        double* sc_irxn = m_stoichCoeffRxnMatrix[irxn];
        for (size_t j = 0; j < m_numComponents; ++j) {
            if (m_molNumSpecies_old[j] > 0.0) {
                tmp = sc_irxn[j] * dx;
                if (-tmp > m_molNumSpecies_old[j]) {
                    retn = 0;
                    dx = std::min(dx, - m_molNumSpecies_old[j] / sc_irxn[j]);
                }
            }
            /*
             * If the component has a zero concentration and is a reactant
             * in the formation reaction, then dx == 0.0, and we just return.
             */
            if (m_molNumSpecies_old[j] <= 0.0) {
                if (sc_irxn[j] < 0.0) {
                    *delta_ptr = 0.0;
                    return 0;
                }
            }
        }
        /*
         * ok, we found a positive dx. implement it.
         */
        *delta_ptr = dx;
        m_molNumSpecies_old[kspec] += dx;
        size_t iph = m_phaseID[kspec];
        m_tPhaseMoles_old[iph] += dx;
        vcs_setFlagsVolPhase(iph, false, VCS_STATECALC_OLD);

        for (size_t j = 0; j < m_numComponents; ++j) {
            tmp = sc_irxn[j] * dx;
            if (tmp != 0.0) {
                iph = m_phaseID[j];
                m_molNumSpecies_old[j] += tmp;
                m_tPhaseMoles_old[iph] += tmp;
                vcs_setFlagsVolPhase(iph, false, VCS_STATECALC_OLD);
                if (m_molNumSpecies_old[j] < 0.0) {
                    m_molNumSpecies_old[j] = 0.0;
                }
            }
        }
    }
    return retn;
}
/*****************************************************************************/

//  Zero out the concentration of a species.
/*
 *     Zero out the concentration of a species. Make sure to conserve
 *  elements and keep track of the total moles in all phases.
 *       w[]
 *       m_tPhaseMoles_old[]
 *
 *  @param kspec  Species index
 *
 *  @return:
 *      1: succeeded
 *      0: failed.
 */
int VCS_SOLVE::vcs_zero_species(const size_t kspec)
{
    int retn = 1;
    /*
     * Calculate a delta that will eliminate the species.
     */
    if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
        double dx = -(m_molNumSpecies_old[kspec]);
        if (dx != 0.0) {
            retn = delta_species(kspec, &dx);
#ifdef DEBUG_MODE
            if (!retn) {
                if (m_debug_print_lvl >= 1) {
                    plogf("vcs_zero_species: Couldn't zero the species %d, "
                          "did delta of %g. orig conc of %g",
                          kspec, dx, m_molNumSpecies_old[kspec] + dx);
                    plogendl();
                }
            }
#endif
        }
    }
    return retn;
}
/**************************************************************************/

// Change a single species from active to inactive status
/*
 *   Rearrange data when species is added or removed. The kspec species is
 *   moved to the back of the active species vector. The back of the species
 *   vector is indicated by the value of  m_numSpeciesRdc, the current
 *   number of active species in the mechanism.
 *
 * @param kspec   Species Index
 * @return
 *     Returns 0 unless.
 *     The return is 1 when the current number of
 *     noncomponent species is equal to zero. A recheck of deleted species
 *     is carried out in the main code.
 */
int VCS_SOLVE::vcs_delete_species(const size_t kspec)
{
    const size_t klast = m_numSpeciesRdc - 1;
    const size_t iph = m_phaseID[kspec];
    vcs_VolPhase* const Vphase = m_VolPhaseList[iph];
    const size_t irxn = kspec - m_numComponents;
    /*
     * Zero the concentration of the species.
     *     -> This zeroes w[kspec] and modifies m_tPhaseMoles_old[]
     */
    const int retn = vcs_zero_species(kspec);
    if (DEBUG_MODE_ENABLED && !retn) {
        plogf("Failed to delete a species!");
        plogendl();
        exit(EXIT_FAILURE);
    }
    /*
     *    Decrement the minor species counter if the current species is
     *    a minor species
     */
    if (m_speciesStatus[kspec] != VCS_SPECIES_MAJOR) {
        --(m_numRxnMinorZeroed);
    }
    m_speciesStatus[kspec] = VCS_SPECIES_DELETED;
    m_deltaGRxn_new[irxn] = 0.0;
    m_deltaGRxn_old[irxn] = 0.0;
    m_feSpecies_new[kspec] = 0.0;
    m_feSpecies_old[kspec] = 0.0;
    m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
    /*
     *    Rearrange the data if the current species isn't the last active
     *    species.
     */
    if (kspec != klast) {
        vcs_switch_pos(true, klast, kspec);
    }
    /*
     *       Adjust the total moles in a phase downwards.
     */
    Vphase->setMolesFromVCSCheck(VCS_STATECALC_OLD,
                                 VCS_DATA_PTR(m_molNumSpecies_old),
                                 VCS_DATA_PTR(m_tPhaseMoles_old));

    /*
     *    Adjust the current number of active species and reactions counters
     */
    --(m_numRxnRdc);
    --(m_numSpeciesRdc);

    /*
     *    Check to see whether we have just annihilated a multispecies phase.
     *    If it is extinct, call the delete_multiphase() function.
     */
    if (! m_SSPhase[klast]) {
        if (Vphase->exists() != VCS_PHASE_EXIST_ALWAYS) {
            bool stillExists = false;
            for (size_t k = 0; k < m_numSpeciesRdc; k++) {
                if (m_speciesUnknownType[k] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
                    if (m_phaseID[k] == iph) {
                        if (m_molNumSpecies_old[k] > 0.0) {
                            stillExists = true;
                            break;
                        }
                    }
                }
            }
            if (!stillExists) {
                vcs_delete_multiphase(iph);
            }
        }
    }
    /*
     *    When the total number of noncomponent species is zero, we
     *    have to signal the calling code
     */
    return (m_numRxnRdc == 0);
}
/***************************************************************************/

/*
 *
 *  reinsert_deleted():
 *
 * irxn = id of the noncomponent species formation reaction for the
 *        species to be added in.
 *
 * We make decisions on the initial mole number, and major-minor status
 * here. We also fix up the total moles in a phase.
 *
 * The algorithm proceeds to implement these decisions in the previous
 * position of the species. Then, vcs_switch_pos is called to move the
 * species into the last active species slot, incrementing the number
 * of active species at the same time.
 *
 * This routine is responsible for the global data manipulation only.
 */
void VCS_SOLVE::vcs_reinsert_deleted(size_t kspec)
{
    size_t k;
    // int irxn = kspec - m_numComponents;
    size_t iph = m_phaseID[kspec];
    double dx;
#ifdef DEBUG_MODE
    if (m_debug_print_lvl >= 2) {
        plogf("   --- Add back a deleted species: %-12s\n", m_speciesName[kspec].c_str());
    }
#endif
    /*
     * Set the species back to minor species status
     *  this adjusts m_molNumSpecies_old[] and m_tPhaseMoles_old[]
     * HKM -> make this a relative mole number!
     */
    dx = m_tPhaseMoles_old[iph] * VCS_RELDELETE_SPECIES_CUTOFF * 10.;
    delta_species(kspec, &dx);
    m_speciesStatus[kspec] = VCS_SPECIES_MINOR;

    if (m_SSPhase[kspec]) {
        m_speciesStatus[kspec] = VCS_SPECIES_MAJOR;
        --(m_numRxnMinorZeroed);
    }

    vcs_VolPhase* Vphase = m_VolPhaseList[iph];
    Vphase->setMolesFromVCSCheck(VCS_STATECALC_OLD,
                                 VCS_DATA_PTR(m_molNumSpecies_old),
                                 VCS_DATA_PTR(m_tPhaseMoles_old));
    /*
     *   We may have popped a multispecies phase back
     *   into existence. If we did, we have to check
     *   the other species in that phase.
     *   Take care of the m_speciesStatus[] flag.
     *   The value of m_speciesStatus[] must change from
     *   VCS_SPECIES_ZEROEDPHASE to VCS_SPECIES_ZEROEDMS
     *   for those other species.
     */
    if (! m_SSPhase[kspec]) {
        if (Vphase->exists() == VCS_PHASE_EXIST_NO) {
            Vphase->setExistence(VCS_PHASE_EXIST_YES);
            for (k = 0; k < m_numSpeciesTot; k++) {
                if (m_phaseID[k] == iph) {
                    if (m_speciesStatus[k] != VCS_SPECIES_DELETED) {
                        m_speciesStatus[k] = VCS_SPECIES_MINOR;
                    }
                }
            }
        }
    } else {
        Vphase->setExistence(VCS_PHASE_EXIST_YES);
    }

    ++(m_numRxnRdc);
    ++(m_numSpeciesRdc);
    ++(m_numRxnMinorZeroed);

    if (kspec != (m_numSpeciesRdc - 1)) {
        /*
         *  Rearrange both the species and the non-component global data
         */
        vcs_switch_pos(true, (m_numSpeciesRdc - 1), kspec);
    }
}
/****************************************************************************/

// This routine handles the bookkeeping involved with the
//  deletion of multiphase phases from the problem.
/*
 *   When they are deleted, all of their species become active
 *   species, even though their mole numbers are set to zero.
 *   This routine does not make the decision to eliminate multiphases.
 *
 *   Note, species in phases with zero mole numbers are still
 *   considered active. Whether the phase pops back into
 *   existence or not is checked as part of the main iteration
 *   loop.
 *
 * @param iph Phase to be deleted
 */
bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph)
{
    size_t kspec, irxn;
    double dx;
    vcs_VolPhase* Vphase = m_VolPhaseList[iph];
    bool successful = true;
    /*
     * set the phase existence flag to dead
     */
    Vphase->setTotalMoles(0.0);
#ifdef DEBUG_MODE
    if (m_debug_print_lvl >= 2) {
        plogf("   --- delete_multiphase %d, %s\n", iph, Vphase->PhaseName.c_str());
    }
#endif

    /*
     * Loop over all of the species in the phase.
     */
    for (kspec = m_numComponents; kspec < m_numSpeciesRdc; ++kspec) {
        if (m_phaseID[kspec] == iph) {
            if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
                /*
                 * calculate an extent of rxn, dx, that zeroes out the species.
                 */
                dx = - (m_molNumSpecies_old[kspec]);
                double dxTent = dx;

                int retn = delta_species(kspec, &dxTent);
                if (retn != 1) {
                    successful = false;
#ifdef DEBUG_MODE
                    if (m_debug_print_lvl >= 2) {
                        plogf("   --- delete_multiphase %d, %s ERROR problems deleting species %s\n",
                              iph, Vphase->PhaseName.c_str(), m_speciesName[kspec].c_str());
                        plogf("   ---     delta  attempted: %g  achieved: %g   "
                              "  Zeroing it manually\n", dx, dxTent);
                    }
#endif
                    m_molNumSpecies_old[kspec]  = 0.0;
                    m_molNumSpecies_new[kspec]  = 0.0;
                    m_deltaMolNumSpecies[kspec] = 0.0;
                    // recover the total phase moles.
                    vcs_tmoles();
                } else {
                    /*
                     * Set the mole number of that species to zero.
                     */
                    m_molNumSpecies_old[kspec]  = 0.0;
                    m_molNumSpecies_new[kspec]  = 0.0;
                    m_deltaMolNumSpecies[kspec] = 0.0;
                }
                /*
                 * Change the status flag of the species to that of an
                 * zeroed phase
                 */
                m_speciesStatus[kspec] = VCS_SPECIES_ZEROEDMS;
            }
        }
    }

    double dj, dxWant, dxPerm = 0.0, dxPerm2 = 0.0;
    for (size_t kcomp = 0; kcomp < m_numComponents; ++kcomp) {
        if (m_phaseID[kcomp] == iph) {
#ifdef DEBUG_MODE
            if (m_debug_print_lvl >= 2) {
                plogf("   --- delete_multiphase   One of the species is a component %d - %s with mole number %g\n",
                      kcomp, m_speciesName[kcomp].c_str(),   m_molNumSpecies_old[kcomp]);
            }
#endif
            if (m_molNumSpecies_old[kcomp] != 0.0) {
                for (kspec = m_numComponents; kspec < m_numSpeciesRdc; ++kspec) {
                    irxn = kspec - m_numComponents;
                    if (m_phaseID[kspec] != iph) {
                        if (m_stoichCoeffRxnMatrix[irxn][kcomp] != 0.0) {
                            dxWant = -m_molNumSpecies_old[kcomp] / m_stoichCoeffRxnMatrix[irxn][kcomp];
                            if (dxWant + m_molNumSpecies_old[kspec]  < 0.0) {
                                dxPerm = -m_molNumSpecies_old[kspec];
                            }
                            for (size_t jcomp = 0;  kcomp < m_numComponents; ++kcomp) {
                                if (jcomp != kcomp) {
                                    if (m_phaseID[jcomp] == iph) {
                                        dxPerm = 0.0;
                                    } else {
                                        dj = dxWant * m_stoichCoeffRxnMatrix[irxn][jcomp];
                                        if (dj + m_molNumSpecies_old[kcomp]  < 0.0) {
                                            dxPerm2 = -m_molNumSpecies_old[kcomp] / m_stoichCoeffRxnMatrix[irxn][jcomp];
                                        }
                                        if (fabs(dxPerm2) < fabs(dxPerm)) {
                                            dxPerm = dxPerm2;
                                        }
                                    }
                                }
                            }
                        }
                        if (dxPerm != 0.0) {
                            delta_species(kspec, &dxPerm);
                        }
                    }
                }

            }
            if (m_molNumSpecies_old[kcomp] != 0.0) {
#ifdef DEBUG_MODE
                if (m_debug_print_lvl >= 2) {
                    plogf("   --- delete_multiphase   One of the species is a component %d - %s still with mole number %g\n",
                          kcomp, m_speciesName[kcomp].c_str(),   m_molNumSpecies_old[kcomp]);
                    plogf("   ---                     zeroing it \n");
                }
#endif
                m_molNumSpecies_old[kcomp] = 0.0;
            }
            m_speciesStatus[kcomp] = VCS_SPECIES_ZEROEDMS;
        }
    }


    /*
     *   Loop over all of the inactive species in the phase:
     *   Right now we reinstate all species in a deleted multiphase.
     *   We may only want to reinstate the "major ones" in the future.
     *   Note, species in phases with zero mole numbers are still
     *   considered active. Whether the phase pops back into
     *   existence or not is checked as part of the main iteration
     *   loop.
     */
    for (kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) {
        if (m_phaseID[kspec] == iph) {
            irxn = kspec - m_numComponents;
            m_molNumSpecies_old[kspec]  = 0.0;
            m_molNumSpecies_new[kspec] = 0.0;
            m_deltaMolNumSpecies[kspec] = 0.0;
            m_speciesStatus[kspec] = VCS_SPECIES_ZEROEDMS;

            ++(m_numRxnRdc);
            ++(m_numSpeciesRdc);
#ifdef DEBUG_MODE
            if (m_debug_print_lvl >= 2) {
                plogf("   ---    Make %s", m_speciesName[kspec].c_str());
                plogf(" an active but zeroed species because its phase "
                      "was zeroed\n");
            }
#endif
            if (kspec != (m_numSpeciesRdc - 1)) {
                /*
                 *  Rearrange both the species and the non-component global data
                 */
                vcs_switch_pos(true, (m_numSpeciesRdc - 1), kspec);
            }
        }
    }

    /*
     * Zero out the total moles counters for the phase
     */
    m_tPhaseMoles_old[iph] = 0.0;
    m_tPhaseMoles_new[iph] = 0.0;
    m_deltaPhaseMoles[iph] = 0.0;
    /*
     * Upload the state to the VP object
     */
    Vphase->setTotalMoles(0.0);

    return successful;
}
/**********************************************************************************/

//    Recheck deleted species in multispecies phases.
/*
 *   We are checking the equation:
 *
 *         sum_u = sum_j_comp [ sigma_i_j * u_j ]
 *               = u_i_O + log((AC_i * W_i)/m_tPhaseMoles_old)
 *
 *   by first evaluating:
 *
 *          DG_i_O = u_i_O - sum_u.
 *
 *   Then, if TL is zero, the phase pops into existence if DG_i_O < 0.
 *   Also, if the phase exists, then we check to see if the species
 *   can have a mole number larger than  VCS_DELETE_SPECIES_CUTOFF
 *   (default value = 1.0E-32).
 *
 */
int VCS_SOLVE::vcs_recheck_deleted()
{
    int npb;
    size_t iph, kspec, irxn;
    double* xtcutoff = VCS_DATA_PTR(m_TmpPhase);
#ifdef DEBUG_MODE
    if (m_debug_print_lvl >= 2) {
        plogf("   --- Start rechecking deleted species in multispec phases\n");
    }
#endif
    if (m_numSpeciesRdc == m_numSpeciesTot) {
        return 0;
    }
    /*
     * Use the standard chemical potentials for the chemical potentials
     * of deleted species. Then, calculate Delta G for
     * for formation reactions.
     * Note: fe[] here includes everything except for the ln(x[i]) term
     */
    for (kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) {
        iph = m_phaseID[kspec];
        m_feSpecies_new[kspec] = (m_SSfeSpecies[kspec] + log(m_actCoeffSpecies_old[kspec])
                                  - m_lnMnaughtSpecies[kspec]
                                  + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iph]);
    }

    /*
     *      Recalculate the DeltaG's of the formation reactions for the
     *      deleted species in the mechanism
     */
    vcs_deltag(0, true, VCS_STATECALC_NEW);

    for (iph = 0; iph < m_numPhases; iph++) {
        if (m_tPhaseMoles_old[iph] > 0.0) {
            xtcutoff[iph] = log(1.0 / VCS_RELDELETE_SPECIES_CUTOFF);
        } else {
            xtcutoff[iph] = 0.0;
        }
    }

    /*
     *
     *   We are checking the equation:
     *
     *         sum_u = sum_j_comp [ sigma_i_j * u_j ]
     *               = u_i_O + log((AC_i * W_i)/m_tPhaseMoles_old)
     *
     *   by first evaluating:
     *
     *          DG_i_O = u_i_O - sum_u.
     *
     *   Then, if TL is zero, the phase pops into existence if DG_i_O < 0.
     *   Also, if the phase exists, then we check to see if the species
     *   can have a mole number larger than  VCS_DELETE_SPECIES_CUTOFF
     *   (default value = 1.0E-32).
     *
     *   HKM:
     *       This seems to be an inconsistency in the algorithm here that needs
     *   correcting. The requirement above may bypass some multiphases which
     *   should exist. The real requirement for the phase to exist is:
     *
     *           sum_i_in_phase [ exp(-DG_i_O) ] >= 1.0
     *
     *   Thus, we need to amend the code. Also nonideal solutions will tend to
     *   complicate matters severely also.
     */
    npb = 0;
    for (irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) {
        kspec = m_indexRxnToSpecies[irxn];
        iph = m_phaseID[kspec];
        if (m_tPhaseMoles_old[iph] == 0.0) {
            if (m_deltaGRxn_new[irxn] < 0.0) {
                vcs_reinsert_deleted(kspec);
                npb++;
            } else {
                m_molNumSpecies_old[kspec] = 0.0;
            }
        } else if (m_tPhaseMoles_old[iph] > 0.0) {
            if (m_deltaGRxn_new[irxn] < xtcutoff[iph]) {
                vcs_reinsert_deleted(kspec);
                npb++;
            }
        }
    }
    return npb;
}
/***********************************************************************************/

bool VCS_SOLVE::recheck_deleted_phase(const int iphase)
{

    // Check first to see if the phase is in fact deleted
    const vcs_VolPhase* Vphase = m_VolPhaseList[iphase];
    if (Vphase->exists() != VCS_PHASE_EXIST_NO) {
        return false;
    }
    if (Vphase->exists() == VCS_PHASE_EXIST_ZEROEDPHASE) {
        return false;
    }
    size_t irxn, kspec;
    if (Vphase->m_singleSpecies) {
        kspec = Vphase->spGlobalIndexVCS(0);
        irxn = kspec + m_numComponents;
        if (m_deltaGRxn_old[irxn] < 0.0) {
            return true;
        }
        return false;
    }

    double phaseDG = 1.0;
    for (size_t kk = 0; kk <  Vphase->nSpecies(); kk++) {
        kspec = Vphase->spGlobalIndexVCS(kk);
        irxn = kspec + m_numComponents;
        if (m_deltaGRxn_old[irxn] >  50.0) {
            m_deltaGRxn_old[irxn] =  50.0;
        }
        if (m_deltaGRxn_old[irxn] < -50.0) {
            m_deltaGRxn_old[irxn] = -50.0;
        }
        phaseDG -= exp(-m_deltaGRxn_old[irxn]);
    }

    if (phaseDG < 0.0) {
        return true;
    }
    return false;
}
//====================================================================================================================
//  Provide an estimate for the deleted species in phases that are not zeroed out
/*
 *  Try to add back in all deleted species. An estimate of the kmol numbers
 *  are obtained and the species is added back into the equation system,
 *  into the old state vector.
 */
size_t VCS_SOLVE::vcs_add_all_deleted()
{
    size_t iph, kspec, retn;
    if (m_numSpeciesRdc == m_numSpeciesTot) {
        return 0;
    }
    /*
     * Use the standard chemical potentials for the chemical potentials of deleted species. Then, calculate Delta G for
     * for formation reactions.
     *     We are relying here on a old saved value of m_actCoeffSpecies_old[kspec]
     *  being sufficiently good. Note, we will recalculate everything at the end of the routine.
     */
    vcs_dcopy(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_molNumSpecies_old), m_numSpeciesTot);

    for (int cits = 0; cits < 3; cits++) {
        for (kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) {
            iph = m_phaseID[kspec];
            vcs_VolPhase* Vphase = m_VolPhaseList[iph];
            if (m_molNumSpecies_new[kspec] == 0.0) {
                m_molNumSpecies_new[kspec] = VCS_DELETE_MINORSPECIES_CUTOFF * 1.0E-10;
            }
            if (!Vphase->m_singleSpecies) {
                Vphase->sendToVCS_ActCoeff(VCS_STATECALC_NEW, VCS_DATA_PTR(m_actCoeffSpecies_new));
            }
            m_feSpecies_new[kspec] = (m_SSfeSpecies[kspec] + log(m_actCoeffSpecies_new[kspec]) - m_lnMnaughtSpecies[kspec]
                                      + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iph]);
        }
        /*
         * Recalculate the DeltaG's of the formation reactions for the deleted species in the mechanism
         */
        vcs_deltag(0, true, VCS_STATECALC_NEW);
        for (size_t irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) {
            kspec = m_indexRxnToSpecies[irxn];
            iph = m_phaseID[kspec];
            if (m_tPhaseMoles_old[iph] > 0.0) {
                double maxDG = std::min(m_deltaGRxn_new[irxn], 690.0);
                double dx = m_tPhaseMoles_old[iph] * exp(- maxDG);
                m_molNumSpecies_new[kspec] = dx;
                if (m_molNumSpecies_new[kspec] > 2 *VCS_DELETE_MINORSPECIES_CUTOFF) {
                    m_molNumSpecies_new[kspec] = 2 * VCS_DELETE_MINORSPECIES_CUTOFF;
                }
            }
        }
    }
    for (size_t irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) {
        kspec = m_indexRxnToSpecies[irxn];
        iph = m_phaseID[kspec];
        if (m_tPhaseMoles_old[iph] > 0.0) {
            double dx =  m_molNumSpecies_new[kspec];
            retn = delta_species(kspec, &dx);
            if (retn == 0) {
#ifdef DEBUG_MODE
                if (m_debug_print_lvl) {
                    plogf("  --- add_deleted(): delta_species() failed for species %s (%d) with mol number %g\n",
                          m_speciesName[kspec].c_str(), kspec, dx);
                }
#endif
                if (dx > 1.0E-50) {
                    dx = 1.0E-50;
                    retn = delta_species(kspec, &dx);
#ifdef DEBUG_MODE
                    if (retn == 0) {
                        if (m_debug_print_lvl) {
                            plogf("  --- add_deleted(): delta_species() failed for species %s (%d) with mol number %g\n",
                                  m_speciesName[kspec].c_str(), kspec, dx);
                        }
                    }
#endif
                }
            }
#ifdef DEBUG_MODE
            if (m_debug_print_lvl >= 2) {
                if (retn != 0) {
                    plogf("  --- add_deleted():  species %s added back in with mol number %g",
                          m_speciesName[kspec].c_str(), dx);
                    plogendl();
                } else {
                    plogf("  --- add_deleted():  species %s failed to be added  back in");
                    plogendl();
                }
            }
#endif
        }
    }
    vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
    vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesTot);
    vcs_deltag(0, true, VCS_STATECALC_OLD);

    retn = 0;
    for (size_t irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) {
        kspec = m_indexRxnToSpecies[irxn];
        iph = m_phaseID[kspec];
        if (m_tPhaseMoles_old[iph] > 0.0) {
            if (fabs(m_deltaGRxn_old[irxn]) > m_tolmin) {
                if (((m_molNumSpecies_old[kspec] * exp(-m_deltaGRxn_old[irxn])) >
                        VCS_DELETE_MINORSPECIES_CUTOFF) ||
                        (m_molNumSpecies_old[kspec] > VCS_DELETE_MINORSPECIES_CUTOFF)) {
                    retn++;
#ifdef DEBUG_MODE
                    if (m_debug_print_lvl >= 2) {
                        plogf("  --- add_deleted():  species %s with mol number %g not converged: DG = %g",
                              m_speciesName[kspec].c_str(), m_molNumSpecies_old[kspec],
                              m_deltaGRxn_old[irxn]);
                        plogendl();
                    }
#endif
                }
            }
        }
    }
    return retn;
}
/***********************************************************************************/

//  This routine optimizes the minimization of the total gibbs free
//  energy by making sure the slope of the Gibbs free energy stays negative
/*
 *  The slope of the following functional is equivalent to the slope of the total
 *  Gibbs free energy of the system:
 *
 *                 d_Gibbs/ds = sum_k( m_deltaGRxn * m_deltaMolNumSpecies[k] )
 *
 *  along the current direction m_deltaMolNumSpecies[], by choosing a value, al: (0<al<1)
 *  such that the a parabola approximation to Gibbs(al) fit to the
 *  end points al = 0 and al = 1 is minimized.
 *      s1 = slope of Gibbs function at al = 0, which is the previous
 *           solution = d(Gibbs)/d(al).
 *      s2 = slope of Gibbs function at al = 1, which is the current
 *           solution = d(Gibbs)/d(al).
 *  Only if there has been an inflection point (i.e., s1 < 0 and s2 > 0),
 *  does this code section kick in. It finds the point on the parabola
 *  where the slope is equal to zero.
 *
 */
bool VCS_SOLVE::vcs_globStepDamp()
{
    double s1, s2, al;
    size_t irxn, kspec, iph;
    double* dptr = VCS_DATA_PTR(m_deltaGRxn_new);

    /* *************************************************** */
    /* **** CALCULATE SLOPE AT END OF THE STEP  ********** */
    /* *************************************************** */
    s2 = 0.0;
    for (irxn = 0; irxn < m_numRxnRdc; ++irxn) {
        kspec = irxn + m_numComponents;
        if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
            s2 += dptr[irxn] * m_deltaMolNumSpecies[kspec];
        }
    }


    /* *************************************************** */
    /* **** CALCULATE ORIGINAL SLOPE ********************* */
    /* ************************************************** */
    s1 = 0.0;
    dptr = VCS_DATA_PTR(m_deltaGRxn_old);
    for (irxn = 0; irxn < m_numRxnRdc; ++irxn) {
        kspec = irxn + m_numComponents;
        if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
            s1 += dptr[irxn] * m_deltaMolNumSpecies[kspec];
        }
    }

#ifdef DEBUG_MODE
    if (m_debug_print_lvl >= 2) {
        plogf("   --- subroutine FORCE: Beginning Slope = %g\n", s1);
        plogf("   --- subroutine FORCE: End Slope       = %g\n", s2);
    }
#endif

    if (s1 > 0.0) {
#ifdef DEBUG_MODE
        if (m_debug_print_lvl >= 2) {
            plogf("   --- subroutine FORCE produced no adjustments,");
            if (s1 < 1.0E-40) {
                plogf(" s1 positive but really small");
            } else {
                plogf(" failed s1 test");
            }
            plogendl();
        }
#endif
        return false;
    }

    if (s2 <= 0.0) {
#ifdef DEBUG_MODE
        if (m_debug_print_lvl >= 2) {
            plogf("   --- subroutine FORCE produced no adjustments, s2 < 0");
            plogendl();
        }
#endif
        return false;
    }

    /* *************************************************** */
    /* **** FIT PCJ2822ARABOLA ********************************* */
    /* *************************************************** */
    al = 1.0;
    if (fabs(s1 -s2) > 1.0E-200) {
        al = s1 / (s1 - s2);
    }
    if (al >= 0.95 || al < 0.0) {
#ifdef DEBUG_MODE
        if (m_debug_print_lvl >= 2) {
            plogf("   --- subroutine FORCE produced no adjustments (al = %g)\n", al);
        }
#endif
        return false;
    }
#ifdef DEBUG_MODE
    if (m_debug_print_lvl >= 2) {
        plogf("   --- subroutine FORCE produced a damping factor = %g\n", al);
    }
#endif

    /* *************************************************** */
    /* **** ADJUST MOLE NUMBERS, CHEM. POT *************** */
    /* *************************************************** */
#ifdef DEBUG_MODE
    if (m_debug_print_lvl >= 2) {
        vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_tmp), VCS_DATA_PTR(m_deltaGRxn_new),
                  m_numRxnRdc);
    }
#endif

    dptr = VCS_DATA_PTR(m_molNumSpecies_new);
    for (kspec = 0; kspec < m_numSpeciesRdc; ++kspec) {
        m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] +
                                     al * m_deltaMolNumSpecies[kspec];
    }
    for (iph = 0; iph < m_numPhases; iph++) {
        m_tPhaseMoles_new[iph] = m_tPhaseMoles_old[iph] + al * m_deltaPhaseMoles[iph];
    }
    vcs_updateVP(VCS_STATECALC_NEW);

#ifdef DEBUG_MODE
    if (m_debug_print_lvl >= 2) {
        plogf("   --- subroutine FORCE adjusted the mole "
              "numbers, AL = %10.3f\n", al);
    }
#endif
    /*
     *           Because we changed the mole numbers, we need to
     *           calculate the chemical potentials again. If a major-
     *           only step is being carried out, then we don't need to
     *           update the minor noncomponents.
     */
    vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW);
    vcs_dfe(VCS_STATECALC_NEW, 0, 0, m_numSpeciesRdc);

    /*
     *           Evaluate DeltaG for all components if ITI=0, and for
     *           major components only if ITI NE 0
     */
    vcs_deltag(0, false, VCS_STATECALC_NEW);

    dptr = VCS_DATA_PTR(m_deltaGRxn_new);
    s2 = 0.0;
    for (irxn = 0; irxn < m_numRxnRdc; ++irxn) {
        kspec = irxn + m_numComponents;
        if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
            s2 += dptr[irxn] * m_deltaMolNumSpecies[kspec];
        }
    }


#ifdef DEBUG_MODE
    if (m_debug_print_lvl >= 2) {
        plogf("   --- subroutine FORCE: Adj End Slope   = %g", s2);
        plogendl();
    }
#endif
    return true;
}
/****************************************************************************************/

int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], double sm[],
                          double ss[], double test, bool* const usedZeroedSpecies)
{
    size_t  j, k, l, i, jl, ml, jr, irxn, kspec;
    bool lindep;
    size_t ncTrial;
    size_t juse = npos;
    size_t jlose = npos;
    double* dptr, *scrxn_ptr;
    Cantera::clockWC tickTock;
#ifdef DEBUG_MODE
    if (m_debug_print_lvl >= 2) {
        plogf("   ");
        for (i=0; i<77; i++) {
            plogf("-");
        }
        plogf("\n");
        plogf("   --- Subroutine BASOPT called to ");
        if (doJustComponents) {
            plogf("calculate the number of components\n");
        } else {
            plogf("reevaluate the components\n");
        }
        if (m_debug_print_lvl >= 2) {
            plogf("\n");
            plogf("   ---       Formula Matrix used in BASOPT calculation\n");
            plogf("   ---      Active |   ");
            for (j = 0; j < m_numElemConstraints; j++) {
                plogf("  %1d  ", m_elementActive[j]);
            }
            plogf("\n");
            plogf("   ---     Species | ");
            for (j = 0; j < m_numElemConstraints; j++) {
                plogf(" ");
                vcs_print_stringTrunc(m_elementName[j].c_str(), 8, 1);
            }
            plogf("\n");
            for (k = 0; k < m_numSpeciesTot; k++) {
                plogf("   --- ");
                vcs_print_stringTrunc(m_speciesName[k].c_str(), 11, 1);
                plogf(" | ");
                for (j = 0; j < m_numElemConstraints; j++) {
                    plogf(" %8.2g", m_formulaMatrix[j][k]);
                }
                plogf("\n");
            }
            plogendl();
        }
    }
#endif

    /*
     *  Calculate the maximum value of the number of components possible
     *     It's equal to the minimum of the number of elements and the
     *     number of total species.
     */
    ncTrial = std::min(m_numElemConstraints, m_numSpeciesTot);
    m_numComponents = ncTrial;
    *usedZeroedSpecies = false;

    /*
     *     Use a temporary work array for the mole numbers, aw[]
     */
    vcs_dcopy(aw, VCS_DATA_PTR(m_molNumSpecies_old), m_numSpeciesTot);
    /*
     * Take out the Voltage unknowns from consideration
     */
    for (k = 0; k < m_numSpeciesTot; k++) {
        if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
            aw[k] = test;
        }
    }

    jr = npos;
    /*
     *   Top of a loop of some sort based on the index JR. JR is the
     *   current number of component species found.
     */
    do {
        ++jr;
        /* - Top of another loop point based on finding a linearly */
        /* - independent species */
        do {
            /*
             *    Search the remaining part of the mole fraction vector, AW,
             *    for the largest remaining species. Return its identity in K.
             *    The first search criteria is always the largest positive
             *    magnitude of the mole number.
             */
            k = vcs_basisOptMax(aw, jr, m_numSpeciesTot);

            /*
             * The fun really starts when you have run out of species that have a significant
             * concentration. It becomes extremely important to make a good choice of which
             * species you want to pick to fill out the basis. Basically, you don't want to
             * use species with elements abundances which aren't pegged to zero. This means
             * that those modes will never be allowed to grow. You want to have the
             * best chance that the component will grow positively.
             *
             * Suppose  you start with CH4, N2, as the only species with nonzero compositions.
             * You have the following abundances:
             *
             *  Abundances:
             * ----------------
             *     C 2.0
             *     N 2.0
             *     H 4.0
             *     O 0.0
             *
             * For example, Make the following choice:
             *
             *   CH4   N2  O choose ->  OH
             * or
             *   CH4   N2  O choose ->  H2
             *
             *  OH and H2 both fill out the basis. They will pass the algorithm. However,
             *  choosing OH as the next species will create a situation where H2 can not
             *  grow in concentration. This happened in practice, btw. The reason is that
             *  the formation reaction for H2 will cause one of the component species
             *  to go negative.
             *
             *  The basic idea here is to pick a simple species whose mole number
             *  can grow according to the element compositions. Candidates are still
             *  filtered according to their linear independence.
             *
             *   Note, if there is electronic charge and the electron species,
             *   you should probably pick the electron as a component, if it
             *   linearly independent. The algorithm below will do this automagically.
             *
             */
            if ((aw[k] != test) && aw[k] < VCS_DELETE_MINORSPECIES_CUTOFF) {
                *usedZeroedSpecies = true;

                double maxConcPossKspec = 0.0;
                double maxConcPoss = 0.0;
                size_t kfound = npos;
                int minNonZeroes = 100000;
                int nonZeroesKspec = 0;
                for (kspec = ncTrial; kspec < m_numSpeciesTot; kspec++) {
                    if (aw[kspec] >= 0.0) {
                        if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
                            maxConcPossKspec = 1.0E10;
                            nonZeroesKspec = 0;
                            for (size_t j = 0; j < m_numElemConstraints; ++j) {
                                if (m_elementActive[j]) {
                                    if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
                                        double nu = m_formulaMatrix[j][kspec];
                                        if (nu != 0.0) {
                                            nonZeroesKspec++;
                                            maxConcPossKspec = std::min(m_elemAbundancesGoal[j] / nu, maxConcPossKspec);
                                        }
                                    }
                                }
                            }
                            if ((maxConcPossKspec >= maxConcPoss) || (maxConcPossKspec > 1.0E-5)) {
                                if (nonZeroesKspec <= minNonZeroes) {
                                    if (kfound == npos || nonZeroesKspec < minNonZeroes) {
                                        kfound = kspec;
                                    } else {
                                        // ok we are sitting pretty equal here decide on the raw ss Gibbs energy
                                        if (m_SSfeSpecies[kspec] <= m_SSfeSpecies[kfound]) {
                                            kfound = kspec;
                                        }
                                    }
                                }
                                if (nonZeroesKspec < minNonZeroes) {
                                    minNonZeroes = nonZeroesKspec;
                                }
                                if (maxConcPossKspec > maxConcPoss) {
                                    maxConcPoss = maxConcPossKspec;
                                }
                            }
                        }
                    }
                }
                if (kfound == npos) {
                    double gmin = 0.0;
                    kfound = k;
                    for (kspec = ncTrial; kspec < m_numSpeciesTot; kspec++) {
                        if (aw[kspec] >= 0.0) {
                            irxn = kspec - ncTrial;
                            if (m_deltaGRxn_new[irxn] < gmin) {
                                gmin = m_deltaGRxn_new[irxn];
                                kfound = kspec;
                            }
                        }
                    }
                }
                k = kfound;
            }


            if (aw[k] == test) {
                m_numComponents = jr;
                ncTrial = m_numComponents;
                size_t numPreDeleted = m_numRxnTot - m_numRxnRdc;
                if (numPreDeleted != (m_numSpeciesTot - m_numSpeciesRdc)) {
                    plogf("vcs_basopt:: we shouldn't be here\n");
                    exit(EXIT_FAILURE);
                }
                m_numRxnTot = m_numSpeciesTot - ncTrial;
                m_numRxnRdc = m_numRxnTot - numPreDeleted;
                m_numSpeciesRdc = m_numSpeciesTot - numPreDeleted;
                for (i = 0; i < m_numSpeciesTot; ++i) {
                    m_indexRxnToSpecies[i] = ncTrial + i;
                }
#ifdef DEBUG_MODE
                if (m_debug_print_lvl >= 2) {
                    plogf("   ---   Total number of components found = %3d (ne = %d)\n ",
                          ncTrial, m_numElemConstraints);
                }
#endif
                goto L_END_LOOP;
            }
            /*
             *  Assign a small negative number to the component that we have
             *  just found, in order to take it out of further consideration.
             */
            aw[k] = test;
            /* *********************************************************** */
            /* **** CHECK LINEAR INDEPENDENCE WITH PREVIOUS SPECIES ****** */
            /* *********************************************************** */
            /*
             *          Modified Gram-Schmidt Method, p. 202 Dalquist
             *          QR factorization of a matrix without row pivoting.
             */
            jl = jr;
            for (j = 0; j < m_numElemConstraints; ++j) {
                sm[j + jr*m_numElemConstraints] = m_formulaMatrix[j][k];
            }
            if (jl > 0) {
                /*
                 *         Compute the coefficients of JA column of the
                 *         the upper triangular R matrix, SS(J) = R_J_JR
                 *         (this is slightly different than Dalquist)
                 *         R_JA_JA = 1
                 */
                for (j = 0; j < jl; ++j) {
                    ss[j] = 0.0;
                    for (i = 0; i < m_numElemConstraints; ++i) {
                        ss[j] += sm[i + jr*m_numElemConstraints] * sm[i + j*m_numElemConstraints];
                    }
                    ss[j] /= sa[j];
                }
                /*
                 *     Now make the new column, (*,JR), orthogonal to the
                 *     previous columns
                 */
                for (j = 0; j < jl; ++j) {
                    for (l = 0; l < m_numElemConstraints; ++l) {
                        sm[l + jr*m_numElemConstraints] -= ss[j] * sm[l + j*m_numElemConstraints];
                    }
                }
            }
            /*
             *        Find the new length of the new column in Q.
             *        It will be used in the denominator in future row calcs.
             */
            sa[jr] = 0.0;
            for (ml = 0; ml < m_numElemConstraints; ++ml) {
                sa[jr] += SQUARE(sm[ml + jr*m_numElemConstraints]);
            }
            /* **************************************************** */
            /* **** IF NORM OF NEW ROW  .LT. 1E-3 REJECT ********** */
            /* **************************************************** */
            lindep = (sa[jr] < 1.0e-6);
        } while (lindep);
        /* ****************************************** */
        /* **** REARRANGE THE DATA ****************** */
        /* ****************************************** */
        if (jr != k) {
#ifdef DEBUG_MODE
            if (m_debug_print_lvl >= 2) {
                plogf("   ---   %-12.12s", (m_speciesName[k]).c_str());
                plogf("(%9.2g) replaces %-12.12s", m_molNumSpecies_old[k],
                      m_speciesName[jr].c_str());
                plogf("(%9.2g) as component %3d\n", m_molNumSpecies_old[jr], jr);
            }
#endif
            vcs_switch_pos(false, jr, k);
            std::swap(aw[jr], aw[k]);
        }
#ifdef DEBUG_MODE
        else {
            if (m_debug_print_lvl >= 2) {
                plogf("   ---   %-12.12s", m_speciesName[k].c_str());
                plogf("(%9.2g) remains            ", m_molNumSpecies_old[k]);
                plogf("              as component %3d\n", jr);
            }
        }
#endif
        /* - entry point from up above */
L_END_LOOP:
        ;
        /*
         *      If we haven't found enough components, go back
         *      and find some more. (nc -1 is used below, because
         *      jr is counted from 0, via the C convention.
         */
    } while (jr < (ncTrial-1));

    if (doJustComponents) {
        goto L_CLEANUP;
    }
    /* ****************************************************** */
    /* **** EVALUATE THE STOICHIOMETRY ********************** */
    /* ****************************************************** */
    /*
     *  Formulate the matrix problem for the stoichiometric
     *  coefficients. CX + B = 0
     *      C will be an nc x nc matrix made up of the formula
     * vectors for the components.
     *      n rhs's will be solved for. Thus, B is an nc x n
     * matrix.
     *
     * BIG PROBLEM 1/21/99:
     *
     *    This algorithm makes the assumption that the
     * first nc rows of the formula matrix aren't rank deficient.
     * However, this might not be the case. For example, assume
     * that the first element in m_formulaMatrix[] is argon. Assume that
     * no species in the matrix problem actually includes argon.
     * Then, the first row in sm[], below will be identically
     * zero. bleh.
     *    What needs to be done is to perform a rearrangement
     * of the ELEMENTS -> i.e. rearrange, m_formulaMatrix, sp,
     * and m_elemAbundancesGoal, such
     * that the first nc elements form in combination with the
     * nc components create an invertible sm[]. not a small
     * project, but very doable.
     *    An alternative would be to turn the matrix problem
     * below into an ne x nc problem, and do QR elimination instead
     * of Gauss-Jordan elimination.
     *    Note the rearrangement of elements need only be done once
     * in the problem. It's actually very similar to the top of
     * this program with ne being the species and nc being the
     * elements!!
     */
    for (j = 0; j < ncTrial; ++j) {
        for (i = 0; i < ncTrial; ++i) {
            sm[i + j*m_numElemConstraints] = m_formulaMatrix[i][j];
        }
    }
    for (i = 0; i < m_numRxnTot; ++i) {
        k = m_indexRxnToSpecies[i];
        for (j = 0; j < ncTrial; ++j) {
            m_stoichCoeffRxnMatrix[i][j] = m_formulaMatrix[j][k];
        }
    }
    /*
     *     Use Gauss-Jordan block elimination to calculate
     *     the reaction matrix, m_stoichCoeffRxnMatrix[][].
     */

    j = vcsUtil_gaussj(sm, m_numElemConstraints, ncTrial, m_stoichCoeffRxnMatrix[0], m_numRxnTot);
    // j = vcsUtil_mlequ(sm, m_numElemConstraints, ncTrial, m_stoichCoeffRxnMatrix[0], m_numRxnTot);


    if (j == 1) {
        plogf("vcs_solve_TP ERROR: mlequ returned an error condition\n");
        return  VCS_FAILED_CONVERGENCE;
    }

    /*
     * NOW, if we have interfacial voltage unknowns, what we did
     * was just wrong -> hopefully it didn't blow up. Redo the problem.
     * Search for inactive E
     */
    juse  = npos;
    jlose = npos;
    for (j = 0; j < m_numElemConstraints; j++) {
        if (!(m_elementActive[j])) {
            if (!strcmp((m_elementName[j]).c_str(), "E")) {
                juse = j;
            }
        }
    }
    for (j = 0; j < m_numElemConstraints; j++) {
        if (m_elementActive[j]) {
            if (!strncmp((m_elementName[j]).c_str(), "cn_", 3)) {
                jlose = j;
            }
        }
    }
    for (k = 0; k < m_numSpeciesTot; k++) {
        if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {

            for (j = 0; j < ncTrial; ++j) {
                for (i = 0; i < ncTrial; ++i) {
                    if (i == jlose) {
                        sm[i + j*m_numElemConstraints] = m_formulaMatrix[juse][j];
                    } else {
                        sm[i + j*m_numElemConstraints] = m_formulaMatrix[i][j];
                    }
                }
            }
            for (i = 0; i < m_numRxnTot; ++i) {
                k = m_indexRxnToSpecies[i];
                for (j = 0; j < ncTrial; ++j) {
                    if (j == jlose) {
                        aw[j] = m_formulaMatrix[juse][k];
                    } else {
                        aw[j] = m_formulaMatrix[j][k];
                    }
                }
            }

            j = vcsUtil_gaussj(sm, m_numElemConstraints, ncTrial, aw, 1);
            // j = vcsUtil_mlequ(sm, m_numElemConstraints, ncTrial, aw, 1);
            if (j == 1) {
                plogf("vcs_solve_TP ERROR: mlequ returned an error condition\n");
                return  VCS_FAILED_CONVERGENCE;
            }
            i = k - ncTrial;
            for (j = 0; j < ncTrial; j++) {
                m_stoichCoeffRxnMatrix[i][j] = aw[j];
            }
        }
    }


    /*
     * Calculate the szTmp array for each formation reaction
     */
    for (i = 0; i < m_numRxnTot; i++) {
        double szTmp = 0.0;
        for (j = 0; j < ncTrial; j++) {
            szTmp += fabs(m_stoichCoeffRxnMatrix[i][j]);
        }
        m_scSize[i] = szTmp;
    }


#ifdef DEBUG_MODE
    if (m_debug_print_lvl >= 2) {
        plogf("   ---                Components:");
        for (j = 0; j < ncTrial; j++) {
            plogf("        %3d", j);
        }
        plogf("\n   ---          Components Moles:");
        for (j = 0; j < ncTrial; j++) {
            plogf(" % -10.3E", m_molNumSpecies_old[j]);
        }
        plogf("\n   ---   NonComponent|   Moles  |");
        for (j = 0; j < ncTrial; j++) {
            plogf(" %10.10s", m_speciesName[j].c_str());
        }
        //plogf("|    m_scSize");
        plogf("\n");
        for (i = 0; i < m_numRxnTot; i++) {
            plogf("   --- %3d ", m_indexRxnToSpecies[i]);
            plogf("%-10.10s", m_speciesName[m_indexRxnToSpecies[i]].c_str());
            plogf("|% -10.3E|", m_molNumSpecies_old[m_indexRxnToSpecies[i]]);
            for (j = 0; j < ncTrial; j++) {
                plogf("    %+7.3f", m_stoichCoeffRxnMatrix[i][j]);
            }
            //plogf(" |  %6.2f", m_scSize[i]);
            plogf("\n");
        }


        /*
         *  Manual check on the satisfaction of the reaction matrix's ability
         *  to conserve elements
         */
        double sum;
        double sumMax = -1.0;
        int iMax = -1;
        int jMax = -1;
        size_t n;
        for (i = 0; i < m_numRxnTot; ++i) {
            k = m_indexRxnToSpecies[i];
            for (j = 0; j < ncTrial; ++j) {
                if (j == jlose) {
                    sum = m_formulaMatrix[juse][k];
                    for (n = 0; n < ncTrial; n++) {
                        double numElements = m_formulaMatrix[juse][n];
                        double coeff =  m_stoichCoeffRxnMatrix[i][n];
                        sum += coeff * numElements;
                    }
                } else {
                    sum = m_formulaMatrix[j][k];
                    for (n = 0; n < ncTrial; n++) {
                        double numElements = m_formulaMatrix[j][n];
                        double coeff =  m_stoichCoeffRxnMatrix[i][n];
                        sum += coeff * numElements;
                    }
                }
                if (fabs(sum) > sumMax) {
                    sumMax = fabs(sum);
                    iMax = i;
                    jMax = j;
                    if (j == jlose) {
                        jMax = juse;
                    }
                }
                if (fabs(sum) > 1.0E-6) {
                    printf("we have a prob\n");
                    exit(-1);
                }
            }
        }
        plogf("   ---               largest error in Stoich coeff = %g at rxn = %d ", sumMax, iMax);
        plogf("%-10.10s", m_speciesName[m_indexRxnToSpecies[iMax]].c_str());
        plogf(" element = %d ", jMax);
        plogf("%-5.5s", m_elementName[jMax].c_str());
        plogf("\n");
        plogf("   ");
        for (i=0; i<77; i++) {
            plogf("-");
        }
        plogf("\n");
    }
#endif
    /* **************************************************** */
    /* **** EVALUATE DELTA N VALUES *********************** */
    /* **************************************************** */
    /*
     *       Evaluate the change in gas and liquid total moles
     *       due to reaction vectors, DNG and DNL.
     */

    /*
     *  Zero out the change of Phase Moles array
     */
    vcs_dzero(m_deltaMolNumPhase[0], (NSPECIES0)*(NPHASE0));
    vcs_izero(m_phaseParticipation[0], (NSPECIES0)*(NPHASE0));
    /*
     *  Loop over each reaction, creating the change in Phase Moles
     *  array, m_deltaMolNumPhase[irxn][iphase],
     *  and the phase participation array, PhaseParticipation[irxn][iphase]
     */
    for (irxn = 0; irxn < m_numRxnTot; ++irxn) {
        scrxn_ptr = m_stoichCoeffRxnMatrix[irxn];
        dptr = m_deltaMolNumPhase[irxn];
        kspec = m_indexRxnToSpecies[irxn];
        size_t iph = m_phaseID[kspec];
        int* pp_ptr = m_phaseParticipation[irxn];
        dptr[iph] = 1.0;
        pp_ptr[iph]++;
        for (j = 0; j < ncTrial; ++j) {
            iph = m_phaseID[j];
            if (fabs(scrxn_ptr[j]) <= 1.0e-6) {
                scrxn_ptr[j] = 0.0;
            } else {
                dptr[iph] += scrxn_ptr[j];
                pp_ptr[iph]++;
            }
        }
    }

L_CLEANUP:
    ;
    double tsecond = tickTock.secondsWC();
    m_VCount->Time_basopt += tsecond;
    (m_VCount->Basis_Opts)++;
    return VCS_SUCCESS;
}
/***************************************************************************************/

//  Choose a species to test for the next component
/*
 *  We make the choice based on testing (molNum[i] * spSize[i]) for its maximum value.
 *  Preference for single species phases is also made.
 *
 *  The factors of 1.01 and 1.001 are placed in this routine for a purpose.
 *  The purpose is to ensure that roundoff errors don't influence major decisions.
 *  This means that the optimized and non-optimized versions of the code remain
 *  close to each other.
 *
 *  ( we try to avoid the logic:   a = b
 *                                 if (a > b) { do this }
 *                                 else       { do something else }
 *    because roundoff error makes a difference in the inequality evaluation)
 *
 *  Mole numbers are frequently equal to each other in equilibrium problems
 *  due to constraints. Swaps are only done if there are a 1% difference in the
 *  mole numbers. Of course this logic isn't foolproof.
 *
 * @param molNum  Mole number vector
 * @param j       index into molNum[] that indicates where the search will start from
 *                Previous successful components are swapped into the front of
 *                molNum[].
 * @param n       Length of molNum[]
 */
size_t VCS_SOLVE::vcs_basisOptMax(const double* const molNum, const size_t j,
                                  const size_t n)
{
    size_t largest = j;
    double big = molNum[j] * m_spSize[j] * 1.01;
    if (m_spSize[j] <= 0.0) {
        throw CanteraError("VCS_SOLVE::vcs_basisOptMax",
                           "spSize is nonpositive");
    }
    for (size_t i = j + 1; i < n; ++i) {
        if (m_spSize[i] <= 0.0) {
            throw CanteraError("VCS_SOLVE::vcs_basisOptMax",
                               "spSize is nonpositive");
        }
        bool doSwap = false;
        if (m_SSPhase[j]) {
            doSwap = (molNum[i] * m_spSize[i]) > (big);
            if (!m_SSPhase[i]) {
                if (doSwap) {
                    doSwap = (molNum[i]) > (molNum[largest] * 1.001);
                }
            }
        } else {
            if (m_SSPhase[i]) {
                doSwap = (molNum[i] * m_spSize[i]) > (big);
                if (!doSwap) {
                    doSwap = (molNum[i]) > (molNum[largest] * 1.001);
                }
            } else {
                doSwap = (molNum[i] * m_spSize[i]) > (big);
            }
        }
        if (doSwap) {
            largest = i;
            big = molNum[i] * m_spSize[i] * 1.01;
        }
    }
    return largest;
}
/**********************************************************************************/

// Evaluate the species category for the indicated species
/*
 *  All evaluations are done using the "old" version of the solution.
 *
 *  @param kspec   Species to be evaluated
 *
 *  @return Returns the calculated species type
 */
int VCS_SOLVE::vcs_species_type(const size_t kspec) const
{

    // ---------- Treat special cases first ---------------------


    if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
        return  VCS_SPECIES_INTERFACIALVOLTAGE;
    }

    size_t iph = m_phaseID[kspec];
    int irxn = int(kspec) - int(m_numComponents);
    vcs_VolPhase* VPhase = m_VolPhaseList[iph];
    int phaseExist = VPhase->exists();

    // ---------- Treat zeroed out species first ----------------

    if (m_molNumSpecies_old[kspec] <= 0.0) {

        if (m_tPhaseMoles_old[iph] <= 0.0) {
            if (!m_SSPhase[kspec]) {
                return VCS_SPECIES_ZEROEDMS;
            }
        }

        /*
         *    see if the species has an element
         *     which is so low that species will always be zero
         *
         */
        for (size_t j = 0; j < m_numElemConstraints; ++j) {
            int elType = m_elType[j];
            if (elType == VCS_ELEM_TYPE_ABSPOS) {
                double atomComp = m_formulaMatrix[j][kspec];
                if (atomComp > 0.0) {
                    double maxPermissible = m_elemAbundancesGoal[j] / atomComp;
                    if (maxPermissible < VCS_DELETE_MINORSPECIES_CUTOFF) {
#ifdef DEBUG_MODE
                        if (m_debug_print_lvl >= 2) {
                            plogf("   ---   %s can not be nonzero because"
                                  " needed element %s is zero\n",
                                  m_speciesName[kspec].c_str(), (m_elementName[j]).c_str());
                        }
#endif
                        if (m_SSPhase[kspec]) {
                            return VCS_SPECIES_ZEROEDSS;
                        } else {
                            return VCS_SPECIES_STOICHZERO;
                        }
                    }
                }
            }
        }

        /*
         *     The Gibbs free energy for this species is such that
         *     it will pop back into existence.
         */
        /*
         *     -> An exception to this is if a needed regular element
         *        is also zeroed out. Then, don't pop the phase or the species back into
         *        existence.
         */
        if (irxn >= 0) {
            for (size_t j = 0; j < m_numComponents; ++j) {
                double stoicC = m_stoichCoeffRxnMatrix[irxn][j];
                if (stoicC != 0.0) {
                    double negChangeComp = - stoicC;
                    if (negChangeComp > 0.0) {
                        if (m_molNumSpecies_old[j] < 1.0E-60) {
#ifdef DEBUG_MODE
                            if (m_debug_print_lvl >= 2) {
                                plogf("   ---   %s is prevented from popping into existence because"
                                      " a needed component to be consumed, %s, has a zero mole number\n",
                                      m_speciesName[kspec].c_str(),  m_speciesName[j].c_str());
                            }
#endif
                            if (m_SSPhase[kspec]) {
                                return VCS_SPECIES_ZEROEDSS;
                            } else {
                                return VCS_SPECIES_STOICHZERO;
                            }
                        }
                    } else if (negChangeComp < 0.0) {
                        size_t jph = m_phaseID[j];
                        vcs_VolPhase* jVPhase = m_VolPhaseList[jph];
                        if (jVPhase->exists() <= 0) {
#ifdef DEBUG_MODE
                            if (m_debug_print_lvl >= 2) {
                                plogf("   ---   %s is prevented from popping into existence because"
                                      " a needed component %s is in a zeroed-phase that would be "
                                      "popped into existence at the same time\n",
                                      m_speciesName[kspec].c_str(),  m_speciesName[j].c_str());
                            }
#endif
                            if (m_SSPhase[kspec]) {
                                return VCS_SPECIES_ZEROEDSS;
                            } else {
                                return VCS_SPECIES_STOICHZERO;
                            }
                        }
                    }
                }
            }
        }

        if (irxn >= 0) {
            if (m_deltaGRxn_old[irxn] >= 0.0) {
                /*
                 *  We are here when the species is or should remain zeroed out
                 */
                if (m_SSPhase[kspec]) {
                    return VCS_SPECIES_ZEROEDSS;
                } else {
                    if (phaseExist >= VCS_PHASE_EXIST_YES) {
                        return VCS_SPECIES_ACTIVEBUTZERO;
                    } else if (phaseExist == VCS_PHASE_EXIST_ZEROEDPHASE) {
                        return VCS_SPECIES_ZEROEDPHASE;
                    } else {
                        return VCS_SPECIES_ZEROEDMS;
                    }
                }
            }
        }
        /*
         * If the current phase already exists,
         */
        if (m_tPhaseMoles_old[iph] > 0.0) {
            if (m_SSPhase[kspec]) {
                return VCS_SPECIES_MAJOR;
            } else {
                return VCS_SPECIES_ACTIVEBUTZERO;
            }
        }

        /*
         *     The Gibbs free energy for this species is such that
         *     it will pop back into existence.
         *
         *     -> Set it to a major species in anticipation.
         *      -> note, if we had an estimate for the emerging mole
         *         fraction of the species in the phase, we could do
         *         better here.
         */
        if (m_tPhaseMoles_old[iph] <= 0.0) {
            if (m_SSPhase[kspec]) {
                return VCS_SPECIES_MAJOR;
            } else {
                return VCS_SPECIES_ZEROEDMS;
            }
        }
    }

    // ---------- Treat species with non-zero mole numbers next  ------------

    /*
     *   Always treat species in single species phases as majors if the
     *   phase exists.
     */
    if (m_SSPhase[kspec]) {
        return VCS_SPECIES_MAJOR;
    }

    /*
     *   Check to see whether the current species is a major component
     *   of its phase. If it is, it is a major component. This is consistent
     *   with the above rule about single species phases. A major component
     *   (i.e., a species with a high mole fraction)
     *   in any phase is always treated as a major species
     */
    if (m_molNumSpecies_old[kspec] > (m_tPhaseMoles_old[iph] * 0.001)) {
        return VCS_SPECIES_MAJOR;
    }

    /*
     *   Main check in the loop:
     *      Check to see if there is a component with a mole number that is
     *      within a factor of 100 of the current species.
     *      If there is and that component is not part of a single species
     *      phase and shares a non-zero stoichiometric coefficient, then
     *      the current species is a major species.
     */
    if (irxn < 0) {
        return VCS_SPECIES_MAJOR;
    } else {
        double szAdj = m_scSize[irxn] * std::sqrt((double)m_numRxnTot);
        for (size_t k = 0; k < m_numComponents; ++k) {
            if (!(m_SSPhase[k])) {
                if (m_stoichCoeffRxnMatrix[irxn][k] != 0.0) {
                    if (m_molNumSpecies_old[kspec] * szAdj >= m_molNumSpecies_old[k] * 0.01) {
                        return VCS_SPECIES_MAJOR;
                    }
                }
            }
        }
    }
    return VCS_SPECIES_MINOR;
}
/*****************************************************************************/

//! We calculate the dimensionless chemical potentials of all species
//! in a single phase.
/*!
 *
 * We calculate the dimensionless chemical potentials of all species
 * in a single phase.
 *
 * Note, for multispecies phases which are currently zeroed out,
 * the chemical potential is filled out with the standard chemical
 * potential.
 *
 * For species in multispecies phases whose concentration is zero,
 * we need to set the mole fraction to a very low value.
 * Its chemical potential
 * is then calculated using the VCS_DELETE_MINORSPECIES_CUTOFF concentration
 * to keep numbers positive.
 *
 * Formula:
 * ---------------
 *
 *     Ideal Mixtures:
 *
 *          m_feSpecies(I) = m_SSfeSpecies(I) + ln(z(I)) - ln(m_tPhaseMoles[iph])
 *                         + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase];
 *
 *
 *              ( This is equivalent to the adding the log of the
 *                mole fraction onto the standard chemical
 *                potential. )
 *
 *     Non-Ideal Mixtures:
 *        ActivityConvention = 0:
 *
 *          m_feSpecies(I) = m_SSfeSpecies(I)
 *                         + ln(ActCoeff[I] * z(I)) - ln(m_tPhaseMoles[iph])
 *                         + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase];
 *
 *              ( This is equivalent to the adding the log of the
 *                mole fraction multiplied by the activity coefficient
 *                onto the standard chemical potential. )
 *
 *        ActivityConvention = 1: -> molality activity formulation
 *
 *          m_feSpecies(I) = m_SSfeSpecies(I)
 *                           + ln(ActCoeff[I] * z(I)) - ln(m_tPhaseMoles[iph])
 *                           - ln(Mnaught * m_units)
 *                           + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase];
 *
 *                  note:   m_SSfeSpecies(I) is the molality based standard state.
 *                          However, ActCoeff[I] is the molar based activity coefficient
 *                          We have used the formulas;
 *
 *                     ActCoeff_M[I] =  ActCoeff[I] / Xmol[N]
 *                              where Xmol[N] is the mole fraction of the solvent
 *                               ActCoeff_M[I] is the molality based act coeff.
 *
 *                  note:  This is equivalent to the "normal" molality formulation:
 *
 *                       m_feSpecies(I) = m_SSfeSpecies(I)
 *                                      + ln(ActCoeff_M[I] * m(I))
 *                                      + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase]
 *                                       where m[I] is the molality of the ith solute
 *
 *                                m[I] = Xmol[I] / ( Xmol[N] * Mnaught * m_units)
 *
 *
 *     note:   z(I)/tPhMoles_ptr[iph] = Xmol[i] is the mole fraction
 *                                     of i in the phase.
 *
 *  NOTE:
 *  As per the discussion in vcs_dfe(), for small species where the mole
 *  fraction is small:
 *
 *           z(i) < VCS_DELETE_MINORSPECIES_CUTOFF
 *
 *   The chemical potential is calculated as:
 *
 *          m_feSpecies(I) = m_SSfeSpecies(I)
 *                         + ln(ActCoeff[i](VCS_DELETE_MINORSPECIES_CUTOFF))
 *
 * Input
 * --------
 *   @param iph        Phase to be calculated
 *   @param molNum     molNum[i] is the number of moles of species i
 *                            (VCS species order)
 *   @param do_deleted  Do species that are deleted (default = false)
 *
 * Output
 * -----------
 *   @param ac        Activity coefficients for species in phase
 *                     (VCS species order)
 *   @param mu_i      Dimensionless chemical potentials for phase species
 *                      (VCS species order)
 *
 */
void VCS_SOLVE::vcs_chemPotPhase(const int stateCalc,
                                 const size_t iph, const double* const molNum,
                                 double* const ac, double* const mu_i,
                                 const bool do_deleted)
{

    vcs_VolPhase* Vphase = m_VolPhaseList[iph];
    size_t nkk = Vphase->nSpecies();
    size_t k, kspec;

#ifdef DEBUG_MODE
    //if (m_debug_print_lvl >= 2) {
    // plogf("   --- Subroutine vcs_chemPotPhase called for phase %d\n",
    //     iph);
    //}
#endif
    double tMoles = TPhInertMoles[iph];
    for (k = 0; k < nkk; k++) {
        kspec = Vphase->spGlobalIndexVCS(k);
        tMoles += molNum[kspec];
    }
    double tlogMoles = 0.0;
    if (tMoles > 0.0) {
        tlogMoles = log(tMoles);
    }

    Vphase->setMolesFromVCS(stateCalc, molNum);
    Vphase->sendToVCS_ActCoeff(stateCalc, ac);

    double phi = Vphase->electricPotential();
    double Faraday_phi = m_Faraday_dim * phi;

    for (k = 0; k < nkk; k++) {
        kspec = Vphase->spGlobalIndexVCS(k);
        if (kspec >= m_numComponents) {
            if (!do_deleted &&
                    (m_speciesStatus[kspec] == VCS_SPECIES_DELETED)) {
                continue;
            }
        }
        if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
#ifdef DEBUG_MODE
            if (molNum[kspec] != phi) {
                plogf("We have an inconsistency!\n");
                exit(EXIT_FAILURE);
            }
            if (m_chargeSpecies[kspec] != -1.0) {
                plogf("We have an unexpected situation!\n");
                exit(EXIT_FAILURE);
            }
#endif
            mu_i[kspec] = m_SSfeSpecies[kspec] + m_chargeSpecies[kspec] * Faraday_phi;
        } else {
            if (m_SSPhase[kspec]) {
                mu_i[kspec] = m_SSfeSpecies[kspec] + m_chargeSpecies[kspec] * Faraday_phi;
            } else if (molNum[kspec] <= VCS_DELETE_MINORSPECIES_CUTOFF) {
                mu_i[kspec] = m_SSfeSpecies[kspec] + log(ac[kspec] * VCS_DELETE_MINORSPECIES_CUTOFF)
                              - tlogMoles - m_lnMnaughtSpecies[kspec] + m_chargeSpecies[kspec] * Faraday_phi;
            } else {
                mu_i[kspec] = m_SSfeSpecies[kspec] + log(ac[kspec] * molNum[kspec])
                              - tlogMoles - m_lnMnaughtSpecies[kspec] + m_chargeSpecies[kspec] * Faraday_phi;
            }
        }
    }
}
/*********************************************************************************/

// Calculate the dimensionless chemical potentials of all species or
// of certain groups of species, at a fixed temperature and pressure.
/*
 * We calculate the dimensionless chemical potentials of all species
 * or certain groups of species here, at a fixed temperature and pressure,
 * for the input mole vector z[] in the parameter list.
 * Nondimensionalization is achieved by division by RT.
 *
 * Note, for multispecies phases which are currently zeroed out,
 * the chemical potential is filled out with the standard chemical
 * potential.
 *
 * For species in multispecies phases whose concentration is zero,
 * we need to set the mole fraction to a very low value.
 * Its chemical potential
 * is then calculated using the VCS_DELETE_MINORSPECIES_CUTOFF concentration
 * to keep numbers positive.
 *
 *
 * Formula:
 * ---------------
 *
 *     Ideal Mixtures:
 *
 *          m_feSpecies(I) = m_SSfeSpecies(I) + ln(z(I)) - ln(m_tPhaseMoles[iph])
 *                            + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase];
 *
 *              ( This is equivalent to the adding the log of the
 *                mole fraction onto the standard chemical
 *                potential. )
 *
 *     Non-Ideal Mixtures:  -> molar activity formulation
 *        ActivityConvention = 0:
 *
 *          m_feSpecies(I) = m_SSfeSpecies(I)
 *                           + ln(ActCoeff[I] * z(I)) - ln(m_tPhaseMoles[iph])
 *                            + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase];
 *
 *              ( This is equivalent to the adding the log of the
 *                mole fraction multiplied by the activity coefficient
 *                onto the standard chemical potential. )
 *
 *                 note:   z(I)/tPhMoles_ptr[iph] = Xmol[i] is the mole fraction
 *                                                  of i in the phase.
 *
 *        ActivityConvention = 1: -> molality activity formulation
 *
 *          m_feSpecies(I) = m_SSfeSpecies(I)
 *                           + ln(ActCoeff[I] * z(I)) - ln(m_tPhaseMoles[iph])
 *                           - ln(Mnaught * m_units)
 *                            + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase];
 *
 *                  note:   m_SSfeSpecies(I) is the molality based standard state.
 *                          However, ActCoeff[I] is the molar based activity coefficient
 *                          We have used the formulas;
 *
 *                                ActCoeff_M[I] =  ActCoeff[I] / Xmol[N]
 *                                      where Xmol[N] is the mole fraction of the solvent
 *                                            ActCoeff_M[I] is the molality based act coeff.
 *
 *                  note:   This is equivalent to the "normal" molality formulation below:
 *
 *                              m_feSpecies(I) = m_SSfeSpecies(I)
 *                                             + ln(ActCoeff_M[I] * m(I))
 *                                             + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase]
 *                                       where m[I] is the molality of the ith solute
 *
 *
 *                                m[I] = Xmol[I] / ( Xmol[N] * Mnaught * m_units)
 *
 *
 *  Handling of Small Species:
 * ------------------------------
 *  As per the discussion above, for small species where the mole
 *  fraction
 *
 *           z(i) < VCS_DELETE_MINORSPECIES_CUTOFF
 *
 *   The chemical potential is calculated as:
 *
 *         m_feSpecies(I) = m_SSfeSpecies(I) + ln(ActCoeff[i](VCS_DELETE_MINORSPECIES_CUTOFF))
 *
 *   Handling of "Species" Representing Interfacial Voltages
 *  ---------------------------------------------------------
 *
 *      These species have species types of  VCS_SPECIES_TYPE_INTERFACIALVOLTAGE
 *      The chemical potentials for these "species" refer to electrons in
 *      metal electrodes. They have the following formula
 *
 *          m_feSpecies(I) = m_SSfeSpecies(I) - F z[I] / RT
 *
 *      F is Faraday's constant.
 *      R = gas constant
 *      T = temperature
 *      V = potential of the interface = phi_electrode - phi_solution
 *
 *      For these species, the solution vector unknown, z[I], is V, the phase voltage,
 *      in volts.
 *
 * Input
 * --------
 * @param ll     Determine which group of species gets updated
 *     ll =  0: Calculate for all species
 *        <  0: calculate for components and for major non-components
 *           1: calculate for components and for minor non-components
 *
 * @param lbot    Restricts the calculation of the chemical potential
 *                to the species between LBOT <= i < LTOP. Usually
 *                 LBOT and LTOP will be equal to 0 and MR, respectively.
 * @param ltop    Top value of the loops
 *
 * @param molNum   molNum[i]   : Number of moles of species i
 *                   -> This can either be the old solution vector
 *                      or the new solution vector depending upon the
 *                      stateCalc value
 *
 * @param stateCalc   Determines whether z is old or new or tmp:
 *            VCS_STATECALC_NEW: Use the tentative values for the total number of
 *                               moles in the phases, i.e., use TG1 instead of TG etc.
 *            VCS_STATECALC_OLD: Use the base values of the total number of
 *                               moles in each system.
 *
 *  Also needed:
 *    m_SSfeSpecies[kspec] : standard state chemical potentials. These are the
 *                           chemical potentials of the standard states at
 *                           the same T and P as the solution.
 *     tg     : Total Number of moles in the phase.
 */
void VCS_SOLVE::vcs_dfe(const int stateCalc,
                        const int ll, const size_t lbot, const size_t ltop)
{
    size_t l1, l2, iph, kspec, irxn;
    size_t iphase;
    double* tPhMoles_ptr=0;
    double* actCoeff_ptr=0;
    double* tlogMoles=0;
    vcs_VolPhase* Vphase;

    double* feSpecies=0;
    double* molNum=0;
    if (stateCalc == VCS_STATECALC_OLD) {
        feSpecies = VCS_DATA_PTR(m_feSpecies_old);
        tPhMoles_ptr = VCS_DATA_PTR(m_tPhaseMoles_old);
        actCoeff_ptr = VCS_DATA_PTR(m_actCoeffSpecies_old);
        molNum = VCS_DATA_PTR(m_molNumSpecies_old);
    } else if (stateCalc == VCS_STATECALC_NEW) {
        feSpecies = VCS_DATA_PTR(m_feSpecies_new);
        tPhMoles_ptr = VCS_DATA_PTR(m_tPhaseMoles_new);
        actCoeff_ptr = VCS_DATA_PTR(m_actCoeffSpecies_new);
        molNum = VCS_DATA_PTR(m_molNumSpecies_new);
    }
#ifdef DEBUG_MODE
    else {
        plogf("vcs_dfe: wrong stateCalc value");
        plogf("   --- Subroutine vcs_dfe called with bad stateCalc value: %d", stateCalc);
        plogendl();
        exit(EXIT_FAILURE);
    }
#endif

#ifdef DEBUG_MODE
    if (m_unitsState == VCS_DIMENSIONAL_G) {
        printf("vcs_dfe: called with wrong units state\n");
        exit(EXIT_FAILURE);
    }
#endif

#ifdef DEBUG_MODE
    if (m_debug_print_lvl >= 2) {
        if (ll == 0) {
            if (lbot != 0) {
                plogf("   --- Subroutine vcs_dfe called for one species: ");
                plogf("%-12.12s", m_speciesName[lbot].c_str());
            } else {
                plogf("   --- Subroutine vcs_dfe called for all species");
            }
        } else if (ll > 0) {
            plogf("   --- Subroutine vcs_dfe called for components and minors");
        } else {
            plogf("   --- Subroutine vcs_dfe called for components and majors");
        }
        if (stateCalc == VCS_STATECALC_NEW) {
            plogf(" using tentative solution");
        }
        plogendl();
    }
#endif

    tlogMoles = VCS_DATA_PTR(m_TmpPhase);
    /*
     * Might as well recalculate the phase mole vector
     * and compare to the stored one. They should be correct.
     */
    double* tPhInertMoles = VCS_DATA_PTR(TPhInertMoles);
    for (iph = 0; iph < m_numPhases; iph++) {
        tlogMoles[iph] = tPhInertMoles[iph];

    }
    for (kspec = 0; kspec < m_numSpeciesTot; kspec++) {
        if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
            iph = m_phaseID[kspec];
            tlogMoles[iph] += molNum[kspec];
        }
    }
#ifdef DEBUG_MODE
    for (iph = 0; iph < m_numPhases; iph++) {
        if (! vcs_doubleEqual(tlogMoles[iph], tPhMoles_ptr[iph])) {
            plogf("phase Moles may be off, iph = %d, %20.14g %20.14g \n",
                  iph, tlogMoles[iph], tPhMoles_ptr[iph]);
            exit(EXIT_FAILURE);
        }
    }
#endif
    vcs_dzero(tlogMoles, m_numPhases);
    for (iph = 0; iph < m_numPhases; iph++) {
        if (tPhMoles_ptr[iph] > 0.0) {
            tlogMoles[iph] = log(tPhMoles_ptr[iph]);
        }
    }


    if (ll != 0) {
        l1 = lbot;
        l2 = m_numComponents;
    } else {
        l1 = lbot;
        l2 = ltop;
    }

    /*
     *  Calculate activity coefficients for all phases that are
     *  not current. Here we also trigger an update check for each
     *  VolPhase to see if its mole numbers are current with vcs
     */
    for (iphase = 0; iphase < m_numPhases; iphase++) {
        Vphase = m_VolPhaseList[iphase];
        Vphase->updateFromVCS_MoleNumbers(stateCalc);
        if (!Vphase->m_singleSpecies) {
            Vphase->sendToVCS_ActCoeff(stateCalc, VCS_DATA_PTR(actCoeff_ptr));
        }
        m_phasePhi[iphase] = Vphase->electricPotential();
    }
    /* ************************************************************** */
    /* **** ALL SPECIES, OR COMPONENTS ****************************** */
    /* ************************************************************** */
    /*
     *     Do all of the species when LL = 0. Then we are done for the routine
     *     When LL ne 0., just do the initial components. We will then
     *     finish up below with loops over either the major noncomponent
     *     species or the minor noncomponent species.
     */
    for (kspec = l1; kspec < l2; ++kspec) {
        iphase = m_phaseID[kspec];
        if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
#ifdef DEBUG_MODE
            if (molNum[kspec] != m_phasePhi[iphase]) {
                plogf("We have an inconsistency!\n");
                exit(EXIT_FAILURE);
            }
            if (m_chargeSpecies[kspec] != -1.0) {
                plogf("We have an unexpected situation!\n");
                exit(EXIT_FAILURE);
            }
#endif
            feSpecies[kspec] = m_SSfeSpecies[kspec]
                               + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
        } else {
            if (m_SSPhase[kspec]) {
                feSpecies[kspec] = m_SSfeSpecies[kspec]
                                   + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
            } else if ((m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDMS) ||
                       (m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDPHASE)) {
                feSpecies[kspec] = m_SSfeSpecies[kspec] - m_lnMnaughtSpecies[kspec]
                                   + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
            } else {
                if (molNum[kspec] <= VCS_DELETE_MINORSPECIES_CUTOFF) {
                    iph = m_phaseID[kspec];
                    if (tPhMoles_ptr[iph] > 0.0) {
                        feSpecies[kspec] = m_SSfeSpecies[kspec]
                                           + log(actCoeff_ptr[kspec] * VCS_DELETE_MINORSPECIES_CUTOFF)
                                           - tlogMoles[m_phaseID[kspec]] - m_lnMnaughtSpecies[kspec]
                                           + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
                    } else {
                        feSpecies[kspec] = m_SSfeSpecies[kspec] - m_lnMnaughtSpecies[kspec]
                                           + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
                    }
                } else {
                    feSpecies[kspec] = m_SSfeSpecies[kspec]
                                       + log(actCoeff_ptr[kspec] * molNum[kspec])
                                       - tlogMoles[m_phaseID[kspec]] - m_lnMnaughtSpecies[kspec]
                                       + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
                }
            }
        }
    }
    /* ************************************************ */
    /* **** MAJORS ONLY ******************************* */
    /* ************************************************ */
    if (ll < 0) {
        for (irxn = 0; irxn < m_numRxnRdc; ++irxn) {
            kspec = m_indexRxnToSpecies[irxn];
            if (m_speciesStatus[kspec] != VCS_SPECIES_MINOR) {
                iphase = m_phaseID[kspec];
                if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
#ifdef DEBUG_MODE
                    if (molNum[kspec] != m_phasePhi[iphase]) {
                        plogf("We have an inconsistency!\n");
                        exit(EXIT_FAILURE);
                    }
                    if (m_chargeSpecies[kspec] != -1.0) {
                        plogf("We have an unexpected situation!\n");
                        exit(EXIT_FAILURE);
                    }
#endif
                    feSpecies[kspec] = m_SSfeSpecies[kspec]
                                       + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
                } else {
                    if (m_SSPhase[kspec]) {
                        feSpecies[kspec] = m_SSfeSpecies[kspec]
                                           + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
                    } else if ((m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDMS) ||
                               (m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDPHASE)) {
                        feSpecies[kspec] = m_SSfeSpecies[kspec] - m_lnMnaughtSpecies[kspec]
                                           + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
                    } else {
                        if (molNum[kspec] <= VCS_DELETE_MINORSPECIES_CUTOFF) {
                            iph = m_phaseID[kspec];
                            if (tPhMoles_ptr[iph] > 0.0) {
                                feSpecies[kspec] = m_SSfeSpecies[kspec]
                                                   + log(actCoeff_ptr[kspec] * VCS_DELETE_MINORSPECIES_CUTOFF)
                                                   - tlogMoles[m_phaseID[kspec]] - m_lnMnaughtSpecies[kspec]
                                                   + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
                            } else {
                                feSpecies[kspec] = m_SSfeSpecies[kspec]  - m_lnMnaughtSpecies[kspec]
                                                   + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
                            }
                        } else {
                            feSpecies[kspec] = m_SSfeSpecies[kspec]
                                               + log(actCoeff_ptr[kspec] * molNum[kspec])
                                               - tlogMoles[m_phaseID[kspec]] - m_lnMnaughtSpecies[kspec]
                                               + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
                        }
                    }
                }
            }
        }
        /* ************************************************ */
        /* **** MINORS ONLY ******************************* */
        /* ************************************************ */
    } else if (ll > 0) {
        for (irxn = 0; irxn < m_numRxnRdc; ++irxn) {
            kspec = m_indexRxnToSpecies[irxn];
            if (m_speciesStatus[kspec] == VCS_SPECIES_MINOR) {
                iphase = m_phaseID[kspec];
                if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
#ifdef DEBUG_MODE
                    if (molNum[kspec] != m_phasePhi[iphase]) {
                        plogf("We have an inconsistency!\n");
                        exit(EXIT_FAILURE);
                    }
                    if (m_chargeSpecies[kspec] != -1.0) {
                        plogf("We have an unexpected situation!\n");
                        exit(EXIT_FAILURE);
                    }
#endif
                    feSpecies[kspec] = m_SSfeSpecies[kspec]
                                       + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase]; ;
                } else {
                    if (m_SSPhase[kspec]) {
                        feSpecies[kspec] = m_SSfeSpecies[kspec]
                                           + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
                    } else if ((m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDMS) ||
                               (m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDPHASE)) {
                        feSpecies[kspec] = m_SSfeSpecies[kspec] - m_lnMnaughtSpecies[kspec]
                                           + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
                    } else {
                        if (molNum[kspec] <= VCS_DELETE_MINORSPECIES_CUTOFF) {
                            iph = m_phaseID[kspec];
                            if (tPhMoles_ptr[iph] > 0.0) {
                                feSpecies[kspec] = m_SSfeSpecies[kspec]
                                                   + log(actCoeff_ptr[kspec] * VCS_DELETE_MINORSPECIES_CUTOFF)
                                                   - tlogMoles[m_phaseID[kspec]] - m_lnMnaughtSpecies[kspec]
                                                   + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
                            } else {
                                feSpecies[kspec] = m_SSfeSpecies[kspec] - m_lnMnaughtSpecies[kspec]
                                                   + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
                            }
                        } else {
                            feSpecies[kspec] = m_SSfeSpecies[kspec]
                                               + log(actCoeff_ptr[kspec] * molNum[kspec])
                                               - tlogMoles[m_phaseID[kspec]] - m_lnMnaughtSpecies[kspec]
                                               + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
                        }
                    }
                }
            }
        }
    }


}
//====================================================================================================================

void  VCS_SOLVE::vcs_printSpeciesChemPot(const int stateCalc) const
{
    double mfValue = 1.0;
    bool zeroedPhase = false;
    size_t kspec;

    const double* molNum = VCS_DATA_PTR(m_molNumSpecies_old);
    const double* actCoeff_ptr = VCS_DATA_PTR(m_actCoeffSpecies_old);
    if (stateCalc == VCS_STATECALC_NEW) {
        actCoeff_ptr = VCS_DATA_PTR(m_actCoeffSpecies_new);
        molNum = VCS_DATA_PTR(m_molNumSpecies_new);
    }

    double* tMoles = VCS_DATA_PTR(m_TmpPhase);
    const double* tPhInertMoles = VCS_DATA_PTR(TPhInertMoles);
    for (size_t iph = 0; iph < m_numPhases; iph++) {
        tMoles[iph] = tPhInertMoles[iph];
    }
    for (kspec = 0; kspec < m_numSpeciesTot; kspec++) {
        if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
            size_t iph = m_phaseID[kspec];
            tMoles[iph] += molNum[kspec];
        }
    }

    double RT = m_temperature * Cantera::GasConstant;
    printf("   ---  CHEMICAL POT TABLE (J/kmol) Name PhID     MolFR     ChemoSS   "
           "   logMF       Gamma       Elect       extra       ElectrChem\n");
    printf("   ");
    vcs_print_line("-", 132);

    for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
        mfValue = 1.0;
        size_t iphase = m_phaseID[kspec];
        const vcs_VolPhase* Vphase = m_VolPhaseList[iphase];
        if ((m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDMS)    ||
                (m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDPHASE) ||
                (m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDSS)) {
            zeroedPhase = true;
        } else {
            zeroedPhase = false;
        }
        if (tMoles[iphase] > 0.0) {
            if (molNum[kspec] <= VCS_DELETE_MINORSPECIES_CUTOFF) {
                mfValue = VCS_DELETE_MINORSPECIES_CUTOFF / tMoles[iphase];
            } else {
                mfValue = molNum[kspec]/tMoles[iphase];
            }
        } else {
            size_t klocal = m_speciesLocalPhaseIndex[kspec];
            mfValue = Vphase->moleFraction(klocal);
        }
        double volts = Vphase->electricPotential();
        double elect = m_chargeSpecies[kspec] * m_Faraday_dim * volts;
        double comb = - m_lnMnaughtSpecies[kspec];
        double total = (m_SSfeSpecies[kspec] + log(mfValue) + elect + log(actCoeff_ptr[kspec]) + comb);

        if (zeroedPhase) {
            printf("   --- ** zp *** ");
        } else {
            printf("   ---           ");
        }
        printf("%-24.24s", m_speciesName[kspec].c_str());
        printf(" %-3s", Cantera::int2str(iphase).c_str());
        printf(" % -12.4e", mfValue);
        printf(" % -12.4e", m_SSfeSpecies[kspec] * RT);
        printf(" % -12.4e", log(mfValue) * RT);
        printf(" % -12.4e", log(actCoeff_ptr[kspec]) * RT);
        printf(" % -12.4e", elect * RT);
        printf(" % -12.4e", comb * RT);
        printf(" % -12.4e\n", total *RT);
    }
    printf("   ");
    vcs_print_line("-", 132);
}
/*****************************************************************************/

#ifdef DEBUG_MODE
//! Print out and check the elemental abundance vector
void VCS_SOLVE::prneav() const
{
    size_t j;
    bool kerr;
    std::vector<double> eav(m_numElemConstraints, 0.0);

    for (j = 0; j < m_numElemConstraints; ++j) {
        for (size_t i = 0; i < m_numSpeciesTot; ++i) {
            if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
                eav[j] += m_formulaMatrix[j][i] * m_molNumSpecies_old[i];
            }
        }
    }
    kerr = false;
    plogf("--------------------------------------------------");
    plogf("ELEMENT ABUNDANCE VECTOR:\n");
    plogf(" Element   Now      Orignal      Deviation          Type\n");
    for (j = 0; j < m_numElemConstraints; ++j) {
        plogf("  ");
        plogf("%-2.2s", (m_elementName[j]).c_str());
        plogf(" = %15.6E     %15.6E     %15.6E  %3d\n",
              eav[j], m_elemAbundancesGoal[j], eav[j] - m_elemAbundancesGoal[j], m_elType[j]);
        if (m_elemAbundancesGoal[j] != 0.) {
            if (fabs(eav[j] - m_elemAbundancesGoal[j]) > m_elemAbundancesGoal[j] * 5.0e-9) {
                kerr = true;
            }
        } else {
            if (fabs(eav[j]) > 1.0e-10) {
                kerr = true;
            }
        }
    }
    if (kerr) {
        plogf("Element abundance check failure\n");
    }
    plogf("--------------------------------------------------");
    plogendl();
}
#endif
/*****************************************************************************/

double VCS_SOLVE::l2normdg(double dgLocal[]) const
{
    double tmp;
    size_t irxn;
    if (m_numRxnRdc <= 0) {
        return 0.0;
    }
    for (irxn = 0, tmp = 0.0; irxn < m_numRxnRdc; ++irxn) {
        size_t kspec = irxn + m_numComponents;
        if (m_speciesStatus[kspec] == VCS_SPECIES_MAJOR || m_speciesStatus[kspec] == VCS_SPECIES_MINOR ||
                dgLocal[irxn] < 0.0) {
            if (m_speciesStatus[kspec] != VCS_SPECIES_ZEROEDMS) {
                tmp += dgLocal[irxn] * dgLocal[irxn];
            }
        }
    }
    return std::sqrt(tmp / m_numRxnRdc);
}
/*****************************************************************************/

//  Calculates the total number of moles of species in all phases.
/*
 *  Calculates the total number of moles in all phases.
 *  Reconciles Phase existence flags with total moles in each phase.
 */
double VCS_SOLVE::vcs_tmoles()
{
    double sum;
    vcs_VolPhase* Vphase;
    for (size_t i = 0; i < m_numPhases; i++) {
        m_tPhaseMoles_old[i] = TPhInertMoles[i];
    }
    for (size_t i = 0; i < m_numSpeciesTot; i++) {
        if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_MOLNUM) {
            m_tPhaseMoles_old[m_phaseID[i]] += m_molNumSpecies_old[i];
        }
    }
    sum = 0.0;
    for (size_t i = 0; i < m_numPhases; i++) {
        sum += m_tPhaseMoles_old[i];
        Vphase = m_VolPhaseList[i];
        // Took out because we aren't updating mole fractions in Vphase
        // Vphase->TMoles = m_tPhaseMoles_old[i];
        if (m_tPhaseMoles_old[i] == 0.0) {
            Vphase->setTotalMoles(0.0);
        } else {
            Vphase->setTotalMoles(m_tPhaseMoles_old[i]);
        }
    }
    m_totalMolNum = sum;
    return m_totalMolNum;
}

#ifdef DEBUG_MODE
void VCS_SOLVE::check_tmoles() const
{
    size_t i;
    double sum = 0.0;
    for (i = 0; i < m_numPhases; i++) {
        double m_tPhaseMoles_old_a = TPhInertMoles[i];

        for (size_t k = 0; k < m_numSpeciesTot; k++) {
            if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_MOLNUM) {
                if (m_phaseID[k] == i) {
                    m_tPhaseMoles_old_a += m_molNumSpecies_old[k];
                }
            }
        }
        sum += m_tPhaseMoles_old_a;

        double denom = m_tPhaseMoles_old[i]+ m_tPhaseMoles_old_a + 1.0E-19;
        if (!vcs_doubleEqual(m_tPhaseMoles_old[i]/denom, m_tPhaseMoles_old_a/denom)) {
            plogf("check_tmoles: we have found a problem with phase  %d: %20.15g, %20.15g\n",
                  i, m_tPhaseMoles_old[i], m_tPhaseMoles_old_a);
        }
    }
}
#endif
/*****************************************************************************/

void VCS_SOLVE::vcs_updateVP(const int vcsState)
{
    vcs_VolPhase* Vphase;
    for (size_t i = 0; i < m_numPhases; i++) {
        Vphase = m_VolPhaseList[i];
        if (vcsState == VCS_STATECALC_OLD) {
            Vphase->setMolesFromVCSCheck(VCS_STATECALC_OLD,
                                         VCS_DATA_PTR(m_molNumSpecies_old),
                                         VCS_DATA_PTR(m_tPhaseMoles_old));
        } else if (vcsState == VCS_STATECALC_NEW) {
            Vphase->setMolesFromVCSCheck(VCS_STATECALC_NEW,
                                         VCS_DATA_PTR(m_molNumSpecies_new),
                                         VCS_DATA_PTR(m_tPhaseMoles_new));
        }
#ifdef DEBUG_MODE
        else {
            plogf("vcs_updateVP ERROR: wrong stateCalc value: %d", vcsState);
            plogendl();
            exit(EXIT_FAILURE);
        }
#endif
    }
}

//  This routine evaluates the species type for all species
/*
 *        kspec
 *                              1 -> Major species  VCS_SPECIES_MAJOR
 *                              0 -> Minor species  VCS_SPECIES_MINOR
 *                             -1 -> The species lies in a multicomponent phase
 *                                   that exists. Its concentration is currently
 *                                   very low, necessitating a different method
 *                                   of calculation.
 *                                   -  VCS_SPECIES_ZEROEDPHASE
 *                             -2 -> The species lies in a multicomponent phase
 *                                   which currently doesn't exist.
 *                                   Its concentration is currently zero.
 *                                   -  VCS_SPECIES_ZEROEDMS
 *                             -3 -> Species lies in a single-species phase which
 *                                   is currently zeroed out.
 *                                   - VCS_SPECIES_ZEROEDSS
 *                             -4 -> Species has such a small mole fraction it is
 *                                   deleted even though its phase may possibly exist.
 *                                   The species is believed to have such a small
 *                                   mole fraction that it best to throw the
 *                                   calculation of it out.  It will be added back in
 *                                   at the end of the calculation.
 *                                   - VCS_SPECIES_DELETED
 *                             -5 ->  Species refers to an electron in the metal
 *                                    The unknown is equal to the interfacial voltage
 *                                    drop across the interface on the SHE (standard
 *                                    hydrogen electrode) scale (volts).
 *                                     - VCS_SPECIES_INTERFACIALVOLTAGE
 *                             -6 ->  Species lies in a multicomponent phase that
 *                                    is zeroed atm  and will stay deleted due to a
 *                                    choice from a higher level.
 *                                    These species will formally always have zero
 *                                    mole numbers in the solution vector.
 *                                      - VCS_SPECIES_ZEROEDPHASE
 *                             -7 ->  The species lies in a multicomponent phase which
 *                                    currently does exist.  Its concentration is currently
 *                                    identically zero, though the phase exists. Note, this
 *                                    is a temporary condition that exists at the start
 *                                    of an equilibrium problem.
 *                                    The species is soon "birthed" or "deleted".
 *                                    - VCS_SPECIES_ACTIVEBUTZERO
 *                             -8 ->  The species lies in a multicomponent phase which
 *                                     currently does exist.  Its concentration is currently
 *                                    identically zero, though the phase exists. This is
 *                                    a permanent condition due to stoich constraints
 *                                    - VCS_SPECIES_STOICHZERO
 *
 */
bool VCS_SOLVE::vcs_evaluate_speciesType()
{
    bool allMinorZeroedSpecies;

    m_numRxnMinorZeroed = 0;
#ifdef DEBUG_MODE
    if (m_debug_print_lvl >= 2) {
        plogf("  --- Species Status decision is reevaluated: All species are minor except for:\n");
    } else if (m_debug_print_lvl >= 5) {
        plogf("  --- Species Status decision is reevaluated");
        plogendl();
    }
#endif
    for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
        m_speciesStatus[kspec] = vcs_species_type(kspec);
#ifdef DEBUG_MODE
        if (m_debug_print_lvl >= 5) {
            plogf("  --- %-16s: ", m_speciesName[kspec].c_str());
            if (kspec < m_numComponents) {
                plogf("(COMP) ");
            } else {
                plogf("       ");
            }
            plogf(" %10.3g ", m_molNumSpecies_old[kspec]);
            const char* sString = vcs_speciesType_string(m_speciesStatus[kspec], 100);
            plogf("%s\n", sString);

        } else if (m_debug_print_lvl >= 2) {
            if (m_speciesStatus[kspec] != VCS_SPECIES_MINOR) {
                switch (m_speciesStatus[kspec]) {
                case VCS_SPECIES_COMPONENT:
                    break;
                case VCS_SPECIES_MAJOR:
                    plogf("  ---      Major Species          : %-s\n", m_speciesName[kspec].c_str());
                    break;
                case VCS_SPECIES_ZEROEDPHASE:
                    plogf("  ---      Purposely Zeroed-Phase Species (not in problem): %-s\n",
                          m_speciesName[kspec].c_str());
                    break;
                case VCS_SPECIES_ZEROEDMS:
                    plogf("  ---      Zeroed-MS Phase Species: %-s\n", m_speciesName[kspec].c_str());
                    break;
                case VCS_SPECIES_ZEROEDSS:
                    plogf("  ---      Zeroed-SS Phase Species: %-s\n", m_speciesName[kspec].c_str());
                    break;
                case VCS_SPECIES_DELETED:
                    plogf("  ---      Deleted-Small Species  : %-s\n", m_speciesName[kspec].c_str());
                    break;
                case VCS_SPECIES_ACTIVEBUTZERO:
                    plogf("  ---      Zeroed Species in an active MS phase (tmp): %-s\n",
                          m_speciesName[kspec].c_str());
                    break;
                case VCS_SPECIES_STOICHZERO:
                    plogf("  ---     Zeroed Species in an active MS phase (Stoich Constraint): %-s\n",
                          m_speciesName[kspec].c_str());
                    break;
                case VCS_SPECIES_INTERFACIALVOLTAGE:
                    plogf("  ---      InterfaceVoltage Species: %-s\n", m_speciesName[kspec].c_str());
                    break;
                default:
                    plogf("  --- Unknown type - ERROR %d\n", m_speciesStatus[kspec]);
                    plogendl();
                    exit(EXIT_FAILURE);
                }
            }
        }
#endif
        if (kspec >= m_numComponents) {
            if (m_speciesStatus[kspec] != VCS_SPECIES_MAJOR) {
                ++m_numRxnMinorZeroed;
            }
        }
    }
#ifdef DEBUG_MODE
    if (m_debug_print_lvl >= 2) {
        plogf("  ---");
        plogendl();
    }
#endif

    allMinorZeroedSpecies = (m_numRxnMinorZeroed >= m_numRxnRdc);

    return allMinorZeroedSpecies;
}
/*****************************************************************************/

//    Switch rows and columns of a square matrix
/*
 *  Switches the row and column of a matrix.
 *  So that after
 *
 *     J[k1][j] = J_old[k2][j]  and J[j][k1] = J_old[j][k2]
 *     J[k2][j] = J_old[k1][j]  and J[j][k2] = J_old[j][k1]
 *
 *  @param Jac         Double pointer to the Jacobian
 *  @param k1          first row/column value to be switched
 *  @param k2          second row/column value to be switched
 */
void VCS_SOLVE::vcs_switch2D(double* const* const Jac,
                             const size_t k1, const size_t k2) const
{
    size_t i;
    if (k1 == k2) {
        return;
    }
    for (i = 0; i < m_numSpeciesTot; i++) {
        std::swap(Jac[k1][i], Jac[k2][i]);
    }
    for (i = 0; i < m_numSpeciesTot; i++) {
        std::swap(Jac[i][k1], Jac[i][k2]);
    }
}
/*****************************************************************************/

static void print_space(size_t num)
{
    size_t j;
    for (j = 0; j < num; j++) {
        plogf(" ");
    }
}
/********************************************************************************/

// This subroutine calculates reaction free energy changes for
// all noncomponent formation reactions.
/*
 *  Formation reactions are
 *  reactions which create each noncomponent species from the component
 *  species. m_stoichCoeffRxnMatrix[irxn][jcomp]  are the stoichiometric
 *  coefficients for these  reactions. A stoichiometric coefficient of
 *  one is assumed for species irxn in this reaction.
 *
 *  INPUT
 *  @param l
 *    L < 0   :  Calculate reactions corresponding to
 *               major noncomponent and zeroed species only
 *    L = 0   :  Do all noncomponent reactions, i, between
 *               0 <= i < irxnl
 *    L > 0   :  Calculate reactions corresponding to
 *               minor noncomponent and zeroed species only
 *
 *    @param doDeleted   Do deleted species
 *    @param stateCalc   Calculate deltaG corresponding to either old or new
 *                       free energies
 *
 *    Note we special case one important issue.
 *    If the component has zero moles, then we do not
 *    allow deltaG < 0.0 for formation reactions which
 *    would lead to the loss of more of that same component.
 *    This dG < 0.0 condition feeds back into the algorithm in several
 *    places, and leads to a infinite loop in at least one case.
 */
void VCS_SOLVE::vcs_deltag(const int l, const bool doDeleted,
                           const int vcsState, const bool alterZeroedPhases)
{
    size_t iph;
    size_t irxn, kspec;
    bool lneed;
    double* dtmp_ptr;
    int icase = 0;
    size_t irxnl = m_numRxnRdc;
    if (doDeleted) {
        irxnl = m_numRxnTot;
    }

    double* deltaGRxn;
    double* feSpecies;
    double* molNumSpecies;
    double* actCoeffSpecies;
    if (vcsState == VCS_STATECALC_NEW) {
        deltaGRxn     = VCS_DATA_PTR(m_deltaGRxn_new);
        feSpecies     = VCS_DATA_PTR(m_feSpecies_new);
        molNumSpecies = VCS_DATA_PTR(m_molNumSpecies_new);
        actCoeffSpecies = VCS_DATA_PTR(m_actCoeffSpecies_new);
    } else if (vcsState == VCS_STATECALC_OLD) {
        deltaGRxn     = VCS_DATA_PTR(m_deltaGRxn_old);
        feSpecies     = VCS_DATA_PTR(m_feSpecies_old);
        molNumSpecies = VCS_DATA_PTR(m_molNumSpecies_old);
        actCoeffSpecies = VCS_DATA_PTR(m_actCoeffSpecies_old);
    } else {
        printf("Error\n");
        exit(EXIT_FAILURE);
    }

#ifdef DEBUG_MODE
    if (m_debug_print_lvl >= 2) {
        plogf("   --- Subroutine vcs_deltag called for ");
        if (l < 0) {
            plogf("major noncomponents\n");
        } else if (l == 0) {
            plogf("all noncomponents\n");
        } else {
            plogf("minor noncomponents\n");
        }
    }
#endif
    /* ************************************************* */
    /* **** MAJORS and ZEROED SPECIES ONLY ************* */
    /* ************************************************* */
    if (l < 0) {
        for (irxn = 0; irxn < m_numRxnRdc; ++irxn) {
            kspec = irxn + m_numComponents;
            if (m_speciesStatus[kspec] != VCS_SPECIES_MINOR) {
                icase = 0;
                deltaGRxn[irxn] = feSpecies[m_indexRxnToSpecies[irxn]];
                dtmp_ptr = m_stoichCoeffRxnMatrix[irxn];
                for (kspec = 0; kspec < m_numComponents; ++kspec) {
                    deltaGRxn[irxn] += dtmp_ptr[kspec] * feSpecies[kspec];
                    if (molNumSpecies[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF && dtmp_ptr[kspec] < 0.0) {
                        icase = 1;
                    }
                }
                if (icase) {
                    deltaGRxn[irxn] = std::max(0.0, deltaGRxn[irxn]);
                }
            }
        }
    } else if (l == 0) {
        /* ************************************************* */
        /* **** ALL REACTIONS ****************************** */
        /* ************************************************* */
        for (irxn = 0; irxn < irxnl; ++irxn) {
            icase = 0;
            deltaGRxn[irxn] = feSpecies[m_indexRxnToSpecies[irxn]];
            dtmp_ptr = m_stoichCoeffRxnMatrix[irxn];
            for (kspec = 0; kspec < m_numComponents; ++kspec) {
                deltaGRxn[irxn] += dtmp_ptr[kspec] * feSpecies[kspec];
                if (molNumSpecies[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF &&
                        dtmp_ptr[kspec] < 0.0) {
                    icase = 1;
                }
            }
            if (icase) {
                deltaGRxn[irxn] = std::max(0.0, deltaGRxn[irxn]);
            }
        }
    } else {
        /* ************************************************* */
        /* **** MINORS AND ZEROED SPECIES ****************** */
        /* ************************************************* */
        for (irxn = 0; irxn < m_numRxnRdc; ++irxn) {
            kspec = irxn + m_numComponents;
            if (m_speciesStatus[kspec] <= VCS_SPECIES_MINOR) {
                icase = 0;
                deltaGRxn[irxn] = feSpecies[m_indexRxnToSpecies[irxn]];
                dtmp_ptr = m_stoichCoeffRxnMatrix[irxn];
                for (kspec = 0; kspec < m_numComponents; ++kspec) {
                    deltaGRxn[irxn] += dtmp_ptr[kspec] * feSpecies[kspec];
                    if (m_molNumSpecies_old[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF &&
                            dtmp_ptr[kspec] < 0.0) {
                        icase = 1;
                    }
                }
                if (icase) {
                    deltaGRxn[irxn] = std::max(0.0, deltaGRxn[irxn]);
                }
            }
        }
    }
    /* ************************************************* */
    /* **** MULTISPECIES PHASES WITH ZERO MOLES ******** */
    /* ************************************************* */
    /*
     *    Massage the free energies for species with zero mole fractions
     *  in multispecies phases.  This section implements the
     *  Equation 3.8-5 in Smith and Missen, p.59.
     *  A multispecies phase will exist iff
     *           1 < sum_i(exp(-dg_i)/AC_i)
     *  If DG is negative then that species wants to be reintroduced into
     *  the calculation.
     *  For small dg_i, the expression below becomes:
     *      1 - sum_i(exp(-dg_i)/AC_i) ~ sum_i((dg_i-1)/AC_i)  + 1
     *
     *  So, what we are doing here is equalizing all DG's in a multispecies
     *  phase whose total mole number has already been zeroed out.
     *  It must have to do with the case where a complete multispecies
     *  phase is currently zeroed out. In that case, when one species
     *  in that phase has a negative DG, then the phase should kick in.
     *  This code section will cause that to happen, because a negative
     *  DG will dominate the calculation of SDEL. Then, DG(I) for all
     *  species in that phase will be forced to be equal and negative.
     *  Thus, all species in that phase will come into being at the
     *  same time.
     *
     *  HKM -> The ratio of mole fractions at the reinstatement
     *         time should be equal to the normalized weighting
     *         of exp(-dg_i) / AC_i. This should be implemented.
     *
     *  HKM -> There is circular logic here. ActCoeff depends on the
     *         mole fractions of a phase that does not exist. In actuality
     *         the proto-mole fractions should be selected from the
     *         solution of a nonlinear problem with NsPhase unknowns
     *
     *              X_i = exp(-dg[irxn]) / ActCoeff_i / denom
     *
     *              where
     *               denom = sum_i[  exp(-dg[irxn]) / ActCoeff_i  ]
     *
     *         This can probably be solved by successive iteration.
     *         This should be implemented.
     */
    //alterZeroedPhases = false;
    if (alterZeroedPhases  && false) {
        for (iph = 0; iph < m_numPhases; iph++) {
            lneed = false;
            vcs_VolPhase* Vphase = m_VolPhaseList[iph];
            if (! Vphase->m_singleSpecies) {
                double sum = 0.0;
                for (size_t k = 0; k < Vphase->nSpecies(); k++) {
                    kspec = Vphase->spGlobalIndexVCS(k);
                    if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
                        sum += molNumSpecies[kspec];
                    }
                    if (sum > 0.0) {
                        break;
                    }
                }
                if (sum == 0.0) {
                    lneed = true;
                }
            }

            if (lneed) {
                double poly = 0.0;
                for (size_t k = 0; k < Vphase->nSpecies(); k++) {
                    kspec = Vphase->spGlobalIndexVCS(k);
                    // We may need to look at deltaGRxn for components!
                    if (kspec >= m_numComponents) {
                        irxn = kspec - m_numComponents;
                        if (deltaGRxn[irxn] >  50.0) {
                            deltaGRxn[irxn] =  50.0;
                        }
                        if (deltaGRxn[irxn] < -50.0) {
                            deltaGRxn[irxn] = -50.0;
                        }
                        poly += exp(-deltaGRxn[irxn])/actCoeffSpecies[kspec];
                    }
                }
                /*
                 *      Calculate deltaGRxn[] for each species in a zeroed multispecies phase.
                 *      All of the m_deltaGRxn_new[]'s will be equal. If deltaGRxn[] is
                 *      negative, then the phase will come back into existence.
                 */
                for (size_t k = 0; k < Vphase->nSpecies(); k++) {
                    kspec = Vphase->spGlobalIndexVCS(k);
                    if (kspec >= m_numComponents) {
                        irxn = kspec - m_numComponents;
                        deltaGRxn[irxn] = 1.0 - poly;
                    }
                }

            }
        }
    }


#ifdef DEBUG_NOT
    for (irxn = 0; irxn < m_numRxnRdc; ++irxn) {
        checkFinite(deltaGRxn[irxn]);
    }
#endif
}
//====================================================================================================================
void  VCS_SOLVE::vcs_printDeltaG(const int stateCalc)
{
    size_t j;
    double* deltaGRxn = VCS_DATA_PTR(m_deltaGRxn_old);
    double* feSpecies = VCS_DATA_PTR(m_feSpecies_old);
    double* molNumSpecies = VCS_DATA_PTR(m_molNumSpecies_old);
    const double* tPhMoles_ptr = VCS_DATA_PTR(m_tPhaseMoles_old);
    const double* actCoeff_ptr = VCS_DATA_PTR(m_actCoeffSpecies_old);
    if (stateCalc == VCS_STATECALC_NEW) {
        deltaGRxn     = VCS_DATA_PTR(m_deltaGRxn_new);
        feSpecies     = VCS_DATA_PTR(m_feSpecies_new);
        molNumSpecies = VCS_DATA_PTR(m_molNumSpecies_new);
        actCoeff_ptr = VCS_DATA_PTR(m_actCoeffSpecies_new);
        tPhMoles_ptr = VCS_DATA_PTR(m_tPhaseMoles_new);
    }
    double RT = m_temperature * Cantera::GasConstant;
    bool zeroedPhase = false;
    if (m_debug_print_lvl >= 2) {
        plogf("   --- DELTA_G TABLE  Components:");
        for (j = 0; j < m_numComponents; j++) {
            plogf("     %3d  ", j);
        }
        plogf("\n   ---          Components Moles:");
        for (j = 0; j < m_numComponents; j++) {
            plogf("%10.3g", m_molNumSpecies_old[j]);
        }
        plogf("\n   ---   NonComponent|   Moles  |       ");
        for (j = 0; j < m_numComponents; j++) {
            plogf("%-10.10s", m_speciesName[j].c_str());
        }
        //plogf("|    m_scSize");
        plogf("\n");
        for (size_t i = 0; i < m_numRxnTot; i++) {
            plogf("   --- %3d ", m_indexRxnToSpecies[i]);
            plogf("%-10.10s", m_speciesName[m_indexRxnToSpecies[i]].c_str());
            plogf("|%10.3g|", m_molNumSpecies_old[m_indexRxnToSpecies[i]]);
            for (j = 0; j < m_numComponents; j++) {
                plogf("     %6.2f", m_stoichCoeffRxnMatrix[i][j]);
            }
            //plogf(" |  %6.2f", m_scSize[i]);
            plogf("\n");
        }
        plogf("   ");
        for (int i=0; i<77; i++) {
            plogf("-");
        }
        plogf("\n");
    }

    printf("   --- DeltaG Table (J/kmol) Name       PhID   MoleNum      MolFR  "
           "  ElectrChemStar ElectrChem    DeltaGStar   DeltaG(Pred)  Stability\n");
    printf("   ");
    vcs_print_line("-", 132);

    for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) {

        size_t irxn = kspec -  m_numComponents;

        double mfValue = 1.0;
        size_t iphase = m_phaseID[kspec];
        const vcs_VolPhase* Vphase = m_VolPhaseList[iphase];
        if ((m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDMS)    ||
                (m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDPHASE) ||
                (m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDSS)) {
            zeroedPhase = true;
        } else {
            zeroedPhase = false;
        }
        if (tPhMoles_ptr[iphase] > 0.0) {
            if (molNumSpecies[kspec] <= VCS_DELETE_MINORSPECIES_CUTOFF) {
                mfValue = VCS_DELETE_MINORSPECIES_CUTOFF / tPhMoles_ptr[iphase];
            } else {
                mfValue = molNumSpecies[kspec] / tPhMoles_ptr[iphase];
            }
        } else {
            size_t klocal = m_speciesLocalPhaseIndex[kspec];
            mfValue = Vphase->moleFraction(klocal);
        }
        if (zeroedPhase) {
            printf("   --- ** zp *** ");
        } else {
            printf("   ---           ");
        }
        double feFull =  feSpecies[kspec];
        if ((m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDMS)    ||
                (m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDPHASE)) {
            feFull += log(actCoeff_ptr[kspec]) + log(mfValue);
        }
        printf("%-24.24s", m_speciesName[kspec].c_str());
        printf(" %-3s", Cantera::int2str(iphase).c_str());
        printf(" % -12.4e", molNumSpecies[kspec]);
        printf(" % -12.4e", mfValue);
        printf(" % -12.4e", feSpecies[kspec] * RT);
        printf(" % -12.4e", feFull * RT);
        if (irxn != npos) {
            printf(" % -12.4e", deltaGRxn[irxn] * RT);
            printf(" % -12.4e", (deltaGRxn[irxn] + feFull - feSpecies[kspec]) * RT);

            if (deltaGRxn[irxn] < 0.0) {
                if (molNumSpecies[kspec] > 0.0) {
                    printf("   growing");
                } else {
                    printf("    stable");
                }
            } else if (deltaGRxn[irxn] > 0.0) {
                if (molNumSpecies[kspec] > 0.0) {
                    printf(" shrinking");
                } else {
                    printf("  unstable");
                }
            } else {
                printf(" balanced");
            }
        }

        printf(" \n");
    }

    printf("   ");
    vcs_print_line("-", 132);

}
//====================================================================================================================
//    Calculate deltag of formation for all species in a single phase.
/*
 *     Calculate deltag of formation for all species in a single
 *     phase. It is assumed that the fe[] is up to date for all species.
 *     However, if the phase is currently zeroed out, a subproblem
 *     is calculated to solve for AC[i] and pseudo-X[i] for that
 *     phase.
 *
 *   @param iphase       phase index of the phase to be calculated
 *   @param doDeleted    boolean indicating whether to do deleted species or not
 *   @param stateCalc    integer describing which set of free energies
 *                       to use and where to stick the results.
 *
 *    NOTE: this is currently not used used anywhere. It may be in the future?
 */
void VCS_SOLVE::vcs_deltag_Phase(const size_t iphase, const bool doDeleted,
                                 const int stateCalc, const bool alterZeroedPhases)
{
    size_t iph;
    size_t  irxn, kspec, kcomp;
    double* dtmp_ptr;

    double* feSpecies=0;
    double* deltaGRxn=0;
    double* actCoeffSpecies=0;
    if (stateCalc == VCS_STATECALC_NEW) {
        feSpecies = VCS_DATA_PTR(m_feSpecies_new);
        deltaGRxn = VCS_DATA_PTR(m_deltaGRxn_new);
        actCoeffSpecies = VCS_DATA_PTR(m_actCoeffSpecies_new);
    } else if (stateCalc == VCS_STATECALC_OLD) {
        feSpecies = VCS_DATA_PTR(m_feSpecies_old);
        deltaGRxn = VCS_DATA_PTR(m_deltaGRxn_old);
        actCoeffSpecies = VCS_DATA_PTR(m_actCoeffSpecies_old);
    }
#ifdef DEBUG_MODE
    else {
        plogf("vcs_deltag_Phase: we shouldn't be here\n");
        plogendl();
        exit(EXIT_FAILURE);
    }
#endif

    size_t irxnl = m_numRxnRdc;
    if (doDeleted) {
        irxnl = m_numRxnTot;
    }
    vcs_VolPhase* vPhase = m_VolPhaseList[iphase];

#ifdef DEBUG_MODE
    if (m_debug_print_lvl >= 2) {
        plogf("   --- Subroutine vcs_deltag_Phase called for phase %d\n",
              iphase);
    }
#endif

    /*
     * Single species Phase
     */
    if (vPhase->m_singleSpecies) {
        kspec = vPhase->spGlobalIndexVCS(0);
#ifdef DEBUG_MODE
        if (iphase != m_phaseID[kspec]) {
            plogf("vcs_deltag_Phase index error\n");
            exit(EXIT_FAILURE);
        }
#endif
        if (kspec >= m_numComponents) {
            irxn = kspec - m_numComponents;
            deltaGRxn[irxn] = feSpecies[kspec];
            dtmp_ptr = m_stoichCoeffRxnMatrix[irxn];
            for (kcomp = 0; kcomp < m_numComponents; ++kcomp) {
                deltaGRxn[irxn] += dtmp_ptr[kcomp] * feSpecies[kcomp];
            }
        }
    }
    /*
     * Multispecies Phase
     */
    else {
        bool zeroedPhase = true;

        for (irxn = 0; irxn < irxnl; ++irxn) {
            kspec = m_indexRxnToSpecies[irxn];
            if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
                iph = m_phaseID[kspec];
                if (iph == iphase) {
                    if (m_molNumSpecies_old[kspec] > 0.0) {
                        zeroedPhase = false;
                    }
                    deltaGRxn[irxn] = feSpecies[kspec];
                    dtmp_ptr = m_stoichCoeffRxnMatrix[irxn];
                    for (kcomp = 0; kcomp < m_numComponents; ++kcomp) {
                        deltaGRxn[irxn] += dtmp_ptr[kcomp] * feSpecies[kcomp];
                    }
                }
            }
        }

        /*
         * special section for zeroed phases
         */
        /* ************************************************* */
        /* **** MULTISPECIES PHASES WITH ZERO MOLES************ */
        /* ************************************************* */
        /*
         *    Massage the free energies for species with zero mole fractions
         *  in multispecies phases.  This section implements the
         *  Equation 3.8-5 in Smith and Missen, p.59.
         *  A multispecies phase will exist iff
         *           1 < sum_i(exp(-dg_i)/AC_i)
         *  If DG is negative then that species wants to be reintroduced into
         *  the calculation.
         *  For small dg_i, the expression below becomes:
         *      1 - sum_i(exp(-dg_i)/AC_i) ~ sum_i((dg_i-1)/AC_i)  + 1
         *
         *
         *  HKM -> The ratio of mole fractions at the reinstatement
         *         time should be equal to the normalized weighting
         *         of exp(-dg_i) / AC_i. This should be implemented.
         *
         *  HKM -> There is circular logic here. ActCoeff depends on the
         *         mole fractions of a phase that does not exist. In actuality
         *         the proto-mole fractions should be selected from the
         *         solution of a nonlinear problem with NsPhase unknowns
         *
         *              X_i = exp(-dg[irxn]) / ActCoeff_i / denom
         *
         *              where
         *               denom = sum_i[  exp(-dg[irxn]) / ActCoeff_i  ]
         *
         *         This can probably be solved by successive iteration.
         *         This should be implemented.
         */
        /*
         *      Calculate dg[] for each species in a zeroed multispecies phase.
         *      All of the dg[]'s will be equal. If dg[] is negative, then
         *      the phase will come back into existence.
         */
        if (alterZeroedPhases) {
            if (zeroedPhase) {
                double phaseDG = 1.0;
                for (irxn = 0; irxn < irxnl; ++irxn) {
                    kspec = m_indexRxnToSpecies[irxn];
                    iph = m_phaseID[kspec];
                    if (iph == iphase) {
                        if (deltaGRxn[irxn] >  50.0) {
                            deltaGRxn[irxn] =  50.0;
                        }
                        if (deltaGRxn[irxn] < -50.0) {
                            deltaGRxn[irxn] = -50.0;
                        }
                        phaseDG -= exp(-deltaGRxn[irxn])/actCoeffSpecies[kspec];
                    }
                }
                /*
                 * Overwrite the individual dg's with the phase DG.
                 */
                for (irxn = 0; irxn < irxnl; ++irxn) {
                    kspec = m_indexRxnToSpecies[irxn];
                    iph = m_phaseID[kspec];
                    if (iph == iphase) {
                        deltaGRxn[irxn] = 1.0 - phaseDG;
                    }
                }
            }
        }
    }
}
/****************************************************************************/

//  Swaps the indices for all of the global data for two species, k1
//  and k2.
/*
 *  @param  ifunc:  If true, switch the species data and the noncomponent reaction
 *                  data. This must be called for a non-component species only.
 *                  If false, switch the species data only. Typically, we use this
 *                  option when determining the component species and at the
 *                  end of the calculation, when we want to return unscrambled
 *                  results. All rxn data will be out-of-date.
 *
 *  @param k1        First species index
 *
 *  @param k2        Second species index
 */
void VCS_SOLVE::vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k2)
{
    size_t i1, i2, iph, kp1, kp2;
    vcs_VolPhase* pv1, *pv2;
    if (k1 == k2) {
        return;
    }
#ifdef DEBUG_MODE
    if (k1 > (m_numSpeciesTot - 1) ||
            k2 > (m_numSpeciesTot - 1)) {
        plogf("vcs_switch_pos: ifunc = 0: inappropriate args: %d %d\n",
              k1, k2);
    }
#endif
    /*
     *     Handle the index pointer in the phase structures first
     */
    pv1 = m_VolPhaseList[m_phaseID[k1]];
    pv2 = m_VolPhaseList[m_phaseID[k2]];

    kp1 = m_speciesLocalPhaseIndex[k1];
    kp2 = m_speciesLocalPhaseIndex[k2];
#ifdef DEBUG_MODE
    if (pv1->spGlobalIndexVCS(kp1) != k1) {
        plogf("Indexing error in program\n");
        exit(EXIT_FAILURE);
    }
    if (pv2->spGlobalIndexVCS(kp2) != k2) {
        plogf("Indexing error in program\n");
        exit(EXIT_FAILURE);
    }
#endif
    pv1->setSpGlobalIndexVCS(kp1, k2);
    pv2->setSpGlobalIndexVCS(kp2, k1);
    //pv1->IndSpecies[kp1] = k2;
    //pv2->IndSpecies[kp2] = k1;
    std::swap(m_speciesName[k1], m_speciesName[k2]);
    std::swap(m_molNumSpecies_old[k1], m_molNumSpecies_old[k2]);
    std::swap(m_speciesUnknownType[k1], m_speciesUnknownType[k2]);
    std::swap(m_molNumSpecies_new[k1], m_molNumSpecies_new[k2]);
    std::swap(m_SSfeSpecies[k1], m_SSfeSpecies[k2]);
    std::swap(m_spSize[k1], m_spSize[k2]);
    std::swap(m_deltaMolNumSpecies[k1], m_deltaMolNumSpecies[k2]);
    std::swap(m_feSpecies_old[k1], m_feSpecies_old[k2]);
    std::swap(m_feSpecies_new[k1], m_feSpecies_new[k2]);
    std::swap(m_SSPhase[k1], m_SSPhase[k2]);
    std::swap(m_phaseID[k1], m_phaseID[k2]);
    std::swap(m_speciesMapIndex[k1], m_speciesMapIndex[k2]);
    std::swap(m_speciesLocalPhaseIndex[k1], m_speciesLocalPhaseIndex[k2]);
    std::swap(m_actConventionSpecies[k1], m_actConventionSpecies[k2]);
    std::swap(m_lnMnaughtSpecies[k1], m_lnMnaughtSpecies[k2]);
    std::swap(m_actCoeffSpecies_new[k1], m_actCoeffSpecies_new[k2]);
    std::swap(m_actCoeffSpecies_old[k1], m_actCoeffSpecies_old[k2]);
    std::swap(m_wtSpecies[k1], m_wtSpecies[k2]);
    std::swap(m_chargeSpecies[k1], m_chargeSpecies[k2]);
    std::swap(m_speciesThermoList[k1], m_speciesThermoList[k2]);
    std::swap(m_PMVolumeSpecies[k1], m_PMVolumeSpecies[k2]);

    for (size_t j = 0; j < m_numElemConstraints; ++j) {
        std::swap(m_formulaMatrix[j][k1], m_formulaMatrix[j][k2]);
    }
    if (m_useActCoeffJac) {
        vcs_switch2D(m_np_dLnActCoeffdMolNum.baseDataAddr(), k1, k2);
    }
    std::swap(m_speciesStatus[k1], m_speciesStatus[k2]);
    /*
     *     Handle the index pointer in the phase structures
     */


    if (ifunc) {
        /*
         * Find the Rxn indices corresponding to the two species
         */
        i1 = k1 - m_numComponents;
        i2 = k2 - m_numComponents;
#ifdef DEBUG_MODE
        if (i1 > (m_numRxnTot - 1) ||
                i2 > (m_numRxnTot - 1)) {
            plogf("switch_pos: ifunc = 1: inappropriate noncomp values: %d %d\n",
                  i1 , i2);
        }
#endif
        for (size_t j = 0; j < m_numComponents; ++j) {
            std::swap(m_stoichCoeffRxnMatrix[i1][j], m_stoichCoeffRxnMatrix[i2][j]);
        }
        std::swap(m_scSize[i1], m_scSize[i2]);
        for (iph = 0; iph < m_numPhases; iph++) {
            std::swap(m_deltaMolNumPhase[i1][iph], m_deltaMolNumPhase[i2][iph]);
            std::swap(m_phaseParticipation[i1][iph],
                      m_phaseParticipation[i2][iph]);
        }
        std::swap(m_deltaGRxn_new[i1],  m_deltaGRxn_new[i2]);
        std::swap(m_deltaGRxn_old[i1], m_deltaGRxn_old[i2]);
        std::swap(m_deltaGRxn_tmp[i1], m_deltaGRxn_tmp[i2]);

        /*
         *   We don't want to swap ir[], because the values of ir should
         *   stay the same after the swap
         *
         * vcs_isw(ir, i1, i2);
         */
    }
}
/*******************************************************************************/


//     Birth guess returns the number of moles of a species
//     that is coming back to life.
/*
 *     Birth guess returns the number of moles of a species
 *     that is coming back to life.
 *     Note, this routine is not applicable if the whole phase is coming
 *     back to life, not just one species in that phase.
 *
 *  Do a minor alt calculation. But, cap the mole numbers at
 *  1.0E-15.
 *  For SS phases use VCS_DELETE_SPECIES_CUTOFF * 100.
 *
 *  The routine makes sure the guess doesn't reduce the concentration
 *  of a component by more than 1/3. Note this may mean that
 *  the vlaue coming back from this routine is zero or a
 *  very small number.
 *
 *
 * @param kspec   Species number that is coming back to life
 *
 * @return      Returns the number of kmol that the species should
 *              have.
 */
double VCS_SOLVE::vcs_birthGuess(const int kspec)
{
    size_t irxn = kspec - m_numComponents;
    double dx = 0.0;
    if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
        return dx;
    }
    double w_kspec = VCS_DELETE_MINORSPECIES_CUTOFF;
#ifdef DEBUG_MODE
    // Check to make sure that species is zero in the solution vector
    // If it isn't, we don't know what's happening
    if (m_molNumSpecies_old[kspec] != 0.0) {
        w_kspec = 0.0;
        plogf("vcs_birthGuess:: we shouldn't be here\n");
        exit(EXIT_FAILURE);
    }
#endif
    int ss = m_SSPhase[kspec];
    if (!ss) {
        /*
         * Logic to handle species in multiple species phases
         *    we cap the moles here at 1.0E-15 kmol.
         */
        bool soldel_ret;
#ifdef DEBUG_MODE
        char ANOTE[32];
        double dxm = vcs_minor_alt_calc(kspec, irxn, &soldel_ret, ANOTE);
#else
        double dxm = vcs_minor_alt_calc(kspec, irxn, &soldel_ret);
#endif
        dx = w_kspec + dxm;
        if (dx > 1.0E-15) {
            dx = 1.0E-15;
        }
    } else {
        /*
         * Logic to handle single species phases
         *      There is no real way to estimate the moles. So
         *      we set it to a small number.
         */
        dx = 1.0E-30;
    }

    /*
     * Check to see if the current value of the components
     * allow the dx just estimated.
     * If we are in danger of zeroing a component,
     * only go 1/3 the way to zeroing the component with
     * this dx. Note, this may mean that dx= 0 coming
     * back from this routine. This evaluation should
     * be respected.
     */
    double* sc_irxn = m_stoichCoeffRxnMatrix[irxn];
    for (size_t j = 0; j < m_numComponents; ++j) {
        // Only loop over element constraints that involve positive def. constraints
        if (m_speciesUnknownType[j] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
            if (m_molNumSpecies_old[j] > 0.0) {
                double tmp = sc_irxn[j] * dx;
                if (3.0*(-tmp) > m_molNumSpecies_old[j]) {
                    dx = std::min(dx, - 0.3333* m_molNumSpecies_old[j] / sc_irxn[j]);
                }
            }
            if (m_molNumSpecies_old[j] <= 0.0) {
                if (sc_irxn[j] < 0.0) {
                    dx = 0.0;
                }
            }
        }
    }
    return dx;
}
/*******************************************************************************/

void VCS_SOLVE::vcs_setFlagsVolPhases(const bool upToDate, const int stateCalc)
{
    vcs_VolPhase* Vphase;
    if (!upToDate) {
        for (size_t iph = 0; iph < m_numPhases; iph++) {
            Vphase = m_VolPhaseList[iph];
            Vphase->setMolesOutOfDate(stateCalc);
        }
    } else {
        for (size_t iph = 0; iph < m_numPhases; iph++) {
            Vphase = m_VolPhaseList[iph];
            Vphase->setMolesCurrent(stateCalc);
        }
    }
}
/*******************************************************************************/

void VCS_SOLVE::vcs_setFlagsVolPhase(const size_t iph, const bool upToDate,
                                     const int stateCalc)
{
    vcs_VolPhase* Vphase = m_VolPhaseList[iph];
    if (!upToDate) {
        Vphase->setMolesOutOfDate(stateCalc);
    } else {
        Vphase->setMolesCurrent(stateCalc);
    }
}
/*******************************************************************************/

// Update all underlying vcs_VolPhase objects
/*
 *  Update the mole numbers and the phase voltages of all phases in the
 *  vcs problem
 *
 *  @param stateCalc Location of the update (either VCS_STATECALC_NEW or
 *         VCS_STATECALC_OLD).
 */
void VCS_SOLVE::vcs_updateMolNumVolPhases(const int stateCalc)
{
    vcs_VolPhase* Vphase;
    for (size_t iph = 0; iph < m_numPhases; iph++) {
        Vphase = m_VolPhaseList[iph];
        Vphase->updateFromVCS_MoleNumbers(stateCalc);
    }
}
/*******************************************************************************/

}