From 64287bbc7399dd7e08b5e4b6e3bff07cb7034d16 Mon Sep 17 00:00:00 2001
From: Harry Moffat <hkmoffa@sandia.gov>
Date: Thu, 23 Oct 2008 15:18:26 +0000
Subject: [PATCH] Made explicit the dependence of the calculation of vcs_dfe()
 on m_speciesStatus[]. If a phase is zeroed explicitly and the status flag is
 set, then the mole fraction term is dropped out of the free energy
 calculation. This will ensure that the deficient free energy of formation is
 evaluated for each species within these phases.

---
 Cantera/src/equil/vcs_solve_TP.cpp | 36 ++++++++++++++++++++----------
 1 file changed, 24 insertions(+), 12 deletions(-)

diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp
index 673a94b4d..ae6a4a8c5 100644
--- a/Cantera/src/equil/vcs_solve_TP.cpp
+++ b/Cantera/src/equil/vcs_solve_TP.cpp
@@ -4125,7 +4125,6 @@ namespace VCSnonideal {
     double *actCoeff_ptr;
     double *tlogMoles;
     vcs_VolPhase *Vphase;
-    VCS_SPECIES_THERMO *st_ptr;
 
     double *feSpecies;
     double * molNum;
@@ -4244,18 +4243,23 @@ namespace VCSnonideal {
 #ifdef DEBUG_MODE
 	if (molNum[kspec] != m_phasePhi[iphase]) {
 	  plogf("We have an inconsistency!\n");
-	  exit(-1);
+	  std::exit(-1);
 	}
 	if (m_chargeSpecies[kspec] != -1.0) {
 	  plogf("We have an unexpected situation!\n");
-	  exit(-1);
+	  std::exit(-1);
 	}
 #endif
-	feSpecies[kspec] =
-	  m_SSfeSpecies[kspec] + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
+	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];
+	  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];
@@ -4300,7 +4304,12 @@ namespace VCSnonideal {
 	      + 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];
+	      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];
@@ -4342,12 +4351,16 @@ namespace VCSnonideal {
 	      std::exit(-1);
 	    }
 #endif
-	    feSpecies[kspec] = 
-	      m_SSfeSpecies[kspec] +
-	      m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase]; ;
+	    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];
+	      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];
@@ -4361,7 +4374,6 @@ namespace VCSnonideal {
 		    + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase]; 
 		}
 	      } else {
-		st_ptr = m_speciesThermoList[kspec];
 		feSpecies[kspec] = m_SSfeSpecies[kspec] 
 		  + log(actCoeff_ptr[kspec] * molNum[kspec]) 
 		  - tlogMoles[m_phaseID[kspec]] - m_lnMnaughtSpecies[kspec]