diff --git a/Cantera/src/thermo/GibbsExcessVPSSTP.cpp b/Cantera/src/thermo/GibbsExcessVPSSTP.cpp
index fd7e460c3..13ba5ce75 100644
--- a/Cantera/src/thermo/GibbsExcessVPSSTP.cpp
+++ b/Cantera/src/thermo/GibbsExcessVPSSTP.cpp
@@ -65,7 +65,8 @@ namespace Cantera {
     moleFractions_       = b.moleFractions_;
     lnActCoeff_Scaled_   = b.lnActCoeff_Scaled_;
     dlnActCoeffdT_Scaled_   = b.dlnActCoeffdT_Scaled_;
-    dlnActCoeffdlnC_Scaled_ = b.dlnActCoeffdlnC_Scaled_;
+    dlnActCoeffdlnX_Scaled_ = b.dlnActCoeffdlnX_Scaled_;
+    dlnActCoeffdlnN_Scaled_ = b.dlnActCoeffdlnN_Scaled_;
     m_pp                 = b.m_pp;
 
     return *this;
@@ -156,7 +157,9 @@ namespace Cantera {
   }
 
   void GibbsExcessVPSSTP::calcDensity() {
-    double *vbar = &m_pp[0];
+    doublereal* vbar = NULL;
+    vbar = new doublereal[m_kk];
+    //    double *vbar = &m_pp[0];
     getPartialMolarVolumes(vbar);
    
     doublereal vtotal = 0.0;
@@ -165,6 +168,7 @@ namespace Cantera {
     }
     doublereal dd = meanMolecularWeight() / vtotal;
     State::setDensity(dd);
+    delete [] vbar;
   }
 
   void GibbsExcessVPSSTP::setState_TP(doublereal t, doublereal p) {
@@ -320,7 +324,8 @@ namespace Cantera {
     moleFractions_.resize(m_kk);
     lnActCoeff_Scaled_.resize(m_kk);
     dlnActCoeffdT_Scaled_.resize(m_kk);
-    dlnActCoeffdlnC_Scaled_.resize(m_kk);
+    dlnActCoeffdlnX_Scaled_.resize(m_kk);
+    dlnActCoeffdlnN_Scaled_.resize(m_kk);
     m_pp.resize(m_kk);
   }
 
diff --git a/Cantera/src/thermo/GibbsExcessVPSSTP.h b/Cantera/src/thermo/GibbsExcessVPSSTP.h
index 185c2be27..3e8b7cd10 100644
--- a/Cantera/src/thermo/GibbsExcessVPSSTP.h
+++ b/Cantera/src/thermo/GibbsExcessVPSSTP.h
@@ -301,6 +301,23 @@ namespace Cantera {
     virtual void getdlnActCoeffdT(doublereal *dlnActCoeffdT) const {
       err("getdlnActCoeffdT");
     }
+
+    //! Get the array of change in the log activity coefficients w.r.t. change in state (change temp, change mole fractions)
+    /*!
+     * This function is a virtual class, but it first appears in GibbsExcessVPSSTP
+     * class and derived classes from GibbsExcessVPSSTP.
+     *
+     * This function is a virtual method.  For ideal mixtures 
+     * (unity activity coefficients), this can gradX/X.  
+     *
+     * @param dT    Input of temperature change
+     * @param dX    Input vector of changes in mole fraction. length = m_kk
+     * @param dlnActCoeff    Output vector of derivatives of the 
+     *                         log Activity Coefficients. length = m_kk
+     */
+    virtual void getdlnActCoeff(const doublereal dT, const doublereal * const dX, doublereal *dlnActCoeff) const {
+      err("getdlnActCoeff");
+    }
  
     //! Get the array of log concentration-like derivatives of the 
     //! log activity coefficients
@@ -317,11 +334,33 @@ namespace Cantera {
      *
      *  units = dimensionless
      *
-     * @param dlnActCoeffdlnC    Output vector of derivatives of the 
+     * @param dlnActCoeffdlnN    Output vector of derivatives of the 
      *                         log Activity Coefficients. length = m_kk
      */
-    virtual void getdlnActCoeffdlnC(doublereal *dlnActCoeffdlnC) const {
-      err("getdlnActCoeffdlnC");
+    virtual void getdlnActCoeffdlnN(doublereal *dlnActCoeffdlnN) const {
+      err("getdlnActCoeffdlnN");
+    }
+
+    //! Get the array of log concentration-like derivatives of the 
+    //! log activity coefficients
+    /*!
+     * This function is a virtual method.  For ideal mixtures 
+     * (unity activity coefficients), this can return zero.  
+     * Implementations should take the derivative of the 
+     * logarithm of the activity coefficient with respect to the 
+     * logarithm of the concentration-like variable (i.e. number of moles in
+     * in a unit volume. ) that represents the standard state.  
+     * This quantity is to be used in conjunction with derivatives of 
+     * that concentration-like variable when the derivative of the chemical 
+     * potential is taken.  
+     *
+     *  units = dimensionless
+     *
+     * @param dlnActCoeffdlnX    Output vector of derivatives of the 
+     *                         log Activity Coefficients. length = m_kk
+     */
+    virtual void getdlnActCoeffdlnX(doublereal *dlnActCoeffdlnX) const {
+      err("getdlnActCoeffdlnX");
     }
  
     //@}
@@ -562,7 +601,12 @@ namespace Cantera {
     //! Storage for the current derivative values of the 
     //! gradients with respect to logarithm of the mole fraction of the 
     //! log of theactivity coefficients of the species
-    mutable std::vector<doublereal> dlnActCoeffdlnC_Scaled_;
+    mutable std::vector<doublereal> dlnActCoeffdlnN_Scaled_;
+
+    //! Storage for the current derivative values of the 
+    //! gradients with respect to logarithm of the mole fraction of the 
+    //! log of theactivity coefficients of the species
+    mutable std::vector<doublereal> dlnActCoeffdlnX_Scaled_;
 
     //! Temporary storage space that is fair game
     mutable std::vector<doublereal> m_pp;
diff --git a/Cantera/src/thermo/IonsFromNeutralVPSSTP.cpp b/Cantera/src/thermo/IonsFromNeutralVPSSTP.cpp
index df90bac16..4e95c4f30 100644
--- a/Cantera/src/thermo/IonsFromNeutralVPSSTP.cpp
+++ b/Cantera/src/thermo/IonsFromNeutralVPSSTP.cpp
@@ -197,7 +197,8 @@ namespace Cantera {
     muNeutralMolecule_          = b.muNeutralMolecule_;
     gammaNeutralMolecule_       = b.gammaNeutralMolecule_;
     dlnActCoeffdT_NeutralMolecule_ = b.dlnActCoeffdT_NeutralMolecule_;
-    dlnActCoeffdlnC_NeutralMolecule_ = b.dlnActCoeffdlnC_NeutralMolecule_;
+    dlnActCoeffdlnX_NeutralMolecule_ = b.dlnActCoeffdlnX_NeutralMolecule_;
+    dlnActCoeffdlnN_NeutralMolecule_ = b.dlnActCoeffdlnN_NeutralMolecule_;
 
     return *this;
   }
@@ -332,6 +333,13 @@ namespace Cantera {
   void IonsFromNeutralVPSSTP::getActivityConcentrations(doublereal* c) const {
     getActivities(c);
   }
+
+  void IonsFromNeutralVPSSTP::getDissociationCoeffs(vector_fp& coeffs,vector_fp& charges){
+    coeffs = fm_neutralMolec_ions_;
+    charges = m_speciesCharge;
+    //for ( int k = 0; k < fm_neutralMolec_ions_[k]; k++ )
+    //  coeffs.push_back(fm_neutralMolec_ions_[k]);
+  }
  
   // Return the standard concentration for the kth species
   /*
@@ -598,22 +606,47 @@ namespace Cantera {
      *
      *  units = dimensionless
      *
-     * @param dlnActCoeffdlnC    Output vector of log(mole fraction)  
+     * @param dlnActCoeffdlnX    Output vector of log(mole fraction)  
      *                 derivatives of the log Activity Coefficients.
      *                 length = m_kk
      */
-  void IonsFromNeutralVPSSTP::getdlnActCoeffdlnC(doublereal *dlnActCoeffdlnC) const {
+  void IonsFromNeutralVPSSTP::getdlnActCoeffdlnX(doublereal *dlnActCoeffdlnX) const {
     s_update_lnActCoeff();
-    s_update_dlnActCoeff_dlnC();
+    s_update_dlnActCoeff_dlnX();
 
     for (int k = 0; k < m_kk; k++) {
-      dlnActCoeffdlnC[k] = dlnActCoeffdlnC_Scaled_[k];
+      dlnActCoeffdlnX[k] = dlnActCoeffdlnX_Scaled_[k];
+    }
+  }  
+  
+  //! Get the array of log concentration-like derivatives of the 
+    //! log activity coefficients
+    /*!
+     * This function is a virtual method.  For ideal mixtures 
+     * (unity activity coefficients), this can return zero.  
+     * Implementations should take the derivative of the 
+     * logarithm of the activity coefficient with respect to the 
+     * logarithm of the concentration-like variable (i.e. moles)
+     * that represents the standard state.  
+     * This quantity is to be used in conjunction with derivatives of 
+     * that concentration-like variable when the derivative of the chemical 
+     * potential is taken.  
+     *
+     *  units = dimensionless
+     *
+     * @param dlnActCoeffdlnN    Output vector of log(mole fraction)  
+     *                 derivatives of the log Activity Coefficients.
+     *                 length = m_kk
+     */
+  void IonsFromNeutralVPSSTP::getdlnActCoeffdlnN(doublereal *dlnActCoeffdlnN) const {
+    s_update_lnActCoeff();
+    s_update_dlnActCoeff_dlnN();
+
+    for (int k = 0; k < m_kk; k++) {
+      dlnActCoeffdlnN[k] = dlnActCoeffdlnN_Scaled_[k];
     }
   }
     
-
-
-
   // This is temporary. We will get rid of this
   void IonsFromNeutralVPSSTP::setTemperature(const doublereal temp) {
     double p = pressure();
@@ -717,7 +750,7 @@ namespace Cantera {
     for (k = 0; k < m_kk; k++) {
      sum += moleFractions_[k];
     }
-    if (fabs(sum) > 1.0E-11) {
+    if (fabs(sum) > 1.0E-11)  {
       throw CanteraError("IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions", 
 			 "molefracts don't sum to one: " + fp2str(sum));
     }
@@ -791,7 +824,129 @@ namespace Cantera {
         sum += NeutralMolecMoleFractions_[k];
       }
       for (k = 0; k < numNeutralMoleculeSpecies_; k++) {
-        NeutralMolecMoleFractions_[k] /= sum;
+	NeutralMolecMoleFractions_[k] /= sum;
+      }
+
+      break;
+
+    case  cIonSolnType_SINGLECATION:
+
+      throw CanteraError("eosType", "Unknown type");
+     
+      break;
+     
+    case  cIonSolnType_MULTICATIONANION:
+
+      throw CanteraError("eosType", "Unknown type");
+      break;
+
+    default:
+
+      throw CanteraError("eosType", "Unknown type");
+      break;
+
+    } 
+  } 
+
+// Calculate neutral molecule mole fractions
+  /*
+   *  This routine calculates the neutral molecule mole
+   *  fraction given the vector of ion mole fractions,
+   *  i.e., the mole fractions from this ThermoPhase.
+   *  Note, this routine basically assumes that there
+   *  is charge neutrality. If there isn't, then it wouldn't
+   *  make much sense. 
+   *
+   *  for the case of  cIonSolnType_SINGLEANION, some slough
+   *  in the charge neutrality is allowed. The cation number
+   *  is followed, while the difference in charge neutrality
+   *  is dumped into the anion mole number to fix the imbalance.
+   */
+  void IonsFromNeutralVPSSTP::getNeutralMoleculeMoleGrads(const doublereal * const dx, doublereal *dy) const {
+    int k, icat, jNeut;
+    doublereal sumCat; 
+    doublereal sumAnion;
+    doublereal fmij;
+    vector_fp y;
+    y.resize(numNeutralMoleculeSpecies_,0.0);
+    doublereal sumy, sumdy;
+
+    //! Zero the vector we are trying to find.
+    for (k = 0; k < numNeutralMoleculeSpecies_; k++) {
+      dy[k] = 0.0;
+    }
+
+
+    // bool fmSimple = true;
+
+    switch (ionSolnType_) {
+
+    case cIonSolnType_PASSTHROUGH:
+
+      for (k = 0; k < m_kk; k++) {
+	dy[k] = dx[k];
+      }
+      break;
+
+    case cIonSolnType_SINGLEANION:
+
+      sumCat = 0.0;
+      sumAnion = 0.0;
+
+      for (k = 0; k < (int) cationList_.size(); k++) {
+	//! Get the id for the next cation
+        icat = cationList_[k];
+	jNeut = fm_invert_ionForNeutral[icat];
+	if (jNeut >= 0) {
+	  fmij =  fm_neutralMolec_ions_[icat + jNeut * m_kk];
+	  AssertTrace(fmij != 0.0);
+	  dy[jNeut] += dx[icat] / fmij;
+	  y[jNeut] += moleFractions_[icat] / fmij;
+	}
+      }
+
+      for (k = 0; k <  numPassThroughSpecies_; k++) {
+	icat = passThroughList_[k];
+	jNeut = fm_invert_ionForNeutral[icat];
+	fmij = fm_neutralMolec_ions_[ icat + jNeut * m_kk];
+	dy[jNeut] += dx[icat] / fmij;
+	y[jNeut] += moleFractions_[icat] / fmij;
+      }
+
+#ifdef DEBUG_MODE
+      for (k = 0; k < m_kk; k++) {
+	moleFractionsTmp_[k] = dx[k];
+      }
+      for (jNeut = 0; jNeut <  numNeutralMoleculeSpecies_; jNeut++) {
+	for (k = 0; k < m_kk; k++) {
+	  fmij =  fm_neutralMolec_ions_[k + jNeut * m_kk];
+	  moleFractionsTmp_[k] -= fmij * dy[jNeut];
+	}
+      }
+      for (k = 0; k < m_kk; k++) {
+	if (fabs(moleFractionsTmp_[k]) > 1.0E-13) {
+	  //! Check to see if we have in fact found the inverse.
+	  if (anionList_[0] != k) {
+	    throw CanteraError("", "neutral molecule calc error");
+	  } else {
+	    //! For the single anion case, we will allow some slippage
+	    if (fabs(moleFractionsTmp_[k]) > 1.0E-5) {
+	      throw CanteraError("", "neutral molecule calc error - anion");
+	    }
+	  }
+	}
+      }
+#endif      
+
+      // Normalize the Neutral Molecule mole fractions
+      sumy = 0.0;
+      sumdy = 0.0;
+      for (k = 0; k < numNeutralMoleculeSpecies_; k++) {
+        sumy += y[k];
+	sumdy += dy[k];
+      }
+      for (k = 0; k < numNeutralMoleculeSpecies_; k++) {
+	dy[k] = dy[k]/sumy - y[k]*sumdy/sumy/sumy;
       }
 
       break;
@@ -815,6 +970,7 @@ namespace Cantera {
     } 
   }
 
+
   void IonsFromNeutralVPSSTP::setMassFractions(const doublereal* const y) {
     GibbsExcessVPSSTP::setMassFractions(y);
     calcNeutralMoleculeMoleFractions();
@@ -836,7 +992,7 @@ namespace Cantera {
   void IonsFromNeutralVPSSTP::setMoleFractions_NoNorm(const doublereal* const x) {
     GibbsExcessVPSSTP::setMoleFractions_NoNorm(x);
     calcNeutralMoleculeMoleFractions();
-    neutralMoleculePhase_->setMoleFractions(DATA_PTR(NeutralMolecMoleFractions_));
+    neutralMoleculePhase_->setMoleFractions_NoNorm(DATA_PTR(NeutralMolecMoleFractions_));
   }
 
 
@@ -1030,7 +1186,8 @@ namespace Cantera {
     muNeutralMolecule_.resize(numNeutralMoleculeSpecies_);
     gammaNeutralMolecule_.resize(numNeutralMoleculeSpecies_);
     dlnActCoeffdT_NeutralMolecule_.resize(numNeutralMoleculeSpecies_);
-    dlnActCoeffdlnC_NeutralMolecule_.resize(numNeutralMoleculeSpecies_);
+    dlnActCoeffdlnX_NeutralMolecule_.resize(numNeutralMoleculeSpecies_);
+    dlnActCoeffdlnN_NeutralMolecule_.resize(numNeutralMoleculeSpecies_);
   }
 
   static double factorOverlap(const std::vector<std::string>&  elnamesVN ,
@@ -1243,7 +1400,7 @@ namespace Cantera {
         icat = cationList_[k];
 	jNeut = fm_invert_ionForNeutral[icat];
 	fmij =  fm_neutralMolec_ions_[icat + jNeut * m_kk];
-        lnActCoeff_Scaled_[icat] = fmij * log(gammaNeutralMolecule_[jNeut]);
+        lnActCoeff_Scaled_[icat] = log(gammaNeutralMolecule_[jNeut])/fmij;
       }
 
       // Do the anion list
@@ -1273,6 +1430,75 @@ namespace Cantera {
   }
 
 
+
+  // get the gradient in the activity coefficients
+
+  void IonsFromNeutralVPSSTP::getdlnActCoeff(const doublereal dT, const doublereal * const dX, doublereal *dlnActCoeff) const {
+    int k, icat, jNeut;
+    doublereal fmij;
+    int numNeutMolSpec;
+    /*
+     * Get the activity coefficients of the neutral molecules
+     */
+    GibbsExcessVPSSTP *geThermo = dynamic_cast<GibbsExcessVPSSTP *>(neutralMoleculePhase_);
+    if (!geThermo) {
+      for ( k = 0; k < m_kk; k++ ){
+	dlnActCoeff[k] = dX[k]/moleFractions_[k];
+      }
+      return;
+    }
+
+    numNeutMolSpec = geThermo->nSpecies();
+    vector_fp dlnActCoeff_NeutralMolecule(numNeutMolSpec);
+    vector_fp dX_NeutralMolecule(numNeutMolSpec);
+
+
+    getNeutralMoleculeMoleGrads(DATA_PTR(dX),DATA_PTR(dX_NeutralMolecule));
+
+    // All mole fractions returned to normal
+ 
+    geThermo->getdlnActCoeff(dT, DATA_PTR(dX_NeutralMolecule), DATA_PTR(dlnActCoeff_NeutralMolecule));
+
+    switch (ionSolnType_) {
+    case cIonSolnType_PASSTHROUGH:
+      break;
+    case cIonSolnType_SINGLEANION:
+   
+      // Do the cation list
+      for (k = 0; k < (int) cationList_.size(); k++) {
+	//! Get the id for the next cation
+        icat = cationList_[k];
+	jNeut = fm_invert_ionForNeutral[icat];
+	fmij =  fm_neutralMolec_ions_[icat + jNeut * m_kk];
+        dlnActCoeff[icat] = dlnActCoeff_NeutralMolecule[jNeut]/fmij;
+      }
+
+      // Do the anion list
+      icat = anionList_[0];
+      jNeut = fm_invert_ionForNeutral[icat];
+      dlnActCoeff[icat]= 0.0;
+
+      // Do the list of neutral molecules
+      for (k = 0; k <  numPassThroughSpecies_; k++) {
+	icat = passThroughList_[k];
+	jNeut = fm_invert_ionForNeutral[icat];
+	dlnActCoeff[icat] = dlnActCoeff_NeutralMolecule[jNeut];
+      }
+      break;
+ 
+    case cIonSolnType_SINGLECATION:
+      throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff", "Unimplemented type");
+      break;
+    case cIonSolnType_MULTICATIONANION:
+      throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff", "Unimplemented type");
+      break;
+    default:
+      throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff", "Unimplemented type");
+      break;
+    }
+
+  }
+
   // Update the temperatture derivative of the ln activity coefficients
   /*
    * This function will be called to update the internally storred
@@ -1303,7 +1529,7 @@ namespace Cantera {
         icat = cationList_[k];
 	jNeut = fm_invert_ionForNeutral[icat];
 	fmij =  fm_neutralMolec_ions_[icat + jNeut * m_kk];
-        dlnActCoeffdT_Scaled_[icat] = fmij * dlnActCoeffdT_NeutralMolecule_[jNeut];
+        dlnActCoeffdT_Scaled_[icat] = dlnActCoeffdT_NeutralMolecule_[jNeut]/fmij;
       }
 
       // Do the anion list
@@ -1336,7 +1562,7 @@ namespace Cantera {
    * This function will be called to update the internally storred
    * temperature derivative of the natural logarithm of the activity coefficients
    */
-  void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnC() const {
+  void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnX() const {
     int k, icat, jNeut;
     doublereal fmij;
     /*
@@ -1344,11 +1570,11 @@ namespace Cantera {
      */
     GibbsExcessVPSSTP *geThermo = dynamic_cast<GibbsExcessVPSSTP *>(neutralMoleculePhase_);
     if (!geThermo) {
-      fvo_zero_dbl_1(dlnActCoeffdlnC_Scaled_, m_kk);
+      fvo_zero_dbl_1(dlnActCoeffdlnX_Scaled_, m_kk);
       return;
     }
 
-    geThermo->getdlnActCoeffdlnC(DATA_PTR(dlnActCoeffdlnC_NeutralMolecule_));
+    geThermo->getdlnActCoeffdlnX(DATA_PTR(dlnActCoeffdlnX_NeutralMolecule_));
 
     switch (ionSolnType_) {
     case cIonSolnType_PASSTHROUGH:
@@ -1361,19 +1587,77 @@ namespace Cantera {
         icat = cationList_[k];
 	jNeut = fm_invert_ionForNeutral[icat];
 	fmij =  fm_neutralMolec_ions_[icat + jNeut * m_kk];
-        dlnActCoeffdlnC_Scaled_[icat] = fmij * dlnActCoeffdlnC_NeutralMolecule_[jNeut];
+        dlnActCoeffdlnX_Scaled_[icat] = dlnActCoeffdlnX_NeutralMolecule_[jNeut]/fmij;
       }
 
       // Do the anion list
       icat = anionList_[0];
       jNeut = fm_invert_ionForNeutral[icat];
-      dlnActCoeffdT_Scaled_[icat]= 0.0;
+      dlnActCoeffdlnX_Scaled_[icat]= 0.0;
 
       // Do the list of neutral molecules
       for (k = 0; k <  numPassThroughSpecies_; k++) {
 	icat = passThroughList_[k];
 	jNeut = fm_invert_ionForNeutral[icat];
-	dlnActCoeffdlnC_Scaled_[icat] = dlnActCoeffdlnC_NeutralMolecule_[jNeut];
+	dlnActCoeffdlnX_Scaled_[icat] = dlnActCoeffdlnX_NeutralMolecule_[jNeut];
+      }
+      break;
+ 
+    case cIonSolnType_SINGLECATION:
+      throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff", "Unimplemented type");
+      break;
+    case cIonSolnType_MULTICATIONANION:
+      throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff", "Unimplemented type");
+      break;
+    default:
+      throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff", "Unimplemented type");
+      break;
+    }
+
+  }
+
+  /*
+   * This function will be called to update the internally storred
+   * temperature derivative of the natural logarithm of the activity coefficients
+   */
+  void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnN() const {
+    int k, icat, jNeut;
+    doublereal fmij;
+    /*
+     * Get the activity coefficients of the neutral molecules
+     */
+    GibbsExcessVPSSTP *geThermo = dynamic_cast<GibbsExcessVPSSTP *>(neutralMoleculePhase_);
+    if (!geThermo) {
+      fvo_zero_dbl_1(dlnActCoeffdlnN_Scaled_, m_kk);
+      return;
+    }
+
+    geThermo->getdlnActCoeffdlnN(DATA_PTR(dlnActCoeffdlnN_NeutralMolecule_));
+
+    switch (ionSolnType_) {
+    case cIonSolnType_PASSTHROUGH:
+      break;
+    case cIonSolnType_SINGLEANION:
+   
+      // Do the cation list
+      for (k = 0; k < (int) cationList_.size(); k++) {
+	//! Get the id for the next cation
+        icat = cationList_[k];
+	jNeut = fm_invert_ionForNeutral[icat];
+	fmij =  fm_neutralMolec_ions_[icat + jNeut * m_kk];
+        dlnActCoeffdlnN_Scaled_[icat] = dlnActCoeffdlnN_NeutralMolecule_[jNeut]/fmij;
+      }
+
+      // Do the anion list
+      icat = anionList_[0];
+      jNeut = fm_invert_ionForNeutral[icat];
+      dlnActCoeffdlnN_Scaled_[icat]= 0.0;
+
+      // Do the list of neutral molecules
+      for (k = 0; k <  numPassThroughSpecies_; k++) {
+	icat = passThroughList_[k];
+	jNeut = fm_invert_ionForNeutral[icat];
+	dlnActCoeffdlnN_Scaled_[icat] = dlnActCoeffdlnN_NeutralMolecule_[jNeut];
       }
       break;
  
diff --git a/Cantera/src/thermo/IonsFromNeutralVPSSTP.h b/Cantera/src/thermo/IonsFromNeutralVPSSTP.h
index 88b05e786..712ae622a 100644
--- a/Cantera/src/thermo/IonsFromNeutralVPSSTP.h
+++ b/Cantera/src/thermo/IonsFromNeutralVPSSTP.h
@@ -404,6 +404,21 @@ namespace Cantera {
      */
     virtual void getPartialMolarEntropies(doublereal* sbar) const;
 
+    //! Get the array of change in the log activity coefficients w.r.t. change in state (change temp, change mole fractions)
+    /*!
+     * This function is a virtual class, but it first appears in GibbsExcessVPSSTP
+     * class and derived classes from GibbsExcessVPSSTP.
+     *
+     * This function is a virtual method.  For ideal mixtures 
+     * (unity activity coefficients), this can gradX/X.  
+     *
+     * @param dT    Input of temperature change
+     * @param dX    Input vector of changes in mole fraction. length = m_kk
+     * @param dlnActCoeff    Output vector of derivatives of the 
+     *                         log Activity Coefficients. length = m_kk
+     */
+     virtual void getdlnActCoeff(const doublereal dT, const doublereal * const dX, doublereal *dlnActCoeff) const;
+
     //! Get the array of log concentration-like derivatives of the 
     //! log activity coefficients
     /*!
@@ -411,19 +426,65 @@ namespace Cantera {
      * (unity activity coefficients), this can return zero.  
      * Implementations should take the derivative of the 
      * logarithm of the activity coefficient with respect to the 
-     * logarithm of the concentration-like variable (i.e. mole fraction,
-     * molality, etc.) that represents the standard state.  
+     * logarithm of the concentration-like variable (i.e. mole fraction) 
+     * that represents the standard state.  
      * This quantity is to be used in conjunction with derivatives of 
      * that concentration-like variable when the derivative of the chemical 
      * potential is taken.  
      *
      *  units = dimensionless
      *
-     * @param dlnActCoeffdlnC    Output vector of log(mole fraction)  
+     * @param dlnActCoeffdlnX    Output vector of log(mole fraction)  
      *                 derivatives of the log Activity Coefficients.
      *                 length = m_kk
      */
-    virtual void getdlnActCoeffdlnC(doublereal *dlnActCoeffdlnC) const;
+    virtual void getdlnActCoeffdlnX(doublereal *dlnActCoeffdlnX) const;
+
+    //! Get the array of log concentration-like derivatives of the 
+    //! log activity coefficients
+    /*!
+     * This function is a virtual method.  For ideal mixtures 
+     * (unity activity coefficients), this can return zero.  
+     * Implementations should take the derivative of the 
+     * logarithm of the activity coefficient with respect to the 
+     * logarithm of the concentration-like variable (i.e. number of moles)
+     * that represents the standard state.  
+     * This quantity is to be used in conjunction with derivatives of 
+     * that concentration-like variable when the derivative of the chemical 
+     * potential is taken.  
+     *
+     *  units = dimensionless
+     *
+     * @param dlnActCoeffdlnN    Output vector of log(mole fraction)  
+     *                 derivatives of the log Activity Coefficients.
+     *                 length = m_kk
+     */
+    virtual void getdlnActCoeffdlnN(doublereal *dlnActCoeffdlnN) const;
+
+
+    virtual void getDissociationCoeffs(vector_fp& coeffs, vector_fp& charges);
+
+    virtual void getNeutralMolecMoleFractions(vector_fp& fracs){fracs=NeutralMolecMoleFractions_;}
+
+    //! Calculate neutral molecule mole fractions
+    /*!
+     *  This routine calculates the neutral molecule mole
+     *  fraction given the vector of ion mole fractions,
+     *  i.e., the mole fractions from this ThermoPhase.
+     *  Note, this routine basically assumes that there
+     *  is charge neutrality. If there isn't, then it wouldn't
+     *  make much sense. 
+     *
+     *  for the case of  cIonSolnType_SINGLEANION, some slough
+     *  in the charge neutrality is allowed. The cation number
+     *  is followed, while the difference in charge neutrality
+     *  is dumped into the anion mole number to fix the imbalance.
+     */
+    virtual void getNeutralMoleculeMoleGrads(const doublereal * const x, doublereal *y) const;
+
+    virtual void getCationList(std::vector<int>& cation){cation=cationList_;}
+    virtual void getAnionList(std::vector<int>& anion){anion=anionList_;}
+    virtual void getSpeciesNames(std::vector<std::string>& names){names=m_speciesNames;}
 
 
     //@}
@@ -680,6 +741,14 @@ namespace Cantera {
      */
     void s_update_dlnActCoeffdT() const;
 
+    //! Update the change in the ln activity coefficients
+    /*!
+     * This function will be called to update the internally storred
+     * change of the natural logarithm of the activity coefficients
+     * w.r.t a change in state (temp, mole fraction, etc)
+     */
+    void s_update_dlnActCoeff() const;
+
     //! Update the derivative of the log of the activity coefficients
     //!  wrt log(mole fraction)
     /*!
@@ -687,7 +756,16 @@ namespace Cantera {
      * derivative of the natural logarithm of the activity coefficients
      * wrt logarithm of the mole fractions.
      */
-    void s_update_dlnActCoeff_dlnC() const;
+    void s_update_dlnActCoeff_dlnX() const;
+
+    //! Update the derivative of the log of the activity coefficients
+    //!  wrt log(number of moles)
+    /*!
+     * This function will be called to update the internally storred
+     * derivative of the natural logarithm of the activity coefficients
+     * wrt logarithm of the number of moles of given species.
+     */
+    void s_update_dlnActCoeff_dlnN() const;
 
 
   private:
@@ -818,8 +896,10 @@ namespace Cantera {
 
     mutable std::vector<doublereal> muNeutralMolecule_;
     mutable std::vector<doublereal> gammaNeutralMolecule_;
+    mutable std::vector<doublereal> dlnActCoeff_NeutralMolecule_;
     mutable std::vector<doublereal> dlnActCoeffdT_NeutralMolecule_;
-    mutable std::vector<doublereal> dlnActCoeffdlnC_NeutralMolecule_;
+    mutable std::vector<doublereal> dlnActCoeffdlnX_NeutralMolecule_;
+    mutable std::vector<doublereal> dlnActCoeffdlnN_NeutralMolecule_;
 
   };
 
diff --git a/Cantera/src/thermo/MargulesVPSSTP.cpp b/Cantera/src/thermo/MargulesVPSSTP.cpp
index 2017d8357..a6340b7ae 100644
--- a/Cantera/src/thermo/MargulesVPSSTP.cpp
+++ b/Cantera/src/thermo/MargulesVPSSTP.cpp
@@ -99,6 +99,12 @@ namespace Cantera {
     m_SE_b_ij                   = b.m_SE_b_ij;
     m_SE_c_ij                   = b.m_SE_c_ij;
     m_SE_d_ij                   = b.m_SE_d_ij;
+    m_VHE_b_ij                  = b.m_VHE_b_ij;
+    m_VHE_c_ij                  = b.m_VHE_c_ij;
+    m_VHE_d_ij                  = b.m_VHE_d_ij;
+    m_VSE_b_ij                  = b.m_VSE_b_ij;
+    m_VSE_c_ij                  = b.m_VSE_c_ij;
+    m_VSE_d_ij                  = b.m_VSE_d_ij;
     m_pSpecies_A_ij             = b.m_pSpecies_A_ij;
     m_pSpecies_B_ij             = b.m_pSpecies_B_ij;
     formMargules_               = b.formMargules_;
@@ -154,11 +160,20 @@ namespace Cantera {
     m_SE_b_ij.resize(1);
     m_SE_c_ij.resize(1);
     m_SE_d_ij.resize(1);
+
+    m_VHE_b_ij.resize(1);
+    m_VHE_c_ij.resize(1);
+    m_VHE_d_ij.resize(1);
+
+    m_VSE_b_ij.resize(1);
+    m_VSE_c_ij.resize(1);
+    m_VSE_d_ij.resize(1);
     
     m_pSpecies_A_ij.resize(1);
     m_pSpecies_B_ij.resize(1);
 
 
+    
     m_HE_b_ij[0] = -17570E3;
     m_HE_c_ij[0] = -377.0E3;
     m_HE_d_ij[0] = 0.0;
@@ -166,6 +181,7 @@ namespace Cantera {
     m_SE_b_ij[0] = -7.627E3;
     m_SE_c_ij[0] =  4.958E3;
     m_SE_d_ij[0] =  0.0;
+    
 
     int iLiCl = speciesIndex("LiCl(L)");
     if (iLiCl < 0) {
@@ -217,7 +233,7 @@ namespace Cantera {
    */
   void MargulesVPSSTP::constructPhaseFile(std::string inputFile, std::string id) {
 
-    if (inputFile.size() == 0) {
+    if ((int) inputFile.size() == 0) {
       throw CanteraError("MargulesVPSSTP:constructPhaseFile",
                          "input file is null");
     }
@@ -274,7 +290,7 @@ namespace Cantera {
    */
   void MargulesVPSSTP::constructPhaseXML(XML_Node& phaseNode, std::string id) {
     string stemp;
-    if (id.size() > 0) {
+    if ((int) id.size() > 0) {
       string idp = phaseNode.id();
       if (idp != id) {
 	throw CanteraError("MargulesVPSSTP::constructPhaseXML", 
@@ -350,7 +366,7 @@ namespace Cantera {
      * take the exp of the internally storred coefficients.
      */
     for (int k = 0; k < m_kk; k++) {
-      ac[k] = exp(lnActCoeff_Scaled_[k]);      
+      ac[k] = exp(lnActCoeff_Scaled_[k]); 
     }
   }
 
@@ -472,7 +488,56 @@ namespace Cantera {
     }
   }
   
-  
+  /*
+   * ------------ Partial Molar Properties of the Solution ------------
+   */
+
+  // Return an array of partial molar volumes for the
+  // species in the mixture. Units: m^3/kmol.
+  /*
+   *  Frequently, for this class of thermodynamics representations,
+   *  the excess Volume due to mixing is zero. Here, we set it as
+   *  a default. It may be overriden in derived classes.
+   *
+   *  @param vbar   Output vector of speciar partial molar volumes.
+   *                Length = m_kk. units are m^3/kmol.
+   */
+  void MargulesVPSSTP::getPartialMolarVolumes(doublereal* vbar) const {
+
+    int iA, iB, iK, delAK, delBK;
+    double XA, XB, XK, g0 , g1;
+    double T = temperature();
+
+    /*
+     * Get the standard state values in m^3 kmol-1
+     */
+    getStandardVolumes(vbar);
+    //cout << "species name(0) = " << speciesName(0) << endl;
+    //cout << "iA = " << speciesName(m_pSpecies_A_ij[0]) << endl;   
+    //cout << "iB = " << speciesName(m_pSpecies_B_ij[0]) << endl;   
+    
+    for ( iK = 0; iK < m_kk; iK++ ){
+      delAK = 0;
+      delBK = 0;
+      XK = moleFractions_[iK];  
+      for (int i = 0; i <  numBinaryInteractions_; i++) {
+    
+	iA =  m_pSpecies_A_ij[i];    
+	iB =  m_pSpecies_B_ij[i];
+
+	if (iA==iK) delAK = 1;
+	else if (iB==iK) delBK = 1;
+	
+	XA = moleFractions_[iA];
+	XB = moleFractions_[iB];
+	
+	g0 = (m_VHE_b_ij[i] - T * m_VSE_b_ij[i]);
+	g1 = (m_VHE_c_ij[i] - T * m_VSE_c_ij[i]);
+	
+	vbar[iK] += XA*XB*(g0+g1*XB)+((delAK-XA)*XB+XA*(delBK-XB))*(g0+g1*XB)+XA*XB*(delBK-XB)*g1;
+      }
+    }
+  }  
 
   doublereal MargulesVPSSTP::err(std::string msg) const {
     throw CanteraError("MargulesVPSSTP","Base class method "
@@ -583,14 +648,54 @@ namespace Cantera {
 
   }
 
+ 
   // Update the activity coefficients
   /*
    * This function will be called to update the internally storred
    * natural logarithm of the activity coefficients
    *
-   *   he = X_A X_B(B + C(X_A - X_B))
+   *   he = X_A X_B(B + C X_B)
    */
+    void MargulesVPSSTP::s_update_lnActCoeff() const {
+    int iA, iB, iK, delAK, delBK;
+    double XA, XB, XK, g0 , g1;
+    double T = temperature();
+    double RT = GasConstant*T;
+
+    fvo_zero_dbl_1(lnActCoeff_Scaled_, m_kk);
+    
+    for ( iK = 0; iK < m_kk; iK++ ){
+
+      XK = moleFractions_[iK];  
+
+      for (int i = 0; i <  numBinaryInteractions_; i++) {
+    
+	iA =  m_pSpecies_A_ij[i];    
+	iB =  m_pSpecies_B_ij[i];
+
+	delAK = 0;
+	delBK = 0;
+
+	if (iA==iK) delAK = 1;
+	else if (iB==iK) delBK = 1;
+	
+	XA = moleFractions_[iA];
+	XB = moleFractions_[iB];
+	
+	g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
+	g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT;
+	
+	lnActCoeff_Scaled_[iK] += (delAK*XB+XA*delBK-XA*XB)*(g0+g1*XB)+XA*XB*(delBK-XB)*g1;
+	//lnActCoeff_Scaled_[iK] += XA*XB*(g0+g1*XB)+((delAK-XA)*XB+XA*(delBK-XB))*(g0+g1*XB)+XA*XB*(delBK-XB)*g1;
+      }
+    }
+  }
+
+
+  /*
+  // Not Right???
   void MargulesVPSSTP::s_update_lnActCoeff() const {
+
     int iA, iB;
     double XA, XB, g0 , g1;
     double T = temperature();
@@ -612,15 +717,52 @@ namespace Cantera {
       lnActCoeff_Scaled_[iB] += XA * XA * g0 +  XA * XB * g1 * (2 * XA);
     }
   }
+  */
 
   // Update the derivative of the log of the activity coefficients wrt T
   /*
    * This function will be called to update the internally storred
    * natural logarithm of the activity coefficients
    *
-   *   he = X_A X_B(B + C(X_A - X_B))
+   *   he = X_A X_B(B + C X_B)
    */
   void MargulesVPSSTP::s_update_dlnActCoeff_dT() const {
+    int iA, iB, iK, delAK, delBK;
+    double XA, XB, XK, g0 , g1;
+    double T = temperature();
+    double RTT = GasConstant*T*T;
+
+    fvo_zero_dbl_1(dlnActCoeffdT_Scaled_, m_kk);
+    
+    for ( iK = 0; iK < m_kk; iK++ ){
+
+      XK = moleFractions_[iK];  
+
+      for (int i = 0; i <  numBinaryInteractions_; i++) {
+    
+	iA =  m_pSpecies_A_ij[i];    
+	iB =  m_pSpecies_B_ij[i];
+	
+	delAK = 0;
+	delBK = 0;
+
+	if (iA==iK) delAK = 1;
+	else if (iB==iK) delBK = 1;
+	
+	XA = moleFractions_[iA];
+	XB = moleFractions_[iB];
+	
+	g0 = -m_HE_b_ij[i] / RTT;
+	g1 = -m_HE_c_ij[i] / RTT;
+	
+	dlnActCoeffdT_Scaled_[iK] += (delAK*XB+XA*delBK-XA*XB)*(g0+g1*XB)+XA*XB*(delBK-XB)*g1;
+      }
+    }
+  }
+
+  /* Not Right???
+  void MargulesVPSSTP::s_update_dlnActCoeff_dT() const {}
+
     int iA, iB;
     doublereal XA, XB, h0 , h1;
     doublereal T = temperature();
@@ -642,6 +784,7 @@ namespace Cantera {
       dlnActCoeffdT_Scaled_[iB] += -(XA * XA * h0 +  XA * XB * h1 * (2 * XA))/RTT;
     }
   }
+  */
 
   void MargulesVPSSTP::getdlnActCoeffdT(doublereal *dlnActCoeffdT) const {
     s_update_dlnActCoeff_dT();
@@ -650,22 +793,123 @@ namespace Cantera {
     }
   }
 
+  // calculate the change of the log of the activity coefficients wrt change in state: dT, dX
+  /*
+   * This function will be called to calculate gradient of the 
+   * logarithm of the activity coefficients based on gradients in temperature and mole fraction.
+   *
+   *   he = X_A X_B(B + C X_B)
+   */
+  void MargulesVPSSTP::getdlnActCoeff(const doublereal dT, const doublereal * const dX, doublereal* dlnActCoeff) const {
+    int iA, iB, iK, delAK, delBK;
+    double XA, XB, XK, g0 , g1, dXA, dXB;
+    double T = temperature();
+    double RT = GasConstant*T;
+
+    //fvo_zero_dbl_1(dlnActCoeff, m_kk);
+    s_update_dlnActCoeff_dT();
+
+    for ( iK = 0; iK < m_kk; iK++ ){
+
+      XK = moleFractions_[iK];  
+      dlnActCoeff[iK] = 0.0;
+
+      for (int i = 0; i <  numBinaryInteractions_; i++) {
+    
+	iA =  m_pSpecies_A_ij[i];    
+	iB =  m_pSpecies_B_ij[i];
+
+	delAK = 0;
+	delBK = 0;
+
+	if (iA==iK) delAK = 1;
+	else if (iB==iK) delBK = 1;
+	
+	XA = moleFractions_[iA];
+	XB = moleFractions_[iB];
+
+	dXA = dX[iA];
+	dXB = dX[iB];
+	
+	g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
+	g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT;
+	
+	dlnActCoeff[iK] += ((delBK-XB)*dXA + (delAK-XA)*dXB)*(g0+2*g1*XB) + (delBK-XB)*2*g1*XA*dXB + dlnActCoeffdT_Scaled_[iK]*dT;
+      }
+    }
+  }
+
   // Update the derivative of the log of the activity coefficients wrt ln(X)
   /*
    * This function will be called to update the internally stored gradients of the 
    * logarithm of the activity coefficients.  These are used in the determination 
    * of the diffusion coefficients.
    *
-   *   he = X_A X_B(B + C(X_A - X_B))
+   *   he = X_A X_B(B + C X_B)
    */
-  void MargulesVPSSTP::s_update_dlnActCoeff_dlnC() const {
+  void MargulesVPSSTP::s_update_dlnActCoeff_dlnN() const {
+    int iA, iB, iK, delAK, delBK;
+    double XA, XB, XK, g0 , g1;
+    double T = temperature();
+    double RT = GasConstant*T;
+
+    fvo_zero_dbl_1(dlnActCoeffdlnN_Scaled_, m_kk);
+    
+    for ( iK = 0; iK < m_kk; iK++ ){
+
+      XK = moleFractions_[iK];  
+
+      for (int i = 0; i <  numBinaryInteractions_; i++) {
+    
+	iA =  m_pSpecies_A_ij[i];    
+	iB =  m_pSpecies_B_ij[i];
+
+	delAK = 0;
+	delBK = 0;
+
+	if (iA==iK) delAK = 1;
+	else if (iB==iK) delBK = 1;
+	
+	XA = moleFractions_[iA];
+	XB = moleFractions_[iB];
+	
+	g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
+	g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT;
+	
+	dlnActCoeffdlnN_Scaled_[iK] += 2*(delBK-XB)*(g0*(delAK-XA)+g1*(2*(delAK-XA)*XB+XA*(delBK-XB)));
+      }
+      dlnActCoeffdlnN_Scaled_[iK] = XK*dlnActCoeffdlnN_Scaled_[iK]-XK;
+    }
+  }
+
+  void MargulesVPSSTP::s_update_dlnActCoeff_dlnX() const {
+
     int iA, iB;
     doublereal XA, XB, g0 , g1;
     doublereal T = temperature();
 
-    fvo_zero_dbl_1(dlnActCoeffdlnC_Scaled_, m_kk);
+    fvo_zero_dbl_1(dlnActCoeffdlnX_Scaled_, m_kk);
 
     doublereal RT = GasConstant * T;
+    
+    
+    for (int i = 0; i <  numBinaryInteractions_; i++) {
+      
+      iA =  m_pSpecies_A_ij[i];    
+      iB =  m_pSpecies_B_ij[i];
+      
+      XA = moleFractions_[iA];
+      XB = moleFractions_[iB];
+      
+      g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
+      g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT;
+      
+      dlnActCoeffdlnX_Scaled_[iA] += XA*XB*(2*g1*-2*g0-6*g1*XB);
+      dlnActCoeffdlnX_Scaled_[iB] += XA*XB*(2*g1*-2*g0-6*g1*XB);
+    }
+    
+    /*
+    // Wrong!!!
     for (int i = 0; i <  numBinaryInteractions_; i++) {
       iA =  m_pSpecies_A_ij[i];    
       iB =  m_pSpecies_B_ij[i];
@@ -676,17 +920,25 @@ namespace Cantera {
       g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT ;
       g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT;
       
-      dlnActCoeffdlnC_Scaled_[iA] += XA * ( ( - 2.0 + 2.0 * XA ) * g0
+      dlnActCoeffdlnX_Scaled_[iA] += XA * ( ( - 2.0 + 2.0 * XA ) * g0
 					    + ( - 4.0 + 10.0 * XA - 6.0 * XA*XA ) * g1 ) ;
-      dlnActCoeffdlnC_Scaled_[iB] += XB * ( ( - 2.0 + 2.0 * XB ) * g0
+      dlnActCoeffdlnX_Scaled_[iB] += XB * ( ( - 2.0 + 2.0 * XB ) * g0
 					    + (   2.0 -  8.0 * XB + 6.0 * XB*XB ) * g1 ) ;
     }
+    */
   }
+  
 
-  void MargulesVPSSTP::getdlnActCoeffdlnC(doublereal *dlnActCoeffdlnC) const {
-    s_update_dlnActCoeff_dlnC();
+  void MargulesVPSSTP::getdlnActCoeffdlnN(doublereal *dlnActCoeffdlnN) const {
+    s_update_dlnActCoeff_dlnN();
     for (int k = 0; k < m_kk; k++) {
-      dlnActCoeffdlnC[k] = dlnActCoeffdlnC_Scaled_[k];
+      dlnActCoeffdlnN[k] = dlnActCoeffdlnN_Scaled_[k];
+    }
+  }
+  void MargulesVPSSTP::getdlnActCoeffdlnX(doublereal *dlnActCoeffdlnX) const {
+    s_update_dlnActCoeff_dlnX();
+    for (int k = 0; k < m_kk; k++) {
+      dlnActCoeffdlnX[k] = dlnActCoeffdlnX_Scaled_[k];
     }
   }
 
@@ -699,6 +951,12 @@ namespace Cantera {
     m_SE_b_ij.resize(num, 0.0);
     m_SE_c_ij.resize(num, 0.0);
     m_SE_d_ij.resize(num, 0.0);
+    m_VHE_b_ij.resize(num, 0.0);
+    m_VHE_c_ij.resize(num, 0.0);
+    m_VHE_d_ij.resize(num, 0.0);
+    m_VSE_b_ij.resize(num, 0.0);
+    m_VSE_c_ij.resize(num, 0.0);
+    m_VSE_d_ij.resize(num, 0.0);
 
     m_pSpecies_A_ij.resize(num, -1);
     m_pSpecies_B_ij.resize(num, -1);
@@ -769,7 +1027,7 @@ namespace Cantera {
         /*
          * Get the string containing all of the values
          */
-        getFloatArray(xmlChild, vParams, true, "", "excessEnthalpy");
+        getFloatArray(xmlChild, vParams, true, "toSI", "excessEnthalpy");
         nParamsFound = vParams.size();
        
 	if (nParamsFound != 2) {
@@ -785,7 +1043,7 @@ namespace Cantera {
         /*
          * Get the string containing all of the values
          */
-        getFloatArray(xmlChild, vParams, true, "", "excessEntropy");
+        getFloatArray(xmlChild, vParams, true, "toSI", "excessEntropy");
         nParamsFound = vParams.size();
        
 	if (nParamsFound != 2) {
@@ -797,6 +1055,38 @@ namespace Cantera {
 	m_SE_c_ij[iSpot] = vParams[1];
       }
 
+      if (nodeName == "excessvolume_enthalpy") {
+        /*
+         * Get the string containing all of the values
+         */
+        getFloatArray(xmlChild, vParams, true, "toSI", "excessVolume_Enthalpy");
+        nParamsFound = vParams.size();
+       
+	if (nParamsFound != 2) {
+	  throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies::excessVolume_Enthalpy for " + ispName
+			     + "::" + jspName,
+			     "wrong number of params found");
+	}
+	m_VHE_b_ij[iSpot] = vParams[0];
+	m_VHE_c_ij[iSpot] = vParams[1];
+      }
+
+      if (nodeName == "excessvolume_entropy") {
+        /*
+         * Get the string containing all of the values
+         */
+        getFloatArray(xmlChild, vParams, true, "toSI", "excessVolume_Entropy");
+        nParamsFound = vParams.size();
+       
+	if (nParamsFound != 2) {
+	  throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies::excessVolume_Entropy for " + ispName
+			     + "::" + jspName,
+			     "wrong number of params found");
+	}
+	m_VSE_b_ij[iSpot] = vParams[0];
+	m_VSE_c_ij[iSpot] = vParams[1];
+      }
+
 
     } 
     
diff --git a/Cantera/src/thermo/MargulesVPSSTP.h b/Cantera/src/thermo/MargulesVPSSTP.h
index e9044b1bb..4d8a40776 100644
--- a/Cantera/src/thermo/MargulesVPSSTP.h
+++ b/Cantera/src/thermo/MargulesVPSSTP.h
@@ -565,6 +565,18 @@ namespace Cantera {
     virtual void getPartialMolarEntropies(doublereal* sbar) const;
 
     
+    //! Return an array of partial molar volumes for the
+    //! species in the mixture. Units: m^3/kmol.
+    /*!
+     *  Frequently, for this class of thermodynamics representations,
+     *  the excess Volume due to mixing is zero. Here, we set it as
+     *  a default. It may be overriden in derived classes.
+     *
+     *  @param vbar   Output vector of speciar partial molar volumes.
+     *                Length = m_kk. units are m^3/kmol.
+     */
+    virtual void getPartialMolarVolumes(doublereal* vbar) const;
+
     //! Get the species electrochemical potentials.
     /*!
      * These are partial molar quantities.
@@ -579,6 +591,19 @@ namespace Cantera {
     void getElectrochemPotentials(doublereal* mu) const;
 
  
+    //! Get the array of change in the log activity coefficients with change in state (change temp, change mole fractions)
+    /*!
+     * This function is a virtual class, but it first appears in GibbsExcessVPSSTP
+     * class and derived classes from GibbsExcessVPSSTP.
+     *
+     *  units = 1/Kelvin
+     *
+     * @param dlnActCoeff    Output vector of temperature derivatives of the 
+     *                         log Activity Coefficients. length = m_kk
+     *
+     */
+    virtual void getdlnActCoeff(const doublereal dT, const doublereal * const dX, doublereal *dlnActCoeffdT) const;
+
     //! Get the array of temperature derivatives of the log activity coefficients
     /*!
      * This function is a virtual class, but it first appears in GibbsExcessVPSSTP
@@ -608,11 +633,12 @@ namespace Cantera {
      *
      *  units = dimensionless
      *
-     * @param dlnActCoeffdlnC    Output vector of log(mole fraction)  
+     * @param dlnActCoeffdlnX    Output vector of log(mole fraction)  
      *                 derivatives of the log Activity Coefficients.
      *                 length = m_kk
      */
-    virtual void getdlnActCoeffdlnC(doublereal *dlnActCoeffdlnC) const;
+    virtual void getdlnActCoeffdlnX(doublereal *dlnActCoeffdlnX) const;
+    virtual void getdlnActCoeffdlnN(doublereal *dlnActCoeffdlnN) const;
 
  
     //@}
@@ -761,7 +787,16 @@ namespace Cantera {
      * derivative of the natural logarithm of the activity coefficients
      * wrt logarithm of the mole fractions.
      */
-    void s_update_dlnActCoeff_dlnC() const;
+    void s_update_dlnActCoeff_dlnX() const;
+
+    //! Update the derivative of the log of the activity coefficients
+    //!  wrt log(moles)
+    /*!
+     * This function will be called to update the internally storred
+     * derivative of the natural logarithm of the activity coefficients
+     * wrt logarithm of the moles.
+     */
+    void s_update_dlnActCoeff_dlnN() const;
 
 
   private:
@@ -802,6 +837,30 @@ namespace Cantera {
     //! Entropy term for the quaternary mole fraction interaction of the
     //! excess gibbs free energy expression
     mutable vector_fp m_SE_d_ij;
+
+    //! Enthalpy term for the binary mole fraction interaction of the
+    //! excess gibbs free energy expression
+    mutable vector_fp m_VHE_b_ij;
+
+    //! Enthalpy term for the ternary mole fraction interaction of the
+    //! excess gibbs free energy expression
+    mutable vector_fp m_VHE_c_ij;
+
+    //! Enthalpy term for the quaternary mole fraction interaction of the
+    //! excess gibbs free energy expression
+    mutable vector_fp m_VHE_d_ij;
+
+    //! Entropy term for the binary mole fraction interaction of the
+    //! excess gibbs free energy expression
+    mutable vector_fp m_VSE_b_ij;
+
+    //! Entropy term for the ternary mole fraction interaction of the
+    //! excess gibbs free energy expression
+    mutable vector_fp m_VSE_c_ij;
+
+    //! Entropy term for the quaternary mole fraction interaction of the
+    //! excess gibbs free energy expression
+    mutable vector_fp m_VSE_d_ij;
     
     //! vector of species indices representing species A in the interaction
     /*!
diff --git a/Cantera/src/thermo/PDSS_SSVol.cpp b/Cantera/src/thermo/PDSS_SSVol.cpp
index 71ae96fd0..4ca9cd46f 100644
--- a/Cantera/src/thermo/PDSS_SSVol.cpp
+++ b/Cantera/src/thermo/PDSS_SSVol.cpp
@@ -134,14 +134,14 @@ namespace Cantera {
       m_constMolarVolume = getFloat(*ss, "molarVolume", "toSI");
     } else if (model == "temperature_polynomial") {
       volumeModel_ = cSSVOLUME_TPOLY;
-      int num = getFloatArray(*ss, TCoeff_, true, "", "volumeTemperaturePolynomial");
+      int num = getFloatArray(*ss, TCoeff_, true, "toSI", "volumeTemperaturePolynomial");
       if (num != 4) {
 	throw CanteraError("PDSS_SSVol::constructPDSSXML",
 			   " Didn't get 4 density polynomial numbers for species " + speciesNode.name());
       }
     } else if (model == "density_temperature_polynomial") {
       volumeModel_ = cSSVOLUME_DENSITY_TPOLY;
-      int num = getFloatArray(*ss, TCoeff_, true, "", "densityTemperaturePolynomial");
+      int num = getFloatArray(*ss, TCoeff_, true, "toSI", "densityTemperaturePolynomial");
       if (num != 4) {
 	throw CanteraError("PDSS_SSVol::constructPDSSXML",
 			   " Didn't get 4 density polynomial numbers for species " + speciesNode.name());
diff --git a/Cantera/src/thermo/State.cpp b/Cantera/src/thermo/State.cpp
index 6191ce6bd..f8b3dbee9 100644
--- a/Cantera/src/thermo/State.cpp
+++ b/Cantera/src/thermo/State.cpp
@@ -196,6 +196,10 @@ namespace Cantera {
     return density()/meanMolecularWeight(); 
   }
 
+  doublereal State::molarVolume() const { 
+    return 1.0/molarDensity(); 
+  }
+
   void State::setConcentrations(const doublereal* const conc) {
     int k;
     doublereal sum = 0.0, norm = 0.0;
diff --git a/Cantera/src/thermo/State.h b/Cantera/src/thermo/State.h
index d05941add..ca2826df9 100644
--- a/Cantera/src/thermo/State.h
+++ b/Cantera/src/thermo/State.h
@@ -318,6 +318,9 @@ namespace Cantera {
     /// Molar density (kmol/m^3).
     doublereal molarDensity() const;
 
+    /// Molar density (kmol/m^3).
+    doublereal molarVolume() const;
+
     //! Set the internally storred density (kg/m^3) of the phase
     /*!
      * Note the density of a phase is an indepedent variable.
diff --git a/Cantera/src/thermo/ThermoPhase.h b/Cantera/src/thermo/ThermoPhase.h
index 3d5ddb577..2df0b8546 100644
--- a/Cantera/src/thermo/ThermoPhase.h
+++ b/Cantera/src/thermo/ThermoPhase.h
@@ -883,6 +883,21 @@ namespace Cantera {
     }
 
 
+    //! Get the change in activity coefficients w.r.t. change in state 
+    //! (temp, mole fraction, etc.)
+    /*!
+     * This function is a virtual method.  For ideal mixtures 
+     * (unity activity coefficients), this can gradX/X.  
+     *
+     * @param dT    Input of temperature change
+     * @param dX    Input vector of changes in mole fraction. length = m_kk
+     * @param dlnActCoeff    Output vector of derivatives of the 
+     *                         log Activity Coefficients. length = m_kk
+     */
+    virtual void getdlnActCoeff(const doublereal dT, const doublereal * const dX, doublereal *dlnActCoeff) const {
+      err("getdlnActCoeff");
+    }
+
     //! Get the array of log concentration-like derivatives of the 
     //! log activity coefficients
     /*!
@@ -890,19 +905,41 @@ namespace Cantera {
      * (unity activity coefficients), this can return zero.  
      * Implementations should take the derivative of the 
      * logarithm of the activity coefficient with respect to the 
-     * logarithm of the concentration-like variable (i.e. mole fraction,
-     * molality, etc.) that represents the standard state.  
+     * logarithm of the concentration-like variable (i.e. mole fraction)
+     * that represents the standard state.  
      * This quantity is to be used in conjunction with derivatives of 
      * that concentration-like variable when the derivative of the chemical 
      * potential is taken.  
      *
      *  units = dimensionless
      *
-     * @param dlnActCoeffdlnC    Output vector of derivatives of the 
+     * @param dlnActCoeffdlnX    Output vector of derivatives of the 
      *                         log Activity Coefficients. length = m_kk
      */
-    virtual void getdlnActCoeffdlnC(doublereal *dlnActCoeffdlnC) const {
-      err("getdlnActCoeffdlnC");
+    virtual void getdlnActCoeffdlnX(doublereal *dlnActCoeffdlnX) const {
+      err("getdlnActCoeffdlnX");
+    }
+
+    //! Get the array of log concentration-like derivatives of the 
+    //! log activity coefficients
+    /*!
+     * This function is a virtual method.  For ideal mixtures 
+     * (unity activity coefficients), this can return zero.  
+     * Implementations should take the derivative of the 
+     * logarithm of the activity coefficient with respect to the 
+     * logarithm of the concentration-like variable (i.e. moles)
+     * that represents the standard state.  
+     * This quantity is to be used in conjunction with derivatives of 
+     * that concentration-like variable when the derivative of the chemical 
+     * potential is taken.  
+     *
+     *  units = dimensionless
+     *
+     * @param dlnActCoeffdlnN    Output vector of derivatives of the 
+     *                         log Activity Coefficients. length = m_kk
+     */
+    virtual void getdlnActCoeffdlnN(doublereal *dlnActCoeffdlnN) const {
+      err("getdlnActCoeffdlnN");
     }
 
 
diff --git a/Cantera/src/thermo/VPSSMgr.cpp b/Cantera/src/thermo/VPSSMgr.cpp
index 64ec7e22c..1204bce2c 100644
--- a/Cantera/src/thermo/VPSSMgr.cpp
+++ b/Cantera/src/thermo/VPSSMgr.cpp
@@ -280,7 +280,8 @@ namespace Cantera {
 
   void 
   VPSSMgr::getStandardVolumes_ref(doublereal *vol) const{
-    err("getStandardVolumes_ref");
+    getStandardVolumes(vol);
+    //err("getStandardVolumes_ref");
   }
 
   /*****************************************************************/
@@ -359,6 +360,7 @@ namespace Cantera {
     m_sss_R.resize(m_kk, 0.0);
     m_Vss.resize(m_kk, 0.0);
 
+
     // Storage used by the PDSS objects to store their
     // answers.
     mPDSS_h0_RT.resize(m_kk, 0.0);
diff --git a/Cantera/src/thermo/VPStandardStateTP.h b/Cantera/src/thermo/VPStandardStateTP.h
index c13412fe4..7990b2cff 100644
--- a/Cantera/src/thermo/VPStandardStateTP.h
+++ b/Cantera/src/thermo/VPStandardStateTP.h
@@ -127,19 +127,41 @@ namespace Cantera {
      * (unity activity coefficients), this can return zero.  
      * Implementations should take the derivative of the 
      * logarithm of the activity coefficient with respect to the 
-     * logarithm of the concentration-like variable (i.e. mole fraction,
-     * molality, etc.) that represents the standard state.  
+     * logarithm of the concentration-like variable (i.e. moles)
+     *  that represents the standard state.  
      * This quantity is to be used in conjunction with derivatives of 
      * that concentration-like variable when the derivative of the chemical 
      * potential is taken.  
      *
      *  units = dimensionless
      *
-     * @param dlnActCoeffdlnC    Output vector of derivatives of the 
+     * @param dlnActCoeffdlnN    Output vector of derivatives of the 
      *                         log Activity Coefficients. length = m_kk
      */
-    virtual void getdlnActCoeffdlnC(doublereal *dlnActCoeffdlnC) const {
-      err("getdlnActCoeffdlnC");
+    virtual void getdlnActCoeffdlnN(doublereal *dlnActCoeffdlnN) const {
+      err("getdlnActCoeffdlnN");
+    }
+
+    //! Get the array of log concentration-like derivatives of the 
+    //! log activity coefficients
+    /*!
+     * This function is a virtual method.  For ideal mixtures 
+     * (unity activity coefficients), this can return zero.  
+     * Implementations should take the derivative of the 
+     * logarithm of the activity coefficient with respect to the 
+     * logarithm of the concentration-like variable (i.e. mole fraction)
+     * that represents the standard state.  
+     * This quantity is to be used in conjunction with derivatives of 
+     * that concentration-like variable when the derivative of the chemical 
+     * potential is taken.  
+     *
+     *  units = dimensionless
+     *
+     * @param dlnActCoeffdlnX    Output vector of derivatives of the 
+     *                         log Activity Coefficients. length = m_kk
+     */
+    virtual void getdlnActCoeffdlnX(doublereal *dlnActCoeffdlnX) const {
+      err("getdlnActCoeffdlnX");
     }