diff --git a/src/equil/ChemEquil.cpp b/src/equil/ChemEquil.cpp
index 6599356c6..81b368388 100644
--- a/src/equil/ChemEquil.cpp
+++ b/src/equil/ChemEquil.cpp
@@ -548,13 +548,8 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
             }
             // update the T estimate
             t0 = t0 + dt;
-            if (t0 <= tminPhase || t0 >= tmaxPhase) {
-                printf("We shouldn't be here\n");
-                exit(EXIT_FAILURE);
-            }
-            if (t0 < 100.) {
-                printf("t0 - we are here %g\n", t0);
-                exit(EXIT_FAILURE);
+            if (t0 <= tminPhase || t0 >= tmaxPhase || t0 < 100.0) {
+                throw CanteraError("ChemEquil::equilibrate", "T out of bounds");
             }
             s.setTemperature(t0);
         }
diff --git a/src/equil/MultiPhaseEquil.cpp b/src/equil/MultiPhaseEquil.cpp
index e6145c515..282250186 100644
--- a/src/equil/MultiPhaseEquil.cpp
+++ b/src/equil/MultiPhaseEquil.cpp
@@ -726,8 +726,7 @@ void MultiPhaseEquil::reportCSV(const std::string& reportFile)
 
     FILE* FP = fopen(reportFile.c_str(), "w");
     if (!FP) {
-        printf("Failure to open file\n");
-        exit(EXIT_FAILURE);
+        throw CanteraError("MultiPhaseEquil::reportCSV", "Failure to open file");
     }
     double Temp = m_mix->temperature();
     double pres = m_mix->pressure();
diff --git a/src/equil/vcs_MultiPhaseEquil.cpp b/src/equil/vcs_MultiPhaseEquil.cpp
index 139753042..2b5dc84fd 100644
--- a/src/equil/vcs_MultiPhaseEquil.cpp
+++ b/src/equil/vcs_MultiPhaseEquil.cpp
@@ -608,8 +608,8 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile)
 
     FILE* FP = fopen(reportFile.c_str(), "w");
     if (!FP) {
-        plogf("Failure to open file\n");
-        exit(EXIT_FAILURE);
+        throw CanteraError("vcs_MultiPhaseEquil::reportCSV",
+                           "Failure to open file");
     }
     double Temp = m_mix->temperature();
     double pres = m_mix->pressure();
@@ -742,8 +742,7 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile)
             if (!vcs_doubleEqual(fe[istart+k], mu[k])) {
                 fprintf(FP,"ERROR: incompatibility!\n");
                 fclose(FP);
-                plogf("ERROR: incompatibility!\n");
-                exit(EXIT_FAILURE);
+                throw CanteraError("vcs_MultiPhaseEquil::reportCSV", "incompatibility!");
             }
         }
 #endif
@@ -874,8 +873,8 @@ int  vcs_Cantera_to_vprob(Cantera::MultiPhase* mphase,
             }
             VolPhase->m_eqnState = VCS_EOS_UNK_CANTERA;
             if (!VolPhase->usingCanteraCalls()) {
-                plogf("vcs functions asked for, but unimplemented\n");
-                exit(EXIT_FAILURE);
+                throw CanteraError("vcs_Cantera_to_vprob",
+                                   "vcs functions asked for, but unimplemented");
             }
             break;
         }
@@ -1028,8 +1027,8 @@ int  vcs_Cantera_to_vprob(Cantera::MultiPhase* mphase,
                 ts_ptr->SS0_Model = VCS_SS0_NOTHANDLED;
                 ts_ptr->SSStar_Model = VCS_SSSTAR_NOTHANDLED;
                 if (!(ts_ptr->UseCanteraCalls)) {
-                    plogf("Cantera calls not being used -> exiting\n");
-                    exit(EXIT_FAILURE);
+                    throw CanteraError("vcs_Cantera_to_vprob",
+                                       "Cantera calls not being used -> aborting");
                 }
             }
 
diff --git a/src/equil/vcs_VolPhase.cpp b/src/equil/vcs_VolPhase.cpp
index 692d716a5..6b6cf31ad 100644
--- a/src/equil/vcs_VolPhase.cpp
+++ b/src/equil/vcs_VolPhase.cpp
@@ -14,10 +14,13 @@
 
 #include "cantera/thermo/ThermoPhase.h"
 #include "cantera/thermo/mix_defs.h"
+#include "cantera/base/stringUtils.h"
 
 #include <sstream>
 #include <cstdio>
 
+using namespace Cantera;
+
 namespace VCSnonideal
 {
 
@@ -197,19 +200,16 @@ void vcs_VolPhase::resize(const size_t phaseNum, const size_t nspecies,
                           const size_t numElem, const char* const phaseName,
                           const double molesInert)
 {
-    if (DEBUG_MODE_ENABLED && nspecies <= 0) {
-        plogf("nspecies Error\n");
-        exit(EXIT_FAILURE);
-    }
+    AssertThrowMsg(nspecies > 0, "vcs_VolPhase::resize", "nspecies Error");
     setTotalMolesInert(molesInert);
     m_phi = 0.0;
     m_phiVarIndex = npos;
 
     if (phaseNum == VP_ID_) {
         if (strcmp(PhaseName.c_str(), phaseName)) {
-            plogf("Strings are different: %s %s :unknown situation\n",
-                  PhaseName.c_str(), phaseName);
-            exit(EXIT_FAILURE);
+            throw CanteraError("vcs_VolPhase::resize",
+                               "Strings are different: " + PhaseName + " " +
+                               phaseName + " :unknown situation");
         }
     } else {
         VP_ID_ = phaseNum;
@@ -421,14 +421,14 @@ void vcs_VolPhase::setMoleFractionsState(const double totalMoles,
         // There are other ways to set the mole fractions when VCS_STATECALC
         // is set to a normal settting.
         if (vcsStateStatus != VCS_STATECALC_TMP) {
-            printf("vcs_VolPhase::setMolesFractionsState: inappropriate usage\n");
-            exit(EXIT_FAILURE);
+            throw CanteraError("vcs_VolPhase::setMolesFractionsState",
+                               "inappropriate usage");
         }
         m_UpToDate = false;
         m_vcsStateStatus = VCS_STATECALC_TMP;
         if (m_existence == VCS_PHASE_EXIST_ZEROEDPHASE) {
-            printf("vcs_VolPhase::setMolesFractionsState: inappropriate usage\n");
-            exit(EXIT_FAILURE);
+            throw CanteraError("vcs_VolPhase::setMolesFractionsState",
+                               "inappropriate usage");
         }
         m_existence = VCS_PHASE_EXIST_YES;
     } else {
@@ -440,9 +440,9 @@ void vcs_VolPhase::setMoleFractionsState(const double totalMoles,
     v_totalMoles = totalMoles;
     if (m_totalMolesInert > 0.0) {
         if (m_totalMolesInert > v_totalMoles) {
-            printf("vcs_VolPhase::setMolesFractionsState: inerts greater than total: %g %g\n",
-                   v_totalMoles,  m_totalMolesInert);
-            exit(EXIT_FAILURE);
+            throw CanteraError("vcs_VolPhase::setMolesFractionsState",
+                 "inerts greater than total: " + fp2str(v_totalMoles) + " " +
+                 fp2str(m_totalMolesInert));
         }
         fractotal = 1.0 - m_totalMolesInert/v_totalMoles;
     }
@@ -452,8 +452,8 @@ void vcs_VolPhase::setMoleFractionsState(const double totalMoles,
         sum += moleFractions[k];
     }
     if (sum == 0.0) {
-        printf("vcs_VolPhase::setMolesFractionsState: inappropriate usage\n");
-        exit(EXIT_FAILURE);
+        throw CanteraError("vcs_VolPhase::setMolesFractionsState",
+                           "inappropriate usage");
     }
     if (sum  != fractotal) {
         for (size_t k = 0; k < m_numSpecies; k++) {
@@ -470,28 +470,22 @@ void vcs_VolPhase::setMolesFromVCS(const int stateCalc,
     v_totalMoles = m_totalMolesInert;
 
     if (molesSpeciesVCS == 0) {
-        if (DEBUG_MODE_ENABLED && m_owningSolverObject == 0) {
-            printf("vcs_VolPhase::setMolesFromVCS  shouldn't be here\n");
-            exit(EXIT_FAILURE);
-        }
+        AssertThrowMsg(m_owningSolverObject, "vcs_VolPhase::setMolesFromVCS",
+                       "shouldn't be here");
         if (stateCalc == VCS_STATECALC_OLD) {
             molesSpeciesVCS = VCS_DATA_PTR(m_owningSolverObject->m_molNumSpecies_old);
         } else if (stateCalc == VCS_STATECALC_NEW) {
             molesSpeciesVCS = VCS_DATA_PTR(m_owningSolverObject->m_molNumSpecies_new);
         } else if (DEBUG_MODE_ENABLED) {
-            printf("vcs_VolPhase::setMolesFromVCS shouldn't be here\n");
-            exit(EXIT_FAILURE);
-        }
+            throw CanteraError("vcs_VolPhase::setMolesFromVCS", "shouldn't be here");        }
     } else if (DEBUG_MODE_ENABLED && m_owningSolverObject) {
         if (stateCalc == VCS_STATECALC_OLD) {
             if (molesSpeciesVCS != VCS_DATA_PTR(m_owningSolverObject->m_molNumSpecies_old)) {
-                printf("vcs_VolPhase::setMolesFromVCS shouldn't be here\n");
-                exit(EXIT_FAILURE);
+                throw CanteraError("vcs_VolPhase::setMolesFromVCS", "shouldn't be here");
             }
         } else if (stateCalc == VCS_STATECALC_NEW) {
             if (molesSpeciesVCS != VCS_DATA_PTR(m_owningSolverObject->m_molNumSpecies_new)) {
-                printf("vcs_VolPhase::setMolesFromVCS shouldn't be here\n");
-                exit(EXIT_FAILURE);
+                throw CanteraError("vcs_VolPhase::setMolesFromVCS", "shouldn't be here");
             }
         }
     }
@@ -570,10 +564,9 @@ void vcs_VolPhase::setMolesFromVCSCheck(const int vcsStateStatus,
         if (vcs_doubleEqual(Tcheck, v_totalMoles)) {
             Tcheck = v_totalMoles;
         } else {
-            plogf("vcs_VolPhase::setMolesFromVCSCheck: "
-                  "We have a consistency problem: %21.16g %21.16g\n",
-                  Tcheck, v_totalMoles);
-            exit(EXIT_FAILURE);
+            throw CanteraError("vcs_VolPhase::setMolesFromVCSCheck",
+                  "We have a consistency problem: " + fp2str(Tcheck) + " " +
+                  fp2str(v_totalMoles));
         }
     }
 }
@@ -718,8 +711,7 @@ double vcs_VolPhase::_updateVolPM() const
             double volI = m_totalMolesInert * Cantera::GasConstant * Temp_ / Pres_;
             m_totalVol += volI;
         } else {
-            printf("unknown situation\n");
-            exit(EXIT_FAILURE);
+            throw CanteraError("vcs_VolPhase::_updateVolPM", "unknown situation");
         }
     }
     m_UpToDate_VolPM = true;
@@ -911,12 +903,10 @@ void vcs_VolPhase::setTotalMoles(const double totalMols)
     v_totalMoles = totalMols;
     if (m_totalMolesInert > 0.0) {
         m_existence = VCS_PHASE_EXIST_ALWAYS;
-        if (DEBUG_MODE_ENABLED && totalMols < m_totalMolesInert) {
-            printf(" vcs_VolPhase::setTotalMoles:: ERROR totalMoles "
-                   "less than inert moles: %g %g\n",
-                   totalMols, m_totalMolesInert);
-            exit(EXIT_FAILURE);
-        }
+        AssertThrowMsg(totalMols >= m_totalMolesInert,
+                       "vcs_VolPhase::setTotalMoles",
+                       "totalMoles less than inert moles: " +
+                       fp2str(totalMols) + " " + fp2str(m_totalMolesInert));
     } else {
         if (m_singleSpecies && (m_phiVarIndex == 0)) {
             m_existence =  VCS_PHASE_EXIST_ALWAYS;
@@ -1019,9 +1009,8 @@ void vcs_VolPhase::setExistence(const int existence)
     if (existence == VCS_PHASE_EXIST_NO || existence == VCS_PHASE_EXIST_ZEROEDPHASE) {
         if (v_totalMoles != 0.0) {
             if (DEBUG_MODE_ENABLED) {
-                plogf("vcs_VolPhase::setExistence setting false existence for phase with moles");
-                plogendl();
-                exit(EXIT_FAILURE);
+                throw CanteraError("vcs_VolPhase::setExistence",
+                                   "setting false existence for phase with moles");
             } else {
                 v_totalMoles = 0.0;
             }
@@ -1029,18 +1018,16 @@ void vcs_VolPhase::setExistence(const int existence)
     } else if (DEBUG_MODE_ENABLED && m_totalMolesInert == 0.0) {
         if (v_totalMoles == 0.0) {
             if (!m_singleSpecies  || m_phiVarIndex != 0) {
-                plogf("vcs_VolPhase::setExistence setting true existence for phase with no moles");
-                plogendl();
-                exit(EXIT_FAILURE);
+                throw CanteraError("vcs_VolPhase::setExistence",
+                        "setting true existence for phase with no moles");
             }
         }
     }
     if (DEBUG_MODE_ENABLED && m_singleSpecies) {
         if (m_phiVarIndex == 0) {
             if (existence == VCS_PHASE_EXIST_NO || existence == VCS_PHASE_EXIST_ZEROEDPHASE) {
-                plogf("vcs_VolPhase::Trying to set existence of an electron phase to false");
-                plogendl();
-                exit(EXIT_FAILURE);
+                throw CanteraError("vcs_VolPhase::setExistence",
+                        "Trying to set existence of an electron phase to false");
             }
         }
     }
diff --git a/src/equil/vcs_elem_rearrange.cpp b/src/equil/vcs_elem_rearrange.cpp
index 8174e9468..38747aaab 100644
--- a/src/equil/vcs_elem_rearrange.cpp
+++ b/src/equil/vcs_elem_rearrange.cpp
@@ -12,6 +12,9 @@
 #include "cantera/equil/vcs_solve.h"
 #include "cantera/equil/vcs_internal.h"
 #include "cantera/equil/vcs_VolPhase.h"
+#include "cantera/base/stringUtils.h"
+
+using namespace Cantera;
 
 namespace VCSnonideal
 {
@@ -76,9 +79,8 @@ int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa,
                 }
             }
             if (k == m_numElemConstraints) {
-                plogf("vcs_elem_rearrange::Shouldn't be here. Algorithm misfired.");
-                plogendl();
-                exit(EXIT_FAILURE);
+                throw CanteraError("vcs_elem_rearrange",
+                        "Shouldn't be here. Algorithm misfired.");
             }
 
             /*
@@ -176,13 +178,9 @@ void VCS_SOLVE::vcs_switch_elem_pos(size_t ipos, size_t jpos)
     if (ipos == jpos) {
         return;
     }
-    if (DEBUG_MODE_ENABLED && (ipos > (m_numElemConstraints - 1) ||
-                               jpos > (m_numElemConstraints - 1))) {
-        plogf("vcs_switch_elem_pos: ifunc = 0: inappropriate args: %d %d\n",
-              ipos, jpos);
-        plogendl();
-        exit(EXIT_FAILURE);
-    }
+    AssertThrowMsg(ipos < m_numElemConstraints && jpos < m_numElemConstraints,
+                   "vcs_switch_elem_pos",
+                   "inappropriate args: " + int2str(ipos) + " " + int2str(jpos));
     /*
      * Change the element Global Index list in each vcs_VolPhase object
      * to reflect the switch in the element positions.
diff --git a/src/equil/vcs_nondim.cpp b/src/equil/vcs_nondim.cpp
index e17e8675a..ac8100e87 100644
--- a/src/equil/vcs_nondim.cpp
+++ b/src/equil/vcs_nondim.cpp
@@ -14,6 +14,8 @@
 #include "cantera/base/stringUtils.h"
 #include "cantera/base/ctexceptions.h"
 
+using namespace Cantera;
+
 namespace VCSnonideal
 {
 double VCS_SOLVE::vcs_nondim_Farad(int mu_units, double TKelvin) const
@@ -32,9 +34,8 @@ double VCS_SOLVE::vcs_nondim_Farad(int mu_units, double TKelvin) const
     case VCS_UNITS_KELVIN:
         return Cantera::ElectronCharge * Cantera::Avogadro/ TKelvin;
     default:
-        plogf("vcs_nondim_Farad error: unknown units: %d\n", mu_units);
-        plogendl();
-        exit(EXIT_FAILURE);
+        throw CanteraError("vcs_nondim_Farad",
+                           "unknown units: " + int2str(mu_units));
     }
 }
 
@@ -55,9 +56,8 @@ double VCS_SOLVE::vcs_nondimMult_TP(int mu_units, double TKelvin) const
     case VCS_UNITS_MKS:
         return TKelvin * Cantera::GasConstant;
     default:
-        plogf("vcs_nondimMult_TP error: unknown units: %d\n", mu_units);
-        plogendl();
-        exit(EXIT_FAILURE);
+        throw CanteraError("vcs_nondimMult_TP",
+                           "unknown units: " + int2str(mu_units));
     }
 }
 
@@ -104,12 +104,9 @@ void VCS_SOLVE::vcs_nondim_TP()
          * reasonable input would be between these two numbers below.
          */
         if (tmole_orig < 1.0E-200 || tmole_orig > 1.0E200) {
-            plogf(" VCS_SOLVE::vcs_nondim_TP ERROR: Total input moles , %g,  is outside the range handled by vcs. exit",
-                  tmole_orig);
-            plogendl();
-            throw Cantera::CanteraError("VCS_SOLVE::vcs_nondim_TP",
-                                        " Total input moles ," + Cantera::fp2str(tmole_orig) +
-                                        "is outside the range handled by vcs.\n");
+            throw CanteraError("VCS_SOLVE::vcs_nondim_TP",
+                               "Total input moles ," + fp2str(tmole_orig) +
+                               "is outside the range handled by vcs.\n");
         }
 
         // Determine the scale of the problem
@@ -212,8 +209,7 @@ void VCS_SOLVE::vcs_printChemPotUnits(int unitsFormat) const
         plogf("J/kmol");
         break;
     default:
-        plogf("unknown units!");
-        exit(EXIT_FAILURE);
+        throw CanteraError("VCS_SOLVE::vcs_printChemPotUnits", "unknown units!");
     }
 }
 
diff --git a/src/equil/vcs_phaseStability.cpp b/src/equil/vcs_phaseStability.cpp
index 33b5b0b45..96291d982 100644
--- a/src/equil/vcs_phaseStability.cpp
+++ b/src/equil/vcs_phaseStability.cpp
@@ -7,6 +7,7 @@
 #include "cantera/equil/vcs_internal.h"
 #include "cantera/equil/vcs_species_thermo.h"
 #include "cantera/equil/vcs_VolPhase.h"
+#include "cantera/base/stringUtils.h"
 
 #include <cstdio>
 #include <algorithm>
@@ -36,13 +37,8 @@ bool VCS_SOLVE::vcs_popPhasePossible(const size_t iphasePop) const
 {
     vcs_VolPhase* Vphase = m_VolPhaseList[iphasePop];
 
-#ifdef DEBUG_MODE
-    int existence = Vphase->exists();
-    if (existence > 0) {
-        printf("ERROR vcs_popPhasePossible called for a phase that exists!");
-        std::exit(-1);
-    }
-#endif
+    AssertThrowMsg(!Vphase->exists(), "VCS_SOLVE::vcs_popPhasePossible",
+                   "called for a phase that exists!");
 
     /*
      * Loop through all of the species in the phase. We say the phase
@@ -52,11 +48,10 @@ bool VCS_SOLVE::vcs_popPhasePossible(const size_t iphasePop) const
      */
     for (size_t k = 0; k < Vphase->nSpecies(); k++) {
         size_t kspec = Vphase->spGlobalIndexVCS(k);
-        if (DEBUG_MODE_ENABLED && m_molNumSpecies_old[kspec] > 0.0) {
-            printf("ERROR vcs_popPhasePossible we shouldn't be here %lu %g > 0.0",
-                   kspec, m_molNumSpecies_old[kspec]);
-            exit(-1);
-        }
+        AssertThrowMsg(m_molNumSpecies_old[kspec] <= 0.0,
+                       "VCS_SOLVE::vcs_popPhasePossible",
+                       "we shouldn't be here " + int2str(kspec) + " "+
+                       fp2str(m_molNumSpecies_old[kspec]) + " > 0.0");
         size_t irxn = kspec - m_numComponents;
         if (kspec >= m_numComponents) {
             bool iPopPossible = true;
@@ -304,10 +299,8 @@ size_t VCS_SOLVE::vcs_popPhaseID(std::vector<size_t> & phasePopPhaseIDs)
                 }
                 if (DEBUG_MODE_ENABLED && Fephase < 0.0) {
                     strcpy(anote," (not stable)");
-                    if (m_tPhaseMoles_old[iph] > 0.0) {
-                        printf("shouldn't be here\n");
-                        exit(-1);
-                    }
+                    AssertThrowMsg(m_tPhaseMoles_old[iph] <= 0.0,
+                        "VCS_SOLVE::vcs_popPhaseID", "shouldn't be here");
                 }
 
                 if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
@@ -379,16 +372,12 @@ int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop)
     //   zeroed out within the algorithm.  We may later adjust the value.
     doublereal tPhaseMoles = 10. * m_totalMolNum * VCS_DELETE_PHASE_CUTOFF;
 
-#ifdef DEBUG_MODE
-    if (Vphase->exists() > 0) {
-        printf("ERROR vcs_popPhaseRxnStepSizes called for a phase that exists!");
-        exit(-1);
-    }
-    if (m_debug_print_lvl >= 2) {
+    AssertThrowMsg(!Vphase->exists(), "VCS_SOLVE::vcs_popPhaseRxnStepSizes",
+                   "called for a phase that exists!");
+    if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
         plogf("  ---  vcs_popPhaseRxnStepSizes() called to pop phase %s %d into existence\n",
               Vphase->PhaseName.c_str(), iphasePop);
     }
-#endif
     // Section for a single-species phase
     if (Vphase->m_singleSpecies) {
         double s = 0.0;
@@ -825,8 +814,7 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph)
 
 
     } else {
-        printf("not done yet\n");
-        exit(-1);
+        throw CanteraError("VCS_SOLVE::vcs_phaseStabilityTest", "not done yet");
     }
     if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
         plogf("  ------------------------------------------------------------"
diff --git a/src/equil/vcs_prob.cpp b/src/equil/vcs_prob.cpp
index 93520fb42..05601e9de 100644
--- a/src/equil/vcs_prob.cpp
+++ b/src/equil/vcs_prob.cpp
@@ -51,22 +51,22 @@ VCS_PROB::VCS_PROB(size_t nsp, size_t nel, size_t nph) :
 {
     NSPECIES0 = nspecies;
     if (nspecies <= 0) {
-        plogf("number of species is zero or neg\n");
-        exit(EXIT_FAILURE);
+        throw CanteraError("VCS_PROB::VCS_PROB",
+                           "number of species is zero or neg");
     }
     NE0       = ne;
     if (ne <= 0) {
-        plogf("number of elements is zero or neg\n");
-        exit(EXIT_FAILURE);
+        throw CanteraError("VCS_PROB::VCS_PROB",
+                           "number of elements is zero or neg");
     }
     NPHASE0   = NPhase;
     if (NPhase <= 0) {
-        plogf("number of phases is zero or neg\n");
-        exit(EXIT_FAILURE);
+        throw CanteraError("VCS_PROB::VCS_PROB",
+                           "number of phases is zero or neg");
     }
     if (nspecies < NPhase) {
-        plogf("number of species is less than number of phases\n");
-        exit(EXIT_FAILURE);
+        throw CanteraError("VCS_PROB::VCS_PROB",
+                           "number of species is less than number of phases");
     }
 
     m_gibbsSpecies.resize(nspecies, 0.0);
@@ -87,8 +87,8 @@ VCS_PROB::VCS_PROB(size_t nsp, size_t nel, size_t nph) :
     for (size_t kspec = 0; kspec < nspecies; kspec++) {
         VCS_SPECIES_THERMO* ts_tmp = new VCS_SPECIES_THERMO(0, 0);
         if (ts_tmp == 0) {
-            plogf("Failed to init a ts struct\n");
-            exit(EXIT_FAILURE);
+            throw CanteraError("VCS_PROB::VCS_PROB",
+                               "Failed to init a ts struct");
         }
         SpeciesThermo[kspec] = ts_tmp;
     }
@@ -132,9 +132,7 @@ void VCS_PROB::resizeSpecies(size_t nsp, int force)
         Charge.resize(nsp, 0.0);
         NSPECIES0 = nsp;
         if (nspecies > NSPECIES0) {
-            nspecies = NSPECIES0;
-            plogf("shouldn't be here\n");
-            exit(EXIT_FAILURE);
+            throw CanteraError("VCS_PROB::resizeSpecies", "shouldn't be here");
         }
     }
 }
@@ -186,8 +184,7 @@ void VCS_PROB::prob_report(int print_lvl)
 
             plogf("\t\tPres = %g atm\n", pres_atm);
         } else {
-            plogf("\tUnknown problem type\n");
-            exit(EXIT_FAILURE);
+            throw CanteraError("VCS_PROB::prob_report", "Unknown problem type");
         }
         plogf("\n");
         plogf("             Phase IDs of species\n");
@@ -321,8 +318,8 @@ void VCS_PROB::addPhaseElements(vcs_VolPhase* volPhase)
 size_t VCS_PROB::addElement(const char* elNameNew, int elType, int elactive)
 {
     if (!elNameNew) {
-        plogf("error: element must have a name\n");
-        exit(EXIT_FAILURE);
+        throw CanteraError("VCS_PROB::addElement",
+                           "error: element must have a name");
     }
     size_t nel = ne + 1;
     resizeElements(nel, 1);
@@ -339,15 +336,13 @@ size_t VCS_PROB::addOnePhaseSpecies(vcs_VolPhase* volPhase, size_t k, size_t kT)
         /*
          * Need to expand the number of species here
          */
-        plogf("Shouldn't be here\n");
-        exit(EXIT_FAILURE);
+        throw CanteraError("VCS_PROB::addOnePhaseSpecies", "Shouldn't be here");
     }
     const Cantera::Array2D& fm = volPhase->getFormulaMatrix();
     for (size_t eVP = 0; eVP < volPhase->nElemConstraints(); eVP++) {
         size_t e = volPhase->elemGlobalIndex(eVP);
-        if (DEBUG_MODE_ENABLED && e == npos) {
-            exit(EXIT_FAILURE);
-        }
+        AssertThrowMsg(e != npos, "VCS_PROB::addOnePhaseSpecies",
+                       "element not found");
         FormulaMatrix(kT,e) = fm(k,eVP);
     }
     /*
@@ -362,8 +357,7 @@ void VCS_PROB::reportCSV(const std::string& reportFile)
 {
     FILE* FP = fopen(reportFile.c_str(), "w");
     if (!FP) {
-        plogf("Failure to open file\n");
-        exit(EXIT_FAILURE);
+        throw CanteraError("VCS_PROB::reportCSV", "Failure to open file");
     }
 
     std::vector<double> volPM(nspecies, 0.0);
@@ -492,10 +486,8 @@ void VCS_PROB::reportCSV(const std::string& reportFile)
             tp->getChemPotentials(VCS_DATA_PTR(m_gibbsSpecies)+istart);
             for (size_t k = 0; k < nSpeciesPhase; k++) {
                 if (!vcs_doubleEqual(m_gibbsSpecies[istart+k], mu[k])) {
-                    fprintf(FP,"ERROR: incompatibility!\n");
                     fclose(FP);
-                    plogf("ERROR: incompatibility!\n");
-                    exit(EXIT_FAILURE);
+                    throw CanteraError("VCS_PROB::reportCSV", "incompatibility");
                 }
             }
         }
diff --git a/src/equil/vcs_report.cpp b/src/equil/vcs_report.cpp
index 4931f630b..a665a86b9 100644
--- a/src/equil/vcs_report.cpp
+++ b/src/equil/vcs_report.cpp
@@ -119,8 +119,7 @@ int VCS_SOLVE::vcs_report(int iconv)
             plogf("        NA         %14.7E    %12.4E", 1.0, m_feSpecies_old[l]);
             plogf("   Voltage = %14.7E", m_molNumSpecies_old[l] * molScale);
         } else {
-            plogf("we have a problem\n");
-            exit(EXIT_FAILURE);
+            throw CanteraError("VCS_SOLVE::vcs_report", "we have a problem");
         }
         plogf("\n");
     }
@@ -228,8 +227,7 @@ int VCS_SOLVE::vcs_report(int iconv)
         totalMoles +=  m_tPhaseMoles_old[iphase];
         if (m_tPhaseMoles_old[iphase] != VPhase->totalMoles()) {
             if (! vcs_doubleEqual(m_tPhaseMoles_old[iphase], VPhase->totalMoles())) {
-                plogf("We have a problem\n");
-                exit(EXIT_FAILURE);
+                throw CanteraError("VCS_SOLVE::vcs_report", "we have a problem");
             }
         }
         vcs_elabPhase(iphase, VCS_DATA_PTR(gaPhase));
@@ -320,8 +318,8 @@ int VCS_SOLVE::vcs_report(int iconv)
         double tmp = m_SSfeSpecies[l] + log(m_actCoeffSpecies_old[l])
                      + lx - m_lnMnaughtSpecies[l] + eContrib;
         if (fabs(m_feSpecies_old[l] - tmp) > 1.0E-7) {
-            plogf("\n\t\twe have a problem - doesn't add up\n");
-            exit(EXIT_FAILURE);
+            throw CanteraError("VCS_SOLVE::vcs_report",
+                               "we have a problem - doesn't add up");
         }
         plogf(" %12.4E |", m_feSpecies_old[l]);
         if (m_lnMnaughtSpecies[l] != 0.0) {
diff --git a/src/equil/vcs_rxnadj.cpp b/src/equil/vcs_rxnadj.cpp
index 505b84115..01fec624b 100644
--- a/src/equil/vcs_rxnadj.cpp
+++ b/src/equil/vcs_rxnadj.cpp
@@ -581,8 +581,8 @@ double VCS_SOLVE::vcs_Hessian_diag_adj(size_t irxn, double hessianDiag_Ideal)
     double diag = hessianDiag_Ideal;
     double hessActCoef = vcs_Hessian_actCoeff_diag(irxn);
     if (hessianDiag_Ideal <= 0.0) {
-        plogf("vcs_Hessian_diag_adj::We shouldn't be here\n");
-        exit(EXIT_FAILURE);
+        throw CanteraError("VCS_SOLVE::vcs_Hessian_diag_adj",
+                           "We shouldn't be here");
     }
     if (hessActCoef >= 0.0) {
         diag += hessActCoef;
diff --git a/src/equil/vcs_solve.cpp b/src/equil/vcs_solve.cpp
index 2184929bc..1b5e42794 100644
--- a/src/equil/vcs_solve.cpp
+++ b/src/equil/vcs_solve.cpp
@@ -10,6 +10,7 @@
 
 #include "cantera/equil/vcs_solve.h"
 #include "cantera/base/ctexceptions.h"
+#include "cantera/base/stringUtils.h"
 #include "cantera/equil/vcs_internal.h"
 #include "cantera/equil/vcs_prob.h"
 
@@ -20,6 +21,7 @@
 #include "cantera/base/clockWC.h"
 
 using namespace std;
+using namespace Cantera;
 
 namespace VCSnonideal
 {
@@ -654,8 +656,8 @@ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub)
         if (!strncmp(m_elementName[i].c_str(), "cn_", 3)) {
             m_elType[i] = VCS_ELEM_TYPE_CHARGENEUTRALITY;
             if (pub->m_elType[i] != VCS_ELEM_TYPE_CHARGENEUTRALITY) {
-                plogf("we have an inconsistency!\n");
-                exit(EXIT_FAILURE);
+                throw CanteraError("VCS_SOLVE::vcs_prob_specifyFully",
+                                   "we have an inconsistency!");
             }
         }
     }
@@ -664,9 +666,10 @@ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub)
         if (m_elType[i] ==  VCS_ELEM_TYPE_CHARGENEUTRALITY) {
             if (m_elemAbundancesGoal[i] != 0.0) {
                 if (fabs(m_elemAbundancesGoal[i]) > 1.0E-9) {
-                    plogf("Charge neutrality condition %s is signicantly nonzero, %g. Giving up\n",
-                          m_elementName[i].c_str(), m_elemAbundancesGoal[i]);
-                    exit(EXIT_FAILURE);
+                    throw CanteraError("VCS_SOLVE::vcs_prob_specifyFully",
+                            "Charge neutrality condition " + m_elementName[i] +
+                            " is signicantly nonzero, " + fp2str(m_elemAbundancesGoal[i]) +
+                            ". Giving up");
                 } else {
                     if (m_debug_print_lvl >= 2) {
                         plogf("Charge neutrality condition %s not zero, %g. Setting it zero\n",
@@ -932,25 +935,26 @@ int VCS_SOLVE::vcs_prob_update(VCS_PROB* pub)
                 k1 = vPhase->spGlobalIndexVCS(k);
                 double tmp = m_molNumSpecies_old[k1];
                 if (! vcs_doubleEqual(pubPhase->electricPotential() , tmp)) {
-                    plogf("We have an inconsistency in voltage, %g, %g\n",
-                          pubPhase->electricPotential(), tmp);
-                    exit(EXIT_FAILURE);
+                    throw CanteraError("VCS_SOLVE::vcs_prob_update",
+                            "We have an inconsistency in voltage, " +
+                            fp2str(pubPhase->electricPotential()) + " " +
+                            fp2str(tmp));
                 }
             }
 
             if (! vcs_doubleEqual(pub->mf[kT], vPhase->molefraction(k))) {
-                plogf("We have an inconsistency in mole fraction, %g, %g\n",
-                      pub->mf[kT], vPhase->molefraction(k));
-                exit(EXIT_FAILURE);
+                throw CanteraError("VCS_SOLVE::vcs_prob_update",
+                        "We have an inconsistency in mole fraction, " +
+                        fp2str(pub->mf[kT]) + " " + fp2str(vPhase->molefraction(k)));
             }
             if (pubPhase->speciesUnknownType(k) != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
                 sumMoles +=  pub->w[kT];
             }
         }
         if (! vcs_doubleEqual(sumMoles, vPhase->totalMoles())) {
-            plogf("We have an inconsistency in total moles, %g %g\n",
-                  sumMoles, pubPhase->totalMoles());
-            exit(EXIT_FAILURE);
+            throw CanteraError("VCS_SOLVE::vcs_prob_update",
+                            "We have an inconsistency in total moles, " +
+                            fp2str(sumMoles) + " " + fp2str(pubPhase->totalMoles()));
         }
 
     }
diff --git a/src/equil/vcs_solve_TP.cpp b/src/equil/vcs_solve_TP.cpp
index 7f68d18ed..36f647286 100644
--- a/src/equil/vcs_solve_TP.cpp
+++ b/src/equil/vcs_solve_TP.cpp
@@ -44,8 +44,8 @@ void VCS_SOLVE::checkDelta1(double* const dsLocal,
     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);
+            throw CanteraError("VCS_SOLVE::checkDelta1",
+                               "we have found a problem");
         }
     }
 }
@@ -558,14 +558,14 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
         }
         if (kspec >= m_numComponents) {
             if (m_molNumSpecies_new[m_numSpeciesTot] != 0.0) {
-                printf("vcs_solve_tp:: we shouldn't be here!\n");
-                exit(EXIT_FAILURE);
+                throw CanteraError("VCS_SOLVE::solve_tp_inner",
+                                   "we shouldn't be here!");
             }
             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);
+                throw CanteraError("VCS_SOLVE::solve_tp_inner",
+                                   "we shouldn't be here!");
             }
             ++m_numRxnMinorZeroed;
             allMinorZeroedSpecies = (m_numRxnMinorZeroed == m_numRxnRdc);
@@ -981,14 +981,12 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
                  *         to the reaction delta that we just computed.
                  *         This should keep the amount of material constant.
                  */
-                if (DEBUG_MODE_ENABLED && 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);
-                }
+                AssertThrowMsg(fabs(m_deltaMolNumSpecies[kspec] -dx) <
+                        1.0E-14*(fabs(m_deltaMolNumSpecies[kspec]) + fabs(dx) + 1.0E-32),
+                        "VCS_SOLVE::solve_tp_inner",
+                        "ds[kspec] = " + fp2str(m_deltaMolNumSpecies[kspec]) +
+                        " dx = " + fp2str(dx) + " , kspec = " + int2str(kspec) +
+                        "\nwe have a problem!");
                 for (size_t k = 0; k < m_numComponents; ++k) {
                     m_deltaMolNumSpecies[k] += sc_irxn[k] * dx;
                 }
@@ -1104,10 +1102,10 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
             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);
+                throw CanteraError("VCS_SOLVE::solve_tp_inner",
+                    "vcs_solve_TP: ERROR on step change wt[" + int2str(kspec) + ":" +
+                    m_speciesName[kspec] + "]: " +
+                    fp2str(m_molNumSpecies_new[kspec]) + " < 0.0");
             }
         }
 
@@ -1330,10 +1328,8 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
                               VCS_DATA_PTR(m_sm), VCS_DATA_PTR(m_ss), test,
                               &usedZeroedSpecies);
         if (retn != VCS_SUCCESS) {
-            if (DEBUG_MODE_ENABLED) {
-                plogf("   --- BASOPT returned with an error condition\n");
-            }
-            exit(EXIT_FAILURE);
+            throw CanteraError("VCS_SOLVE::solve_tp_inner",
+                               "BASOPT returned with an error condition");
         }
         vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
         vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc);
@@ -1818,11 +1814,8 @@ int VCS_SOLVE::delta_species(const size_t kspec, double* const delta_ptr)
     size_t irxn = kspec - m_numComponents;
     int retn = 1;
     double delta = *delta_ptr;
-    if (DEBUG_MODE_ENABLED && kspec < m_numComponents) {
-        plogf("  --- delete_species() ERROR: called for a component %d", kspec);
-        plogendl();
-        exit(EXIT_FAILURE);
-    }
+    AssertThrowMsg(kspec >= m_numComponents, "VCS_SOLVE::delta_species",
+        "delete_species() ERROR: called for a component " + int2str(kspec));
     if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
         /*
          * Attempt the given dx. If it doesn't work, try to see if a smaller
@@ -1904,11 +1897,8 @@ int VCS_SOLVE::vcs_delete_species(const size_t kspec)
      *     -> 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);
-    }
+    AssertThrowMsg(retn, "VCS_SOLVE::vcs_delete_species",
+        "Failed to delete a species!");
     /*
      *    Decrement the minor species counter if the current species is
      *    a minor species
@@ -2717,8 +2707,7 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[],
                 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);
+                    throw CanteraError("VCS_SOLVE::vcs_basopt", "we shouldn't be here");
                 }
                 m_numRxnTot = m_numSpeciesTot - ncTrial;
                 m_numRxnRdc = m_numRxnTot - numPreDeleted;
@@ -3022,8 +3011,7 @@ L_END_LOOP:
                     }
                 }
                 if (fabs(sum) > 1.0E-6) {
-                    printf("we have a prob\n");
-                    exit(-1);
+                    throw CanteraError("VCS_SOLVE::vcs_basopt", "we have a prob");
                 }
             }
         }
@@ -3361,16 +3349,10 @@ void VCS_SOLVE::vcs_chemPotPhase(const int stateCalc,
             }
         }
         if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
-            if (DEBUG_MODE_ENABLED) {
-                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);
-                }
-            }
+            AssertThrowMsg(molNum[kspec] == phi, "VCS_SOLVE::vcs_chemPotPhase",
+                "We have an inconsistency!");
+            AssertThrowMsg(m_chargeSpecies[kspec] == -1.0, "VCS_SOLVE::vcs_chemPotPhase",
+                "We have an unexpected situation!");
             mu_i[kspec] = m_SSfeSpecies[kspec] + m_chargeSpecies[kspec] * Faraday_phi;
         } else {
             if (m_SSPhase[kspec]) {
@@ -3404,16 +3386,13 @@ void VCS_SOLVE::vcs_dfe(const int stateCalc,
         actCoeff_ptr = VCS_DATA_PTR(m_actCoeffSpecies_new);
         molNum = VCS_DATA_PTR(m_molNumSpecies_new);
     } else if (DEBUG_MODE_ENABLED) {
-        plogf("vcs_dfe: wrong stateCalc value");
-        plogf("   --- Subroutine vcs_dfe called with bad stateCalc value: %d", stateCalc);
-        plogendl();
-        exit(EXIT_FAILURE);
+        throw CanteraError("VCS_SOLVE::vcs_dfe",
+                "Subroutine vcs_dfe called with bad stateCalc value: "+
+                int2str(stateCalc));
     }
 
-    if (DEBUG_MODE_ENABLED && m_unitsState == VCS_DIMENSIONAL_G) {
-        printf("vcs_dfe: called with wrong units state\n");
-        exit(EXIT_FAILURE);
-    }
+    AssertThrowMsg(m_unitsState != VCS_DIMENSIONAL_G, "VCS_SOLVE::vcs_dfe",
+        "called with wrong units state");
 
     if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
         if (ll == 0) {
@@ -3450,14 +3429,11 @@ void VCS_SOLVE::vcs_dfe(const int stateCalc,
             tlogMoles[iph] += molNum[kspec];
         }
     }
-    if (DEBUG_MODE_ENABLED) {
-        for (size_t 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);
-            }
-        }
+    for (size_t iph = 0; iph < m_numPhases; iph++) {
+        AssertThrowMsg(vcs_doubleEqual(tlogMoles[iph], tPhMoles_ptr[iph]),
+            "VCS_SOLVE::vcs_dfe",
+            "phase Moles may be off, iph = " + int2str(iph) + ", " +
+            fp2str(tlogMoles[iph]) + " " + fp2str(tPhMoles_ptr[iph]));
     }
     m_TmpPhase.assign(m_TmpPhase.size(), 0.0);
     for (size_t iph = 0; iph < m_numPhases; iph++) {
@@ -3500,15 +3476,10 @@ void VCS_SOLVE::vcs_dfe(const int stateCalc,
     for (size_t kspec = l1; kspec < l2; ++kspec) {
         size_t iphase = m_phaseID[kspec];
         if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
-            if (DEBUG_MODE_ENABLED) {
-                if (molNum[kspec] != m_phasePhi[iphase]) {
-                    plogf("We have an inconsistency!\n");
-                    exit(EXIT_FAILURE);
-                } else if (m_chargeSpecies[kspec] != -1.0) {
-                    plogf("We have an unexpected situation!\n");
-                    exit(EXIT_FAILURE);
-                }
-            }
+            AssertThrowMsg(molNum[kspec] == m_phasePhi[iphase], "VCS_SOLVE::vcs_dfe",
+                "We have an inconsistency!");
+            AssertThrowMsg(m_chargeSpecies[kspec] == -1.0, "VCS_SOLVE::vcs_dfe",
+                "We have an unexpected situation!");
             feSpecies[kspec] = m_SSfeSpecies[kspec]
                                + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
         } else {
@@ -3549,15 +3520,10 @@ void VCS_SOLVE::vcs_dfe(const int stateCalc,
             if (m_speciesStatus[kspec] != VCS_SPECIES_MINOR) {
                 size_t iphase = m_phaseID[kspec];
                 if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
-                    if (DEBUG_MODE_ENABLED) {
-                        if (molNum[kspec] != m_phasePhi[iphase]) {
-                            plogf("We have an inconsistency!\n");
-                            exit(EXIT_FAILURE);
-                        } else if (m_chargeSpecies[kspec] != -1.0) {
-                            plogf("We have an unexpected situation!\n");
-                            exit(EXIT_FAILURE);
-                        }
-                    }
+                    AssertThrowMsg(molNum[kspec] == m_phasePhi[iphase], "VCS_SOLVE::vcs_dfe",
+                                   "We have an inconsistency!");
+                    AssertThrowMsg(m_chargeSpecies[kspec] == -1.0, "VCS_SOLVE::vcs_dfe",
+                                   "We have an unexpected situation!");
                     feSpecies[kspec] = m_SSfeSpecies[kspec]
                                        + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
                 } else {
@@ -3599,15 +3565,10 @@ void VCS_SOLVE::vcs_dfe(const int stateCalc,
             if (m_speciesStatus[kspec] == VCS_SPECIES_MINOR) {
                 size_t iphase = m_phaseID[kspec];
                 if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
-                    if (DEBUG_MODE_ENABLED) {
-                        if (molNum[kspec] != m_phasePhi[iphase]) {
-                            plogf("We have an inconsistency!\n");
-                            exit(EXIT_FAILURE);
-                        } else if (m_chargeSpecies[kspec] != -1.0) {
-                            plogf("We have an unexpected situation!\n");
-                            exit(EXIT_FAILURE);
-                        }
-                    }
+                    AssertThrowMsg(molNum[kspec] == m_phasePhi[iphase], "VCS_SOLVE::vcs_dfe",
+                        "We have an inconsistency!");
+                    AssertThrowMsg(m_chargeSpecies[kspec] == -1.0, "VCS_SOLVE::vcs_dfe",
+                        "We have an unexpected situation!");
                     feSpecies[kspec] = m_SSfeSpecies[kspec]
                                        + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase]; ;
                 } else {
@@ -3839,9 +3800,8 @@ void VCS_SOLVE::vcs_updateVP(const int vcsState)
                                          VCS_DATA_PTR(m_molNumSpecies_new),
                                          VCS_DATA_PTR(m_tPhaseMoles_new));
         } else if (DEBUG_MODE_ENABLED) {
-            plogf("vcs_updateVP ERROR: wrong stateCalc value: %d", vcsState);
-            plogendl();
-            exit(EXIT_FAILURE);
+            throw CanteraError("VCS_SOLVE::vcs_updateVP",
+                               "wrong stateCalc value: " + int2str(vcsState));
         }
     }
 }
@@ -3901,9 +3861,8 @@ bool VCS_SOLVE::vcs_evaluate_speciesType()
                     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);
+                    throw CanteraError("VCS_SOLVE::vcs_evaluate_speciesType",
+                                       "Unknown type: " + int2str(m_speciesStatus[kspec]));
                 }
             }
         }
@@ -3944,8 +3903,7 @@ void VCS_SOLVE::vcs_deltag(const int l, const bool doDeleted,
         molNumSpecies = VCS_DATA_PTR(m_molNumSpecies_old);
         actCoeffSpecies = VCS_DATA_PTR(m_actCoeffSpecies_old);
     } else {
-        printf("Error\n");
-        exit(EXIT_FAILURE);
+        throw CanteraError("VCS_SOLVE::vcs_deltag", "bad vcsState");
     }
 
     if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
@@ -4265,9 +4223,7 @@ void VCS_SOLVE::vcs_deltag_Phase(const size_t iphase, const bool doDeleted,
         deltaGRxn = VCS_DATA_PTR(m_deltaGRxn_old);
         actCoeffSpecies = VCS_DATA_PTR(m_actCoeffSpecies_old);
     } else if (DEBUG_MODE_ENABLED) {
-        plogf("vcs_deltag_Phase: we shouldn't be here\n");
-        plogendl();
-        exit(EXIT_FAILURE);
+        throw CanteraError("VCS_SOLVE::vcs_deltag_Phase", "bad stateCalc");
     }
 
     size_t irxnl = m_numRxnRdc;
@@ -4286,10 +4242,8 @@ void VCS_SOLVE::vcs_deltag_Phase(const size_t iphase, const bool doDeleted,
      */
     if (vPhase->m_singleSpecies) {
         size_t kspec = vPhase->spGlobalIndexVCS(0);
-        if (DEBUG_MODE_ENABLED && iphase != m_phaseID[kspec]) {
-            plogf("vcs_deltag_Phase index error\n");
-            exit(EXIT_FAILURE);
-        }
+        AssertThrowMsg(iphase == m_phaseID[kspec], "VCS_SOLVE::vcs_deltag_Phase",
+            "index error");
         if (kspec >= m_numComponents) {
             size_t irxn = kspec - m_numComponents;
             deltaGRxn[irxn] = feSpecies[kspec];
@@ -4401,15 +4355,10 @@ void VCS_SOLVE::vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k
 
     size_t kp1 = m_speciesLocalPhaseIndex[k1];
     size_t kp2 = m_speciesLocalPhaseIndex[k2];
-    if (DEBUG_MODE_ENABLED) {
-        if (pv1->spGlobalIndexVCS(kp1) != k1) {
-            plogf("Indexing error in program\n");
-            exit(EXIT_FAILURE);
-        } else if (pv2->spGlobalIndexVCS(kp2) != k2) {
-            plogf("Indexing error in program\n");
-            exit(EXIT_FAILURE);
-        }
-    }
+    AssertThrowMsg(pv1->spGlobalIndexVCS(kp1) == k1, "VCS_SOLVE::vcs_switch_pos",
+        "Indexing error");
+    AssertThrowMsg(pv2->spGlobalIndexVCS(kp2) == k2, "VCS_SOLVE::vcs_switch_pos",
+        "Indexing error");
     pv1->setSpGlobalIndexVCS(kp1, k2);
     pv2->setSpGlobalIndexVCS(kp2, k1);
     std::swap(m_speciesName[k1], m_speciesName[k2]);
@@ -4487,11 +4436,8 @@ double VCS_SOLVE::vcs_birthGuess(const int kspec)
 
     // Check to make sure that species is zero in the solution vector
     // If it isn't, we don't know what's happening
-    if (DEBUG_MODE_ENABLED && m_molNumSpecies_old[kspec] != 0.0) {
-        w_kspec = 0.0;
-        plogf("vcs_birthGuess:: we shouldn't be here\n");
-        exit(EXIT_FAILURE);
-    }
+    AssertThrowMsg(m_molNumSpecies_old[kspec] == 0.0,
+            "VCS_SOLVE::vcs_birthGuess", "we shouldn't be here");
     int ss = m_SSPhase[kspec];
     if (!ss) {
         /*
diff --git a/src/equil/vcs_species_thermo.cpp b/src/equil/vcs_species_thermo.cpp
index 621516a71..67d61b657 100644
--- a/src/equil/vcs_species_thermo.cpp
+++ b/src/equil/vcs_species_thermo.cpp
@@ -17,6 +17,7 @@
 #include "cantera/equil/vcs_internal.h"
 
 using namespace std;
+using namespace Cantera;
 
 namespace VCSnonideal
 {
@@ -123,7 +124,6 @@ VCS_SPECIES_THERMO* VCS_SPECIES_THERMO::duplMyselfAsVCS_SPECIES_THERMO()
 double VCS_SPECIES_THERMO::GStar_R_calc(size_t kglob, double TKelvin,
                                         double pres)
 {
-    char yo[] = "VCS_SPECIES_THERMO::GStar_R_calc ";
     double fe = G0_R_calc(kglob, TKelvin);
     double T = TKelvin;
     if (UseCanteraCalls) {
@@ -145,8 +145,8 @@ double VCS_SPECIES_THERMO::GStar_R_calc(size_t kglob, double TKelvin,
             fe += T * log(pres/ pref);
             break;
         default:
-            plogf("%sERROR: unknown SSStar model\n", yo);
-            exit(EXIT_FAILURE);
+            throw CanteraError("VCS_SPECIES_THERMO::GStar_R_calc",
+                               "unknown SSStar model");
         }
     }
     return fe;
@@ -155,7 +155,6 @@ double VCS_SPECIES_THERMO::GStar_R_calc(size_t kglob, double TKelvin,
 double VCS_SPECIES_THERMO::VolStar_calc(size_t kglob, double TKelvin, 
                                         double presPA)
 {
-    char yo[] = "VCS_SPECIES_THERMO::VStar_calc ";
     double vol;
 
     double T = TKelvin;
@@ -176,8 +175,8 @@ double VCS_SPECIES_THERMO::VolStar_calc(size_t kglob, double TKelvin,
             vol= Cantera::GasConstant * T / presPA;
             break;
         default:
-            plogf("%sERROR: unknown SSVol model\n", yo);
-            exit(EXIT_FAILURE);
+            throw CanteraError("VCS_SPECIES_THERMO::VolStar_calc",
+                               "unknown SSVol model");
         }
     }
     return vol;
@@ -213,8 +212,8 @@ double VCS_SPECIES_THERMO::G0_R_calc(size_t kglob, double TKelvin)
             fe = H - TKelvin * S;
             break;
         default:
-            plogf("VS_SPECIES_THERMO::G0_R_calc ERROR: unknown model\n");
-            exit(EXIT_FAILURE);
+            throw CanteraError("VCS_SPECIES_THERMO::G0_R_calc",
+                               "unknown model");
         }
     }
     SS0_feSave = fe;
@@ -240,8 +239,7 @@ double VCS_SPECIES_THERMO::eval_ac(size_t kglob)
             ac = 1.0;
             break;
         default:
-            plogf("VCS_SPECIES_THERMO::eval_ac ERROR: unknown model\n");
-            exit(EXIT_FAILURE);
+            throw CanteraError("VCS_SPECIES_THERMO::eval_ac" ,"unknown model");
         }
     }
     return ac;
diff --git a/src/equil/vcs_util.cpp b/src/equil/vcs_util.cpp
index 99eef7689..4c7eee469 100644
--- a/src/equil/vcs_util.cpp
+++ b/src/equil/vcs_util.cpp
@@ -9,9 +9,12 @@
  */
 
 #include "cantera/equil/vcs_internal.h"
+#include "cantera/base/stringUtils.h"
+#include "cantera/base/ctexceptions.h"
 #include <cassert>
 
 using namespace std;
+using namespace Cantera;
 
 namespace VCSnonideal
 {
@@ -83,9 +86,8 @@ double vcsUtil_gasConstant(int mu_units)
         /* joules / kg-mol K = kg m2 / s2 kg-mol K */
         return Cantera::GasConstant;
     default:
-        plogf("vcs_gasConstant error: uknown units: %d\n",
-              mu_units);
-        exit(EXIT_FAILURE);
+        throw CanteraError("vcsUtil_gasConstant",
+                           "uknown units: " + int2str(mu_units));
     }
 }
 
diff --git a/src/kinetics/solveSP.cpp b/src/kinetics/solveSP.cpp
index 23233b0ee..3fef75415 100644
--- a/src/kinetics/solveSP.cpp
+++ b/src/kinetics/solveSP.cpp
@@ -823,8 +823,8 @@ void solveSP::print_header(int ioflag, int ifunc, doublereal time_scale,
             printf("\n   SOLVESP Called to integrate surface in time\n");
             printf("           for a total of %9.3e sec\n", time_scale);
         } else {
-            fprintf(stderr,"Unknown ifunc flag = %d\n", ifunc);
-            exit(EXIT_FAILURE);
+            throw CanteraError("solveSP::print_header",
+                               "Unknown ifunc flag = " + int2str(ifunc));
         }
 
         if (m_bulkFunc == BULK_DEPOSITION) {
@@ -832,8 +832,8 @@ void solveSP::print_header(int ioflag, int ifunc, doublereal time_scale,
         } else if (m_bulkFunc == BULK_ETCH) {
             printf("     Bulk Phases have fixed compositions\n");
         } else {
-            fprintf(stderr,"Unknown bulkFunc flag = %d\n", m_bulkFunc);
-            exit(EXIT_FAILURE);
+            throw CanteraError("solveSP::print_header",
+                               "Unknown bulkFunc flag = " + int2str(m_bulkFunc));
         }
 
         if (damping) {
diff --git a/src/numerics/solveProb.cpp b/src/numerics/solveProb.cpp
index e52dda964..c6b056fda 100644
--- a/src/numerics/solveProb.cpp
+++ b/src/numerics/solveProb.cpp
@@ -626,8 +626,8 @@ void solveProb::print_header(int ioflag, int ifunc, doublereal time_scale,
             printf("\n   SOLVEPROB Called to integrate surface in time\n");
             printf("           for a total of %9.3e sec\n", time_scale);
         } else {
-            fprintf(stderr,"Unknown ifunc flag = %d\n", ifunc);
-            exit(EXIT_FAILURE);
+            throw CanteraError("solveProb::print_header",
+                               "Unknown ifunc flag = " + int2str(ifunc));
         }
 
 
diff --git a/src/thermo/DebyeHuckel.cpp b/src/thermo/DebyeHuckel.cpp
index b4d92be29..7bbe3f682 100644
--- a/src/thermo/DebyeHuckel.cpp
+++ b/src/thermo/DebyeHuckel.cpp
@@ -644,12 +644,12 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
         if (formString != "") {
             if (formString == "unity") {
                 m_formGC = 0;
-                printf("exit standardConc = unity not done\n");
-                exit(EXIT_FAILURE);
+                throw CanteraError("DebyeHuckel::initThermoXML",
+                                   "standardConc = unity not done");
             } else if (formString == "molar_volume") {
                 m_formGC = 1;
-                printf("exit standardConc = molar_volume not done\n");
-                exit(EXIT_FAILURE);
+                throw CanteraError("DebyeHuckel::initThermoXML",
+                                   "standardConc = molar_volume not done");
             } else if (formString == "solvent_volume") {
                 m_formGC = 2;
             } else {
@@ -1106,8 +1106,7 @@ double DebyeHuckel::A_Debye_TP(double tempArg, double presArg) const
         m_A_Debye = A;
         break;
     default:
-        printf("shouldn't be here\n");
-        exit(EXIT_FAILURE);
+        throw CanteraError("DebyeHuckel::A_Debye_TP", "shouldn't be here");
     }
     return A;
 }
@@ -1131,8 +1130,7 @@ double DebyeHuckel::dA_DebyedT_TP(double tempArg, double presArg) const
         dAdT = m_waterProps->ADebye(T, P, 1);
         break;
     default:
-        printf("shouldn't be here\n");
-        exit(EXIT_FAILURE);
+        throw CanteraError("DebyeHuckel::dA_DebyedT_TP", "shouldn't be here");
     }
     return dAdT;
 }
@@ -1156,8 +1154,7 @@ double DebyeHuckel::d2A_DebyedT2_TP(double tempArg, double presArg) const
         d2AdT2 = m_waterProps->ADebye(T, P, 2);
         break;
     default:
-        printf("shouldn't be here\n");
-        exit(EXIT_FAILURE);
+        throw CanteraError("DebyeHuckel::d2A_DebyedT2_TP", "shouldn't be here");
     }
     return d2AdT2;
 }
@@ -1181,8 +1178,7 @@ double DebyeHuckel::dA_DebyedP_TP(double tempArg, double presArg) const
         dAdP = m_waterProps->ADebye(T, P, 3);
         break;
     default:
-        printf("shouldn't be here\n");
-        exit(EXIT_FAILURE);
+        throw CanteraError("DebyeHuckel::dA_DebyedP_TP", "shouldn't be here");
     }
     return dAdP;
 }
@@ -1510,8 +1506,7 @@ void DebyeHuckel::s_update_lnMolalityActCoeff() const
         break;
 
     default:
-        printf("ERROR\n");
-        exit(EXIT_FAILURE);
+        throw CanteraError("DebyeHuckel::s_update_lnMolalityActCoeff", "ERROR");
     }
     /*
      * Above, we calculated the ln(activitySolvent). Translate that
@@ -1643,9 +1638,8 @@ void DebyeHuckel::s_update_dlnMolalityActCoeff_dT() const
         break;
 
     default:
-        printf("ERROR\n");
-        exit(EXIT_FAILURE);
-        break;
+        throw CanteraError("DebyeHuckel::s_update_dlnMolalityActCoeff_dT",
+                           "ERROR");
     }
 
 
@@ -1767,9 +1761,8 @@ void DebyeHuckel::s_update_d2lnMolalityActCoeff_dT2() const
         break;
 
     default:
-        printf("ERROR\n");
-        exit(EXIT_FAILURE);
-        break;
+        throw CanteraError("DebyeHuckel::s_update_d2lnMolalityActCoeff_dT2",
+                           "ERROR");
     }
 }
 
@@ -1893,9 +1886,8 @@ void DebyeHuckel::s_update_dlnMolalityActCoeff_dP() const
         break;
 
     default:
-        printf("ERROR\n");
-        exit(EXIT_FAILURE);
-        break;
+        throw CanteraError("DebyeHuckel::s_update_dlnMolalityActCoeff_dP",
+                           "ERROR");
     }
 }
 
diff --git a/src/thermo/HMWSoln.cpp b/src/thermo/HMWSoln.cpp
index 6a5bf8bd6..1db86e514 100644
--- a/src/thermo/HMWSoln.cpp
+++ b/src/thermo/HMWSoln.cpp
@@ -439,8 +439,7 @@ HMWSoln::HMWSoln(int testProb) :
     m_debugCalc(0)
 {
     if (testProb != 1) {
-        printf("unknown test problem\n");
-        exit(EXIT_FAILURE);
+        throw CanteraError("HMWSoln::HMWSoln", "unknown test problem");
     }
 
     initThermoFile("HMW_NaCl.xml", "");
@@ -1016,8 +1015,7 @@ double HMWSoln::A_Debye_TP(double tempArg, double presArg) const
         m_A_Debye = A;
         break;
     default:
-        printf("shouldn't be here\n");
-        exit(EXIT_FAILURE);
+        throw CanteraError("HMWSoln::A_Debye_TP", "shouldn't be here");
     }
     return A;
 }
@@ -1041,8 +1039,7 @@ double HMWSoln::dA_DebyedT_TP(double tempArg, double presArg) const
         dAdT = m_waterProps->ADebye(T, P, 1);
         break;
     default:
-        printf("shouldn't be here\n");
-        exit(EXIT_FAILURE);
+        throw CanteraError("HMWSoln::dA_DebyedT_TP", "shouldn't be here");
     }
     return dAdT;
 }
@@ -1074,8 +1071,7 @@ double HMWSoln::dA_DebyedP_TP(double tempArg, double presArg) const
         }
         break;
     default:
-        printf("shouldn't be here\n");
-        exit(EXIT_FAILURE);
+        throw CanteraError("HMWSoln::dA_DebyedP_TP", "shouldn't be here");
     }
     return dAdP;
 }
@@ -1133,8 +1129,7 @@ double HMWSoln::d2A_DebyedT2_TP(double tempArg, double presArg) const
         d2AdT2 = m_waterProps->ADebye(T, P, 2);
         break;
     default:
-        printf("shouldn't be here\n");
-        exit(EXIT_FAILURE);
+        throw CanteraError("HMWSoln::d2A_DebyedT2_TP", "shouldn't be here");
     }
     return d2AdT2;
 }
@@ -1841,8 +1836,8 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
      *        species 0.
      */
     if (m_indexSolvent != 0) {
-        printf("Wrong index solvent value!\n");
-        exit(EXIT_FAILURE);
+        throw CanteraError("HMWSoln::s_updatePitzer_lnMolalityActCoeff",
+                           "Wrong index solvent value!");
     }
 
     std::string sni,  snj, snk;
@@ -2598,8 +2593,8 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
             for (size_t k = j+1; k < m_kk; k++) {
                 if (j == (m_kk-1)) {
                     // we should never reach this step
-                    printf("logic error 1 in Step 9 of hmw_act");
-                    exit(EXIT_FAILURE);
+                    throw CanteraError("HMWSoln::s_updatePitzer_lnMolalityActCoeff",
+                                       "logic error 1 in Step 9 of hmw_act");
                 }
                 if (charge(k) > 0.0) {
                     /*
@@ -2628,8 +2623,8 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
             for (size_t k = j+1; k < m_kk; k++) {
                 if (j == m_kk-1) {
                     // we should never reach this step
-                    printf("logic error 2 in Step 9 of hmw_act");
-                    exit(EXIT_FAILURE);
+                    throw CanteraError("HMWSoln::s_updatePitzer_lnMolalityActCoeff",
+                                       "logic error 2 in Step 9 of hmw_act");
                 }
                 if (charge(k) < 0) {
                     /*
@@ -2778,8 +2773,8 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
     m_debugCalc = 0;
 #endif
     if (m_indexSolvent != 0) {
-        printf("Wrong index solvent value!\n");
-        exit(EXIT_FAILURE);
+        throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT",
+                           "Wrong index solvent value!");
     }
 
     std::string sni, snj, snk;
@@ -3392,8 +3387,8 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
             for (size_t k = j+1; k < m_kk; k++) {
                 if (j == (m_kk-1)) {
                     // we should never reach this step
-                    printf("logic error 1 in Step 9 of hmw_act");
-                    exit(EXIT_FAILURE);
+                    throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT",
+                                       "logic error 1 in Step 9 of hmw_act");
                 }
                 if (charge(k) > 0.0) {
                     /*
@@ -3422,8 +3417,8 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
             for (size_t k = j+1; k < m_kk; k++) {
                 if (j == m_kk-1) {
                     // we should never reach this step
-                    printf("logic error 2 in Step 9 of hmw_act");
-                    exit(EXIT_FAILURE);
+                    throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT",
+                                       "logic error 2 in Step 9 of hmw_act");
                 }
                 if (charge(k) < 0) {
                     /*
@@ -3559,8 +3554,8 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
     m_debugCalc = 0;
 #endif
     if (m_indexSolvent != 0) {
-        printf("Wrong index solvent value!\n");
-        exit(EXIT_FAILURE);
+        throw CanteraError("HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2",
+                           "Wrong index solvent value!");
     }
 
     std::string sni, snj, snk;
@@ -4173,8 +4168,8 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
             for (size_t k = j+1; k < m_kk; k++) {
                 if (j == (m_kk-1)) {
                     // we should never reach this step
-                    printf("logic error 1 in Step 9 of hmw_act");
-                    exit(EXIT_FAILURE);
+                    throw CanteraError("HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2",
+                                       "logic error 1 in Step 9 of hmw_act");
                 }
                 if (charge(k) > 0.0) {
                     /*
@@ -4203,8 +4198,8 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
             for (size_t k = j+1; k < m_kk; k++) {
                 if (j == m_kk-1) {
                     // we should never reach this step
-                    printf("logic error 2 in Step 9 of hmw_act");
-                    exit(EXIT_FAILURE);
+                    throw CanteraError("HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2",
+                                       "logic error 2 in Step 9 of hmw_act");
                 }
                 if (charge(k) < 0) {
                     /*
@@ -4333,8 +4328,8 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
     m_debugCalc = 0;
 #endif
     if (m_indexSolvent != 0) {
-        printf("Wrong index solvent value!\n");
-        exit(EXIT_FAILURE);
+        throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP",
+                           "Wrong index solvent value!");
     }
 
     std::string sni, snj, snk;
@@ -4949,8 +4944,8 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
             for (size_t k = j+1; k < m_kk; k++) {
                 if (j == (m_kk-1)) {
                     // we should never reach this step
-                    printf("logic error 1 in Step 9 of hmw_act");
-                    exit(EXIT_FAILURE);
+                    throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP",
+                                       "logic error 1 in Step 9 of hmw_act");
                 }
                 if (charge(k) > 0.0) {
                     /*
@@ -4980,8 +4975,8 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
             for (size_t k = j+1; k < m_kk; k++) {
                 if (j == m_kk-1) {
                     // we should never reach this step
-                    printf("logic error 2 in Step 9 of hmw_act");
-                    exit(EXIT_FAILURE);
+                    throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP",
+                                       "logic error 2 in Step 9 of hmw_act");
                 }
                 if (charge(k) < 0) {
                     /*
@@ -5153,15 +5148,10 @@ void HMWSoln::calc_thetas(int z1, int z2,
     i = abs(z1);
     j = abs(z2);
 
-    if (DEBUG_MODE_ENABLED && (i > 4 || j > 4)) {
-        printf("we shouldn't be here\n");
-        exit(EXIT_FAILURE);
-    }
-
-    if ((i == 0) || (j == 0)) {
-        printf("ERROR calc_thetas called with one species being neutral\n");
-        exit(EXIT_FAILURE);
-    }
+    AssertThrowMsg(i <= 4 && j <= 4, "HMWSoln::calc_thetas",
+                   "we shouldn't be here");
+    AssertThrowMsg(i != 0 && j != 0, "HMWSoln::calc_thetas",
+                   "called with one species being neutral");
 
     /*
      *  Check to see if the charges are of opposite sign. If they are of
diff --git a/src/thermo/HMWSoln_input.cpp b/src/thermo/HMWSoln_input.cpp
index e43eb5a97..3633ffecf 100644
--- a/src/thermo/HMWSoln_input.cpp
+++ b/src/thermo/HMWSoln_input.cpp
@@ -1035,12 +1035,12 @@ void HMWSoln::constructPhaseXML(XML_Node& phaseNode, std::string id_)
         if (formString != "") {
             if (formString == "unity") {
                 m_formGC = 0;
-                printf("exit standardConc = unity not done\n");
-                exit(EXIT_FAILURE);
+                throw CanteraError("HMWSoln::constructPhaseXML",
+                                   "standardConc = unity not done");
             } else if (formString == "molar_volume") {
                 m_formGC = 1;
-                printf("exit standardConc = molar_volume not done\n");
-                exit(EXIT_FAILURE);
+                throw CanteraError("HMWSoln::constructPhaseXML",
+                   "standardConc = molar_volume not done");
             } else if (formString == "solvent_volume") {
                 m_formGC = 2;
             } else {
@@ -1162,12 +1162,12 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
         if (formString != "") {
             if (formString == "unity") {
                 m_formGC = 0;
-                printf("exit standardConc = unity not done\n");
-                exit(EXIT_FAILURE);
+                throw CanteraError("HMWSoln::initThermoXML",
+                                   "standardConc = unity not done");
             } else if (formString == "molar_volume") {
                 m_formGC = 1;
-                printf("exit standardConc = molar_volume not done\n");
-                exit(EXIT_FAILURE);
+                throw CanteraError("HMWSoln::initThermoXML",
+                                   "standardConc = molar_volume not done");
             } else if (formString == "solvent_volume") {
                 m_formGC = 2;
             } else {