diff --git a/src/equil/vcs_phaseStability.cpp b/src/equil/vcs_phaseStability.cpp
index 4c7f8cb0e..3e6f1ccd7 100644
--- a/src/equil/vcs_phaseStability.cpp
+++ b/src/equil/vcs_phaseStability.cpp
@@ -398,7 +398,6 @@ int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop)
         printf("ERROR vcs_popPhaseRxnStepSizes called for a phase that exists!");
         exit(-1);
     }
-    char anote[256];
     if (m_debug_print_lvl >= 2) {
         plogf("  ---  vcs_popPhaseRxnStepSizes() called to pop phase %s %d into existence\n",
               Vphase->PhaseName.c_str(), iphasePop);
@@ -426,12 +425,6 @@ int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop)
         if (s != 0.0) {
             double s_old = s;
             s = vcs_Hessian_diag_adj(irxn, s_old);
-#ifdef DEBUG_MODE
-            if (s_old != s) {
-                sprintf(anote, "Normal calc: diag adjusted from %g "
-                        "to %g due to act coeff",  s_old, s);
-            }
-#endif
             m_deltaMolNumSpecies[kspec] = -m_deltaGRxn_new[irxn] / s;
         } else {
             // Ok, s is equal to zero. We can not apply a sophisticated theory
@@ -449,18 +442,8 @@ int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop)
                     double negChangeComp = - stoicC * m_deltaMolNumSpecies[kspec];
                     if (negChangeComp > m_molNumSpecies_old[j]) {
                         if (m_molNumSpecies_old[j] > 0.0) {
-#ifdef DEBUG_MODE
-                            sprintf(anote, "Delta damped from %g "
-                                    "to %g due to component %lu (%10s) going neg", m_deltaMolNumSpecies[kspec],
-                                    -m_molNumSpecies_old[j]/stoicC, j,  m_speciesName[j].c_str());
-#endif
                             m_deltaMolNumSpecies[kspec] = - 0.5 * m_molNumSpecies_old[j] / stoicC;
                         } else {
-#ifdef DEBUG_MODE
-                            sprintf(anote, "Delta damped from %g "
-                                    "to %g due to component %lu (%10s) zero", m_deltaMolNumSpecies[kspec],
-                                    -m_molNumSpecies_old[j]/stoicC, j,  m_speciesName[j].c_str());
-#endif
                             m_deltaMolNumSpecies[kspec] = 0.0;
                         }
                     }
@@ -469,11 +452,6 @@ int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop)
         }
         // Implement a damping term that limits m_deltaMolNumSpecies to the size of the mole number
         if (-m_deltaMolNumSpecies[kspec] > m_molNumSpecies_old[kspec]) {
-#ifdef DEBUG_MODE
-            sprintf(anote, "Delta damped from %g "
-                    "to %g due to %s going negative", m_deltaMolNumSpecies[kspec],
-                    -m_molNumSpecies_old[kspec],  m_speciesName[kspec].c_str());
-#endif
             m_deltaMolNumSpecies[kspec] = -m_molNumSpecies_old[kspec];
         }
 
diff --git a/src/thermo/HMWSoln.cpp b/src/thermo/HMWSoln.cpp
index 3d7be737c..60619e8cc 100644
--- a/src/thermo/HMWSoln.cpp
+++ b/src/thermo/HMWSoln.cpp
@@ -1848,12 +1848,6 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
         exit(EXIT_FAILURE);
     }
 
-#ifdef DEBUG_MODE
-    int printE = 0;
-    if (temperature() == 323.15) {
-        printE = 0;
-    }
-#endif
     std::string sni,  snj, snk;
 
     /*
@@ -2066,17 +2060,6 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
             n = m_kk*i + j;
             counterIJ = m_CounterIJ[n];
 
-#ifdef DEBUG_MODE
-            if (printE) {
-                if (counterIJ == 2) {
-                    printf("%s %s\n", speciesName(i).c_str(),
-                           speciesName(j).c_str());
-                    printf("beta0MX[%d] = %g\n", (int) counterIJ, beta0MX[counterIJ]);
-                    printf("beta1MX[%d] = %g\n", (int) counterIJ, beta1MX[counterIJ]);
-                    printf("beta2MX[%d] = %g\n", (int) counterIJ, beta2MX[counterIJ]);
-                }
-            }
-#endif
             /*
              * both species have a non-zero charge, and one is positive
              * and the other is negative
@@ -2143,15 +2126,6 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
             } else {
                 CMX[counterIJ] = 0.0;
             }
-#ifdef DEBUG_MODE
-            if (printE) {
-                if (counterIJ == 2) {
-                    printf("%s %s\n", speciesName(i).c_str(),
-                           speciesName(j).c_str());
-                    printf("CphiMX[%d] = %g\n", (int) counterIJ, CphiMX[counterIJ]);
-                }
-            }
-#endif
 #ifdef DEBUG_MODE
             if (m_debugCalc) {
                 sni = speciesName(i);
@@ -2220,11 +2194,6 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
     Aphi = A_Debye_TP() / 3.0;
     F = -Aphi * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
                  + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
-#ifdef DEBUG_MODE
-    if (printE) {
-        printf("Aphi = %20.13g\n", Aphi);
-    }
-#endif
 #ifdef DEBUG_MODE
     if (m_debugCalc) {
         printf(" initial value of F = %10.6f \n", F);
@@ -2821,11 +2790,6 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
     } else {
         osmotic_coef = 1.0;
     }
-#ifdef DEBUG_MODE
-    if (printE) {
-        printf("OsmCoef - 1 = %20.13g\n", osmotic_coef - 1.0);
-    }
-#endif
 #ifdef DEBUG_MODE
     if (m_debugCalc) {
         printf(" term1=%10.6f sum1=%10.6f sum2=%10.6f "