diff --git a/include/cantera/equil/ChemEquil.h b/include/cantera/equil/ChemEquil.h
index 885f1e8b8..d9e1a4d31 100644
--- a/include/cantera/equil/ChemEquil.h
+++ b/include/cantera/equil/ChemEquil.h
@@ -48,10 +48,9 @@ public:
     int propertyPair;
 
     /**
-     * Continuation flag. Set true if the calculation should be
-     * initialized from the last calculation. Otherwise, the
-     * calculation will be started from scratch and the initial
-     * composition and element potentials estimated.
+     * Continuation flag. Set true if the calculation should be initialized from
+     * the last calculation. Otherwise, the calculation will be started from
+     * scratch and the initial composition and element potentials estimated.
      */
     bool contin;
 };
@@ -64,26 +63,22 @@ class PropertyCalculator;
  */
 
 /**
- *  Class ChemEquil implements a chemical equilibrium solver for
- *  single-phase solutions. It is a "non-stoichiometric" solver in
- *  the terminology of Smith and Missen, meaning that every
- *  intermediate state is a valid chemical equilibrium state, but
- *  does not necessarily satisfy the element constraints. In
- *  contrast, the solver implemented in class MultiPhaseEquil uses
- *  a "stoichiometric" algorithm, in which each intermediate state
- *  satisfies the element constraints but is not a state of
- *  chemical equilibrium. Non-stoichiometric methods are faster
- *  when they converge, but stoichiometric ones tend to be more
- *  robust and can be used also for problems with multiple
- *  condensed phases. As expected, the ChemEquil solver is faster
- *  than MultiPhaseEquil for many single-phase equilibrium
- *  problems (particularly if there are only a few elements but
- *  very many species), but can be less stable. Problem
- *  situations include low temperatures where only a few species
- *  have non-zero mole fractions, precisely stoichiometric
- *  compositions (e.g. 2 H2 + O2). In general, if speed is
- *  important, this solver should be tried first, and if it fails
- *  then use MultiPhaseEquil.
+ * Class ChemEquil implements a chemical equilibrium solver for single-phase
+ * solutions. It is a "non-stoichiometric" solver in the terminology of Smith
+ * and Missen, meaning that every intermediate state is a valid chemical
+ * equilibrium state, but does not necessarily satisfy the element constraints.
+ * In contrast, the solver implemented in class MultiPhaseEquil uses a
+ * "stoichiometric" algorithm, in which each intermediate state satisfies the
+ * element constraints but is not a state of chemical equilibrium. Non-
+ * stoichiometric methods are faster when they converge, but stoichiometric ones
+ * tend to be more robust and can be used also for problems with multiple
+ * condensed phases. As expected, the ChemEquil solver is faster than
+ * MultiPhaseEquil for many single-phase equilibrium problems (particularly if
+ * there are only a few elements but very many species), but can be less stable.
+ * Problem situations include low temperatures where only a few species have
+ * non-zero mole fractions, precisely stoichiometric compositions (e.g. 2 H2 +
+ * O2). In general, if speed is important, this solver should be tried first,
+ * and if it fails then use MultiPhaseEquil.
  * @ingroup equil
  */
 class ChemEquil
@@ -93,32 +88,33 @@ public:
 
     //! Constructor combined with the initialization function
     /*!
-     *  This constructor initializes the ChemEquil object with everything it
-     *  needs to start solving equilibrium problems.
-     *   @param s ThermoPhase object that will be used in the equilibrium calls.
+     * This constructor initializes the ChemEquil object with everything it
+     * needs to start solving equilibrium problems.
+     *
+     * @param s ThermoPhase object that will be used in the equilibrium calls.
      */
     ChemEquil(thermo_t& s);
 
     virtual ~ChemEquil();
 
     /*!
-     * Equilibrate a phase, holding the elemental composition fixed
-     * at the initial value found within the ThermoPhase object *s*.
+     * Equilibrate a phase, holding the elemental composition fixed at the
+     * initial value found within the ThermoPhase object *s*.
      *
-     * The value of 2 specified properties are obtained by querying the
-     * ThermoPhase object. The properties must be already contained
-     * within the current thermodynamic state of the system.
+     * The value of two specified properties are obtained by querying the
+     * ThermoPhase object. The properties must be already contained within the
+     * current thermodynamic state of the system.
      */
     int equilibrate(thermo_t& s, const char* XY,
                     bool useThermoPhaseElementPotentials = false, int loglevel = 0);
 
     /*!
-     * Compute the equilibrium composition for 2 specified
-     * properties and the specified element moles.
+     * Compute the equilibrium composition for two specified properties and the
+     * specified element moles.
      *
-     * The 2 specified properties are obtained by querying the
-     * ThermoPhase object. The properties must be already contained
-     * within the current thermodynamic state of the system.
+     * The two specified properties are obtained by querying the ThermoPhase
+     * object. The properties must be already contained within the current
+     * thermodynamic state of the system.
      *
      * @param s phase object to be equilibrated
      * @param XY property pair to hold constant
@@ -146,10 +142,9 @@ public:
 protected:
     //! Pointer to the ThermoPhase object used to initialize this object.
     /*!
-     *  This ThermoPhase object must be compatible with the ThermoPhase
-     *  objects input from the equilibrate function. Currently, this
-     *  means that the 2 ThermoPhases have to have consist of the same
-     *  species and elements.
+     * This ThermoPhase object must be compatible with the ThermoPhase objects
+     * input from the equilibrate function. Currently, this means that the 2
+     * ThermoPhases have to have consist of the same species and elements.
      */
     thermo_t* m_phase;
 
@@ -159,14 +154,14 @@ protected:
     }
 
     /*!
-     *  Prepare for equilibrium calculations.
-     *  @param s object representing the solution phase.
+     * Prepare for equilibrium calculations.
+     * @param s object representing the solution phase.
      */
     void initialize(thermo_t& s);
 
     /*!
-     * Set mixture to an equilibrium state consistent with specified
-     * element potentials and temperature.
+     * Set mixture to an equilibrium state consistent with specified element
+     * potentials and temperature.
      *
      * @param s mixture to be updated
      * @param x vector of non-dimensional element potentials
@@ -267,51 +262,36 @@ protected:
     size_t m_kk; //!< number of species in the phase
     size_t m_skip;
 
-    /**
-     * This is equal to the rank of the stoichiometric coefficient
-     * matrix when it is computed. It's initialized to #m_mm.
-     */
+    //! This is equal to the rank of the stoichiometric coefficient matrix when
+    //! it is computed. It's initialized to #m_mm.
     size_t m_nComponents;
 
     std::unique_ptr<PropertyCalculator<thermo_t> > m_p1, m_p2;
 
-    /**
-     * Current value of the mole fractions in the single phase.
-     * -> length = #m_kk.
-     */
+    //! Current value of the mole fractions in the single phase. length = #m_kk.
     vector_fp m_molefractions;
-    /**
-     * Current value of the dimensional element potentials
-     * -> length = #m_mm
-     */
+
+    //! Current value of the dimensional element potentials. length = #m_mm
     vector_fp m_lambda;
 
-    /*
-     * Current value of the sum of the element abundances given the
-     * current element potentials.
-     */
+    //! Current value of the sum of the element abundances given the current
+    //! element potentials.
     doublereal m_elementTotalSum;
-    /*
-     * Current value of the element mole fractions. Note these aren't
-     * the goal element mole fractions.
-     */
+
+    //! Current value of the element mole fractions. Note these aren't the goal
+    //! element mole fractions.
     vector_fp m_elementmolefracs;
     vector_fp m_reswork;
     vector_fp m_jwork1;
     vector_fp m_jwork2;
 
-    /*
-     * Storage of the element compositions
-     *      natom(k,m) = m_comp[k*m_mm+ m];
-     */
+    //! Storage of the element compositions. natom(k,m) = m_comp[k*m_mm+ m];
     vector_fp m_comp;
     doublereal m_temp, m_dens;
     doublereal m_p0;
 
-    /**
-     * Index of the element id corresponding to the electric charge of each
-     * species. Equal to -1 if there is no such element id.
-     */
+    //! Index of the element id corresponding to the electric charge of each
+    //! species. Equal to -1 if there is no such element id.
     size_t m_eloc;
 
     vector_fp m_startSoln;
@@ -319,11 +299,9 @@ protected:
     vector_fp m_grt;
     vector_fp m_mu_RT;
 
-    /**
-     * Dimensionless values of the Gibbs free energy for the
-     * standard state of each species, at the temperature and
-     * pressure of the solution (the star standard state).
-     */
+    //! Dimensionless values of the Gibbs free energy for the standard state of
+    //! each species, at the temperature and pressure of the solution (the star
+    //! standard state).
     vector_fp m_muSS_RT;
     std::vector<size_t> m_component;
 
diff --git a/include/cantera/equil/MultiPhase.h b/include/cantera/equil/MultiPhase.h
index 70a90da75..245e81f52 100644
--- a/include/cantera/equil/MultiPhase.h
+++ b/include/cantera/equil/MultiPhase.h
@@ -145,7 +145,7 @@ public:
     /*!
      * @param kGlob   global species index
      * @param mGlob   global element index
-     * @return        returns the number of atoms.
+     * @returns the number of atoms.
      */
     doublereal nAtoms(const size_t kGlob, const size_t mGlob) const;
 
@@ -176,9 +176,7 @@ public:
     //! Returns the index, given the phase name
     /*!
      * @param pName Name of the phase
-     *
-     * @return returns the index. A value of -1 means
-     *         the phase isn't in the object.
+     * @returns the index. A value of -1 means the phase isn't in the object.
      */
     int phaseIndex(const std::string& pName) const;
 
@@ -236,9 +234,10 @@ public:
      * @param speciesName    Species Name
      * @param phaseName      Phase Name
      *
-     * @return returns the global index
-     *  If the species or phase name is not recognized, this routine throws
-     *  a CanteraError.
+     * @returns the global index
+     *
+     * If the species or phase name is not recognized, this routine throws a
+     * CanteraError.
      */
     size_t speciesIndex(const std::string& speciesName, const std::string& phaseName);
 
@@ -273,14 +272,14 @@ public:
      */
     doublereal elementMoles(size_t m) const;
 
-    //!  Returns a vector of Chemical potentials.
+    //! Returns a vector of Chemical potentials.
     /*!
-     *  Write into array \a mu the chemical potentials of all species
-     *  [J/kmol]. The chemical potentials are related to the activities by
+     * Write into array \a mu the chemical potentials of all species
+     * [J/kmol]. The chemical potentials are related to the activities by
      *
-     *  \f$
+     * \f$
      *          \mu_k = \mu_k^0(T, P) + RT \ln a_k.
-     *  \f$.
+     * \f$.
      *
      * @param mu Chemical potential vector. Length = num global species. Units
      *           = J/kmol.
@@ -372,17 +371,17 @@ public:
 
     //! Set the state of the underlying ThermoPhase objects in one call
     /*!
-     *   @param T      Temperature of the system (kelvin)
-     *   @param Pres   pressure of the system (pascal)
+     * @param T      Temperature of the system (kelvin)
+     * @param Pres   pressure of the system (pascal)
      */
     void setState_TP(const doublereal T, const doublereal Pres);
 
     //! Set the state of the underlying ThermoPhase objects in one call
     /*!
-     *   @param T      Temperature of the system (kelvin)
-     *   @param Pres   pressure of the system (pascal)
-     *   @param Moles  Vector of mole numbers of all the species in all the phases
-     *                 (kmol)
+     * @param T      Temperature of the system (kelvin)
+     * @param Pres   pressure of the system (pascal)
+     * @param Moles  Vector of mole numbers of all the species in all the phases
+     *               (kmol)
      */
     void setState_TPMoles(const doublereal T, const doublereal Pres, const doublereal* Moles);
 
@@ -393,8 +392,8 @@ public:
 
     /// The total mixture volume [m^3].
     /*!
-     * Returns the cumulative sum of the volumes of all the
-     * phases in the mixture.
+     * Returns the cumulative sum of the volumes of all the phases in the
+     * mixture.
      */
     doublereal volume() const;
 
@@ -428,8 +427,8 @@ public:
         return m_np;
     }
 
-    //! Return true is species \a kGlob is a species in a
-    //! multicomponent solution phase.
+    //! Return true is species \a kGlob is a species in a multicomponent
+    //! solution phase.
     /*!
      * @param kGlob   index of the global species
      */
@@ -438,7 +437,7 @@ public:
     //! Returns the phase index of the Kth "global" species
     /*!
      * @param kGlob Global species index.
-     * @return Returns the index of the owning phase.
+     * @returns the index of the owning phase.
      */
     size_t speciesPhaseIndex(const size_t kGlob) const;
 
@@ -495,8 +494,8 @@ public:
 
     //! Adds moles of a certain species to the mixture
     /*!
-     *   @param indexS   Index of the species in the MultiPhase object
-     *   @param addedMoles   Value of the moles that are added to the species.
+     * @param indexS      Index of the species in the MultiPhase object
+     * @param addedMoles  Value of the moles that are added to the species.
      */
     void addSpeciesMoles(const int indexS, const doublereal addedMoles);
 
@@ -508,8 +507,8 @@ public:
      */
     void getElemAbundances(doublereal* elemAbundances) const;
 
-    //! Return true if the phase \a p has valid thermo data for
-    //! the current temperature.
+    //! Return true if the phase \a p has valid thermo data for the current
+    //! temperature.
     /*!
      * @param p  Index of the phase.
      */
@@ -517,31 +516,30 @@ public:
 
     // These methods are meant for internal use.
 
-    //! Update the locally-stored composition within this object
-    //! to match the current compositions of the phase objects.
+    //! Update the locally-stored composition within this object to match the
+    //! current compositions of the phase objects.
     /*!
-     *  Query the underlying ThermoPhase objects for their mole fractions and
-     *  fill in the mole fraction vector of this current object. Adjust
-     *  element compositions within this object to match.
+     * Query the underlying ThermoPhase objects for their mole fractions and
+     * fill in the mole fraction vector of this current object. Adjust element
+     * compositions within this object to match.
      *
-     *  This is an upload operation in the sense that we are taking downstream
-     *  information (ThermoPhase object info) and applying it to an upstream
-     *  object (MultiPhase object).
+     * This is an upload operation in the sense that we are taking downstream
+     * information (ThermoPhase object info) and applying it to an upstream
+     * object (MultiPhase object).
      */
     void uploadMoleFractionsFromPhases();
 
     //! Set the states of the phase objects to the locally-stored
     //! state within this MultiPhase object.
     /*!
-     *  This method sets each phase to the mixture temperature and pressure,
-     *  and sets the phase mole fractions based on the mixture mole numbers.
+     * This method sets each phase to the mixture temperature and pressure,
+     * and sets the phase mole fractions based on the mixture mole numbers.
      *
-     *  This is an download operation in the sense that we are taking
-     *  upstream object information (MultiPhase object) and
-     *  applying it to downstream objects (ThermoPhase object information)
+     * This is an download operation in the sense that we are taking upstream
+     * object information (MultiPhase object) and applying it to downstream
+     * objects (ThermoPhase object information)
      *
-     *  Therefore, the term, "update", is appropriate for a downstream
-     *  operation.
+     * Therefore, the term, "update", is appropriate for a downstream operation.
      */
     void updatePhases() const;
 
@@ -552,15 +550,14 @@ private:
     //! Set the mixture to a state of chemical equilibrium using the
     //! MultiPhaseEquil solver.
     /*!
-     *  @param XY   Integer flag specifying properties to hold fixed.
-     *  @param err  Error tolerance for \f$\Delta \mu/RT \f$ for all reactions.
-     *              Also used as the relative error tolerance for the outer
-     *              loop.
-     *  @param maxsteps Maximum number of steps to take in solving the fixed
-     *                  TP problem.
-     *  @param maxiter Maximum number of "outer" iterations for problems holding
-     *                 fixed something other than (T,P).
-     *  @param loglevel Level of diagnostic output
+     * @param XY   Integer flag specifying properties to hold fixed.
+     * @param err  Error tolerance for \f$\Delta \mu/RT \f$ for all reactions.
+     *             Also used as the relative error tolerance for the outer loop.
+     * @param maxsteps Maximum number of steps to take in solving the fixed TP
+     *                 problem.
+     * @param maxiter Maximum number of "outer" iterations for problems holding
+     *                fixed something other than (T,P).
+     * @param loglevel Level of diagnostic output
      */
     double equilibrate_MultiPhaseEquil(int XY, doublereal err, int maxsteps,
                                        int maxiter, int loglevel);
@@ -576,16 +573,14 @@ private:
 
     //! Global Stoichiometric Coefficient array
     /*!
-     *  This is a two dimensional array m_atoms(m, k). The first index is the
-     *  global element index. The second index, k, is the global species
-     *  index. The value is the number of atoms of type m in species k.
+     * This is a two dimensional array m_atoms(m, k). The first index is the
+     * global element index. The second index, k, is the global species index.
+     * The value is the number of atoms of type m in species k.
      */
     DenseMatrix m_atoms;
 
-    /**
-     * Locally stored vector of mole fractions of all species
-     * comprising the MultiPhase object.
-     */
+    //! Locally stored vector of mole fractions of all species comprising the
+    //! MultiPhase object.
     vector_fp m_moleFractions;
 
     //! Mapping between the global species number and the phase ID
@@ -603,23 +598,15 @@ private:
      */
     std::vector<size_t> m_spstart;
 
-    //! String names of the global elements
-    /*!
-     *  This has a length equal to the number of global elements.
-     */
+    //! String names of the global elements. This has a length equal to the
+    //! number of global elements.
     std::vector<std::string> m_enames;
 
-    //! Atomic number of each element
-    /*!
-     *  This is the atomic number of each global element.
-     */
+    //! Atomic number of each global element.
     vector_int m_atomicNumber;
 
-    //! Vector of species names in the problem
-    /*!
-     *   Vector is over all species defined in the object,
-     *   the global species index.
-     */
+    //! Vector of species names in the problem. Vector is over all species
+    //! defined in the object, the global species index.
     std::vector<std::string> m_snames;
 
     //! Returns the global element index, given the element string name
@@ -646,10 +633,8 @@ private:
     //! True if the init() routine has been called, and the MultiPhase frozen
     bool m_init;
 
-    //! Global ID of the element corresponding to the electronic charge.
-    /*!
-     * If there is none, then this is equal to -1
-     */
+    //! Global ID of the element corresponding to the electronic charge. If
+    //! there is none, then this is equal to -1
     size_t m_eloc;
 
     //! Vector of bools indicating whether temperatures are ok for phases.
@@ -659,17 +644,12 @@ private:
      */
     mutable std::vector<bool> m_temp_OK;
 
-    //! Minimum temperature for which thermo parameterizations are valid
-    /*!
-     *  Stoichiometric phases are ignored in this determination.
-     *  units Kelvin
-     */
+    //! Minimum temperature for which thermo parameterizations are valid.
+    //! Stoichiometric phases are ignored in this determination. units Kelvin
     doublereal m_Tmin;
 
-    //! Minimum temperature for which thermo parameterizations are valid
-    /*!
-     *  Stoichiometric phases are ignored in this determination. units Kelvin
-     */
+    //! Minimum temperature for which thermo parameterizations are valid.
+    //! Stoichiometric phases are ignored in this determination. units Kelvin
     doublereal m_Tmax;
 
     //! Vector of element abundances
@@ -711,9 +691,8 @@ inline std::ostream& operator<<(std::ostream& s, MultiPhase& x)
  * currently a linear combination of the previous components. Then, calculate
  * the stoichiometric coefficient matrix for that basis.
  *
- * Calculates the identity of the component species in the mechanism.
- * Rearranges the solution data to put the component data at the
- * front of the species list.
+ * Calculates the identity of the component species in the mechanism. Rearranges
+ * the solution data to put the component data at the front of the species list.
  *
  * Then, calculates SC(J,I) the formation reactions for all noncomponent
  * species in the mechanism.
diff --git a/src/equil/BasisOptimize.cpp b/src/equil/BasisOptimize.cpp
index f36b8b681..a4f200fd1 100644
--- a/src/equil/BasisOptimize.cpp
+++ b/src/equil/BasisOptimize.cpp
@@ -31,20 +31,16 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
     size_t j, jj, k=0, kk, l, i, jl, ml;
     std::string ename;
     std::string sname;
-    /*
-     * Get the total number of elements defined in the multiphase object
-     */
+
+    // Get the total number of elements defined in the multiphase object
     size_t ne = mphase->nElements();
-    /*
-     * Get the total number of species in the multiphase object
-     */
+
+    // Get the total number of species in the multiphase object
     size_t nspecies = mphase->nSpecies();
     doublereal tmp;
     doublereal const USEDBEFORE = -1;
 
-    /*
-     * Perhaps, initialize the element ordering
-     */
+    // Perhaps, initialize the element ordering
     if (orderVectorElements.size() < ne) {
         orderVectorElements.resize(ne);
         for (j = 0; j < ne; j++) {
@@ -52,9 +48,7 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
         }
     }
 
-    /*
-     * Perhaps, initialize the species ordering
-     */
+    // Perhaps, initialize the species ordering
     if (orderVectorSpecies.size() != nspecies) {
         orderVectorSpecies.resize(nspecies);
         for (k = 0; k < nspecies; k++) {
@@ -101,27 +95,20 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
         }
     }
 
-    /*
-     *  Calculate the maximum value of the number of components possible
-     *     It's equal to the minimum of the number of elements and the
-     *     number of total species.
-     */
+    // Calculate the maximum value of the number of components possible. It's
+    // equal to the minimum of the number of elements and the number of total
+    // species.
     size_t nComponents = std::min(ne, nspecies);
     size_t nNonComponents = nspecies - nComponents;
-    /*
-     * Set this return variable to false
-     */
+
+    // Set this return variable to false
     *usedZeroedSpecies = false;
 
-    /*
-     * Create an array of mole numbers
-     */
+    // Create an array of mole numbers
     vector_fp molNum(nspecies,0.0);
     mphase->getMoles(molNum.data());
 
-    /*
-     * Other workspace
-     */
+    // Other workspace
     vector_fp sm(ne*ne, 0.0);
     vector_fp ss(ne, 0.0);
     vector_fp sa(ne, 0.0);
@@ -129,9 +116,7 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
         formRxnMatrix.resize(nspecies*ne, 0.0);
     }
 
-    /*
-     * For debugging purposes keep an unmodified copy of the array.
-     */
+    // For debugging purposes keep an unmodified copy of the array.
     vector_fp molNumBase;
     if (DEBUG_MODE_ENABLED) {
         molNumBase = molNum;
@@ -139,19 +124,16 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
     double molSave = 0.0;
 
     size_t jr = 0;
-    /*
-     *   Top of a loop of some sort based on the index JR. JR is the
-     *   current number of component species found.
-     */
+
+    // Top of a loop of some sort based on the index JR. JR is the current
+    // number of component species found.
     while (jr < nComponents) {
         // Top of another loop point based on finding a linearly independent
         // species
         while (true) {
-            /*
-             *    Search the remaining part of the mole number vector, molNum
-             *    for the largest remaining species. Return its identity.
-             *    kk is the raw number. k is the orderVectorSpecies index.
-             */
+            // Search the remaining part of the mole number vector, molNum for
+            // the largest remaining species. Return its identity. kk is the raw
+            // number. k is the orderVectorSpecies index.
             kk = max_element(molNum.begin(), molNum.end()) - molNum.begin();
             for (j = 0; j < nspecies; j++) {
                 if (orderVectorSpecies[j] == kk) {
@@ -166,40 +148,33 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
             if (molNum[kk] == 0.0) {
                 *usedZeroedSpecies = true;
             }
-            /*
-             * If the largest molNum is negative, then we are done.
-             */
+            // If the largest molNum is negative, then we are done.
             if (molNum[kk] == USEDBEFORE) {
                 nComponents = jr;
                 nNonComponents = nspecies - nComponents;
                 break;
             }
-            /*
-             *  Assign a small negative number to the component that we have
-             *  just found, in order to take it out of further consideration.
-             */
+
+            // Assign a small negative number to the component that we have
+            // just found, in order to take it out of further consideration.
 #ifdef DEBUG_MODE
             molSave = molNum[kk];
 #endif
             molNum[kk] = USEDBEFORE;
 
-            // **** CHECK LINEAR INDEPENDENCE WITH PREVIOUS SPECIES ******
-            /*
-             *          Modified Gram-Schmidt Method, p. 202 Dalquist
-             *          QR factorization of a matrix without row pivoting.
-             */
+            // CHECK LINEAR INDEPENDENCE WITH PREVIOUS SPECIES
+
+            // Modified Gram-Schmidt Method, p. 202 Dalquist
+            // QR factorization of a matrix without row pivoting.
             jl = jr;
             for (j = 0; j < ne; ++j) {
                 jj = orderVectorElements[j];
                 sm[j + jr*ne] = mphase->nAtoms(kk,jj);
             }
             if (jl > 0) {
-                /*
-                 *         Compute the coefficients of JA column of the
-                 *         the upper triangular R matrix, SS(J) = R_J_JR
-                 *         (this is slightly different than Dalquist)
-                 *         R_JA_JA = 1
-                 */
+                // Compute the coefficients of JA column of the the upper
+                // triangular R matrix, SS(J) = R_J_JR (this is slightly
+                // different than Dalquist) R_JA_JA = 1
                 for (j = 0; j < jl; ++j) {
                     ss[j] = 0.0;
                     for (i = 0; i < ne; ++i) {
@@ -207,31 +182,31 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
                     }
                     ss[j] /= sa[j];
                 }
-                /*
-                 *     Now make the new column, (*,JR), orthogonal to the
-                 *     previous columns
-                 */
+
+                // Now make the new column, (*,JR), orthogonal to the previous
+                // columns
                 for (j = 0; j < jl; ++j) {
                     for (l = 0; l < ne; ++l) {
                         sm[l + jr*ne] -= ss[j] * sm[l + j*ne];
                     }
                 }
             }
-            /*
-             *        Find the new length of the new column in Q.
-             *        It will be used in the denominator in future row calcs.
-             */
+
+            // Find the new length of the new column in Q.
+            // It will be used in the denominator in future row calcs.
             sa[jr] = 0.0;
             for (ml = 0; ml < ne; ++ml) {
                 tmp = sm[ml + jr*ne];
                 sa[jr] += tmp * tmp;
             }
-            // **** IF NORM OF NEW ROW  .LT. 1E-3 REJECT **********
+
+            // IF NORM OF NEW ROW  .LT. 1E-3 REJECT
             if (sa[jr] > 1.0e-6) {
                 break;
             }
         }
-        // **** REARRANGE THE DATA ******************
+
+        // REARRANGE THE DATA
         if (jr != k) {
             if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl >= 1) {
                 kk = orderVectorSpecies[k];
@@ -252,38 +227,36 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
         return nComponents;
     }
 
-    // **** EVALUATE THE STOICHIOMETRY **********************
-    /*
-     *  Formulate the matrix problem for the stoichiometric
-     *  coefficients. CX + B = 0
-     *      C will be an nc x nc matrix made up of the formula
-     * vectors for the components. Each component's formula
-     * vector is a column. The rows are the elements.
-     *      n RHS's will be solved for. Thus, B is an nc x n
-     * matrix.
-     *
-     * BIG PROBLEM 1/21/99:
-     *
-     *    This algorithm makes the assumption that the
-     * first nc rows of the formula matrix aren't rank deficient.
-     * However, this might not be the case. For example, assume
-     * that the first element in FormulaMatrix[] is argon. Assume that
-     * no species in the matrix problem actually includes argon.
-     * Then, the first row in sm[], below will be identically
-     * zero. bleh.
-     *    What needs to be done is to perform a rearrangement
-     * of the ELEMENTS -> i.e. rearrange, FormulaMatrix, sp, and gai, such
-     * that the first nc elements form in combination with the
-     * nc components create an invertible sm[]. not a small
-     * project, but very doable.
-     *    An alternative would be to turn the matrix problem
-     * below into an ne x nc problem, and do QR elimination instead
-     * of Gauss-Jordan elimination.
-     *    Note the rearrangement of elements need only be done once
-     * in the problem. It's actually very similar to the top of
-     * this program with ne being the species and nc being the
-     * elements!!
-     */
+    // EVALUATE THE STOICHIOMETRY
+    //
+    // Formulate the matrix problem for the stoichiometric
+    // coefficients. CX + B = 0
+    //
+    // C will be an nc x nc matrix made up of the formula vectors for the
+    // components. Each component's formula vector is a column. The rows are the
+    // elements.
+    //
+    // n RHS's will be solved for. Thus, B is an nc x n matrix.
+    //
+    // BIG PROBLEM 1/21/99:
+    //
+    // This algorithm makes the assumption that the first nc rows of the formula
+    // matrix aren't rank deficient. However, this might not be the case. For
+    // example, assume that the first element in FormulaMatrix[] is argon.
+    // Assume that no species in the matrix problem actually includes argon.
+    // Then, the first row in sm[], below will be identically zero. bleh.
+    //
+    // What needs to be done is to perform a rearrangement of the ELEMENTS ->
+    // i.e. rearrange, FormulaMatrix, sp, and gai, such that the first nc
+    // elements form in combination with the nc components create an invertible
+    // sm[]. not a small project, but very doable.
+    //
+    // An alternative would be to turn the matrix problem below into an ne x nc
+    // problem, and do QR elimination instead of Gauss-Jordan elimination.
+    //
+    // Note the rearrangement of elements need only be done once in the problem.
+    // It's actually very similar to the top of this program with ne being the
+    // species and nc being the elements!!
     for (k = 0; k < nComponents; ++k) {
         kk = orderVectorSpecies[k];
         for (j = 0; j < nComponents; ++j) {
@@ -337,9 +310,7 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
             writelogf("   --- %3d (%3d) ", k, kk);
             writelogf("%-10.10s", mphase->speciesName(kk));
             writelogf("|%10.3g|", molNumBase[kk]);
-            /*
-             * Print the negative of formRxnMatrix[]; it's easier to interpret.
-             */
+            // Print the negative of formRxnMatrix[]; it's easier to interpret.
             for (j = 0; j < nComponents; j++) {
                 writelogf("     %6.2f", - formRxnMatrix[j + i * ne]);
             }
@@ -355,21 +326,17 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
     return nComponents;
 } // basopt()
 
+/**
+ * Print a string within a given space limit. This routine limits the amount of
+ * the string that will be printed to a maximum of "space" characters.
+ *
+ * str = String -> must be null terminated.
+ * space = space limit for the printing.
+ * alignment = 0 centered
+ *             1 right aligned
+ *             2 left aligned
+ */
 static void print_stringTrunc(const char* str, int space, int alignment)
-
-/***********************************************************************
- *  vcs_print_stringTrunc():
- *
- *     Print a string within a given space limit. This routine
- *     limits the amount of the string that will be printed to a
- *     maximum of "space" characters.
- *
- *     str = String -> must be null terminated.
- *     space = space limit for the printing.
- *     alignment = 0 centered
- *           1 right aligned
- *           2 left aligned
- ***********************************************************************/
 {
     int i, ls=0, rs=0;
     int len = static_cast<int>(strlen(str));
@@ -408,9 +375,7 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
     size_t j, k, l, i, jl, ml, jr, ielem, jj, kk=0;
     size_t nelements = mphase->nElements();
     std::string ename;
-    /*
-     * Get the total number of species in the multiphase object
-     */
+    // Get the total number of species in the multiphase object
     size_t nspecies = mphase->nSpecies();
 
     double test = -1.0E10;
@@ -425,9 +390,7 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
         writelog("   ---    and to rearrange the element ordering once\n");
     }
 
-    /*
-     * Perhaps, initialize the element ordering
-     */
+    // Perhaps, initialize the element ordering
     if (orderVectorElements.size() < nelements) {
         orderVectorElements.resize(nelements);
         for (j = 0; j < nelements; j++) {
@@ -435,11 +398,8 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
         }
     }
 
-    /*
-     * Perhaps, initialize the species ordering. However, this is
-     * dangerous, as this ordering is assumed to yield the
-     * component species for the problem
-     */
+    // Perhaps, initialize the species ordering. However, this is dangerous, as
+    // this ordering is assumed to yield the component species for the problem
     if (orderVectorSpecies.size() != nspecies) {
         orderVectorSpecies.resize(nspecies);
         for (k = 0; k < nspecies; k++) {
@@ -447,12 +407,9 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
         }
     }
 
-    /*
-     * If the elementAbundances aren't input, just create a fake one
-     * based on summing the column of the stoich matrix.
-     * This will force elements with zero species to the
-     * end of the element ordering.
-     */
+    // If the elementAbundances aren't input, just create a fake one based on
+    // summing the column of the stoich matrix. This will force elements with
+    // zero species to the end of the element ordering.
     vector_fp eAbund(nelements,0.0);
     if (elementAbundances.size() != nelements) {
         for (j = 0; j < nelements; j++) {
@@ -470,23 +427,16 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
     vector_fp ss(nelements,0.0);
     vector_fp sm(nelements*nelements,0.0);
 
-    /*
-     *        Top of a loop of some sort based on the index JR. JR is the
-     *       current number independent elements found.
-     */
+    // Top of a loop of some sort based on the index JR. JR is the current
+    // number independent elements found.
     jr = 0;
     while (jr < nComponents) {
-        /*
-         *     Top of another loop point based on finding a linearly
-         *     independent element
-         */
+        // Top of another loop point based on finding a linearly independent
+        // element
         while (true) {
-            /*
-             *    Search the element vector. We first locate elements that
-             *    are present in any amount. Then, we locate elements that
-             *    are not present in any amount.
-             *    Return its identity in K.
-             */
+            // Search the element vector. We first locate elements that are
+            // present in any amount. Then, we locate elements that are not
+            // present in any amount. Return its identity in K.
             k = nelements;
             for (ielem = jr; ielem < nelements; ielem++) {
                 kk = orderVectorElements[ielem];
@@ -512,37 +462,30 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
                 throw CanteraError("ElemRearrange", "Required number of elements not found.");
             }
 
-            /*
-             *  Assign a large negative number to the element that we have
-             *  just found, in order to take it out of further consideration.
-             */
+            // Assign a large negative number to the element that we have
+            // just found, in order to take it out of further consideration.
             eAbund[kk] = test;
 
-            // **** CHECK LINEAR INDEPENDENCE OF CURRENT FORMULA MATRIX
-            // **** LINE WITH PREVIOUS LINES OF THE FORMULA MATRIX  ******
-            /*
-             *          Modified Gram-Schmidt Method, p. 202 Dalquist
-             *          QR factorization of a matrix without row pivoting.
-             */
+            // CHECK LINEAR INDEPENDENCE OF CURRENT FORMULA MATRIX
+            // LINE WITH PREVIOUS LINES OF THE FORMULA MATRIX
+
+            // Modified Gram-Schmidt Method, p. 202 Dalquist
+            // QR factorization of a matrix without row pivoting.
             jl = jr;
-            /*
-             *   Fill in the row for the current element, k, under consideration
-             *   The row will contain the Formula matrix value for that element
-             *   with respect to the vector of component species.
-             *   (note j and k indices are flipped compared to the previous routine)
-             */
+
+            // Fill in the row for the current element, k, under consideration
+            // The row will contain the Formula matrix value for that element
+            // with respect to the vector of component species. (note j and k
+            // indices are flipped compared to the previous routine)
             for (j = 0; j < nComponents; ++j) {
                 jj = orderVectorSpecies[j];
                 kk = orderVectorElements[k];
                 sm[j + jr*nComponents] = mphase->nAtoms(jj,kk);
             }
             if (jl > 0) {
-                /*
-                 *         Compute the coefficients of JA column of the
-                 *         the upper triangular R matrix, SS(J) = R_J_JR
-                 *         (this is slightly different than Dalquist)
-                 *         R_JA_JA = 1
-                 */
+                // Compute the coefficients of JA column of the the upper
+                // triangular R matrix, SS(J) = R_J_JR (this is slightly
+                // different than Dalquist) R_JA_JA = 1
                 for (j = 0; j < jl; ++j) {
                     ss[j] = 0.0;
                     for (i = 0; i < nComponents; ++i) {
@@ -550,10 +493,9 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
                     }
                     ss[j] /= sa[j];
                 }
-                /*
-                 *     Now make the new column, (*,JR), orthogonal to the
-                 *     previous columns
-                 */
+
+                // Now make the new column, (*,JR), orthogonal to the
+                // previous columns
                 for (j = 0; j < jl; ++j) {
                     for (l = 0; l < nComponents; ++l) {
                         sm[l + jr*nComponents] -= ss[j] * sm[l + j*nComponents];
@@ -561,21 +503,19 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
                 }
             }
 
-            /*
-             *        Find the new length of the new column in Q.
-             *        It will be used in the denominator in future row calcs.
-             */
+            // Find the new length of the new column in Q.
+            // It will be used in the denominator in future row calcs.
             sa[jr] = 0.0;
             for (ml = 0; ml < nComponents; ++ml) {
                 double tmp = sm[ml + jr*nComponents];
                 sa[jr] += tmp * tmp;
             }
-            // **** IF NORM OF NEW ROW  .LT. 1E-6 REJECT **********
+            // IF NORM OF NEW ROW  .LT. 1E-6 REJECT
             if (sa[jr] > 1.0e-6) {
                 break;
             }
         }
-        // **** REARRANGE THE DATA ******************
+        // REARRANGE THE DATA
         if (jr != k) {
             if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl > 0) {
                 kk = orderVectorElements[k];
diff --git a/src/equil/ChemEquil.cpp b/src/equil/ChemEquil.cpp
index 62ae84e58..444cb2f69 100644
--- a/src/equil/ChemEquil.cpp
+++ b/src/equil/ChemEquil.cpp
@@ -201,11 +201,9 @@ int ChemEquil::setInitialMoles(thermo_t& s, vector_fp& elMoleGoal,
         for (size_t m = 0; m < m_nComponents; m++) {
             m_component[m] = e.componentIndex(m);
         }
-        /*
-         * Update the current values of the temp, density, and
-         * mole fraction, and element abundance vectors kept
-         * within the ChemEquil object.
-         */
+
+        // Update the current values of the temp, density, and mole fraction,
+        // and element abundance vectors kept within the ChemEquil object.
         update(s);
 
         if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) {
@@ -347,9 +345,7 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
     vector_fp state;
     s.saveState(state);
 
-    /*
-     * Check Compatibility
-     */
+    // Check Compatibility
     if (m_mm != s.nElements() || m_kk != s.nSpecies()) {
         throw CanteraError("ChemEquil::equilibrate ERROR",
                            "Input ThermoPhase is incompatible with initialization");
@@ -406,11 +402,8 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
         }
     }
 
-    /*
-     * Before we do anything to change the ThermoPhase object,
-     * we calculate and store the two specified thermodynamic
-     * properties that we are after.
-     */
+    // Before we do anything to change the ThermoPhase object, we calculate and
+    // store the two specified thermodynamic properties that we are after.
     xval = m_p1->value(s);
     yval = m_p2->value(s);
 
@@ -420,13 +413,10 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
     vector_fp x(nvar, -102.0); // solution vector
     vector_fp res_trial(nvar, 0.0); // residual
 
-    /*
-     * Replace one of the element abundance fraction equations
-     * with the specified property calculation.
-     *
-     * We choose the equation of the element with the highest element
-     * abundance.
-     */
+    // Replace one of the element abundance fraction equations with the
+    // specified property calculation.
+    //
+    // We choose the equation of the element with the highest element abundance.
     size_t m;
     tmp = -1.0;
     for (size_t im = 0; im < m_nComponents; im++) {
@@ -441,20 +431,17 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
                            "Element Abundance Vector is zeroed");
     }
 
-    // start with a composition with everything non-zero. Note
-    // that since we have already save the target element moles,
-    // changing the composition at this point only affects the
-    // starting point, not the final solution.
+    // start with a composition with everything non-zero. Note that since we
+    // have already save the target element moles, changing the composition at
+    // this point only affects the starting point, not the final solution.
     vector_fp xmm(m_kk, 0.0);
     for (size_t k = 0; k < m_kk; k++) {
         xmm[k] = s.moleFraction(k) + 1.0E-32;
     }
     s.setMoleFractions(xmm.data());
 
-    /*
-     * Update the internally stored values of m_temp,
-     * m_dens, and the element mole fractions.
-     */
+    // Update the internally stored values of m_temp, m_dens, and the element
+    // mole fractions.
     update(s);
 
     doublereal tmaxPhase = s.maxTemp();
@@ -472,10 +459,9 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
 
         doublereal slope, phigh, plow, pval, dt;
 
-        // first get the property values at the upper and lower
-        // temperature limits. Since p1 (h, s, or u) is monotonic
-        // in T, these values determine the upper and lower
-        // bounnds (phigh, plow) for p1.
+        // first get the property values at the upper and lower temperature
+        // limits. Since p1 (h, s, or u) is monotonic in T, these values
+        // determine the upper and lower bounnds (phigh, plow) for p1.
 
         s.setTemperature(tmax);
         setInitialMoles(s, elMolesGoal, loglevel - 1);
@@ -495,10 +481,10 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
             setInitialMoles(s, elMolesGoal, loglevel - 1);
             pval = m_p1->value(s);
 
-            // If this value of p1 is greater than the specified
-            // property value, then the current temperature is too
-            // high. Use it as the new upper bound. Otherwise, it
-            // is too low, so use it as the new lower bound.
+            // If this value of p1 is greater than the specified property value,
+            // then the current temperature is too high. Use it as the new upper
+            // bound. Otherwise, it is too low, so use it as the new lower
+            // bound.
             if (pval > xval) {
                 tmax = t0;
                 phigh = pval;
@@ -534,11 +520,8 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
 
     setInitialMoles(s, elMolesGoal,loglevel);
 
-    /*
-     * If requested, get the initial estimate for the
-     * chemical potentials from the ThermoPhase object
-     * itself. Or else, create our own estimate.
-     */
+    // If requested, get the initial estimate for the chemical potentials from
+    // the ThermoPhase object itself. Or else, create our own estimate.
     if (useThermoPhaseElementPotentials) {
         bool haveEm = s.getElementPotentials(x.data());
         if (haveEm) {
@@ -552,44 +535,32 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
             estimateElementPotentials(s, x, elMolesGoal);
         }
     } else {
-        /*
-         * Calculate initial estimates of the element potentials.
-         * This algorithm uese the MultiPhaseEquil object's
-         * initialization capabilities to calculate an initial
-         * estimate of the mole fractions for a set of linearly
-         * independent component species. Then, the element
-         * potentials are solved for based on the chemical
-         * potentials of the component species.
-         */
+        // Calculate initial estimates of the element potentials. This algorithm
+        // uese the MultiPhaseEquil object's initialization capabilities to
+        // calculate an initial estimate of the mole fractions for a set of
+        // linearly independent component species. Then, the element potentials
+        // are solved for based on the chemical potentials of the component
+        // species.
         estimateElementPotentials(s, x, elMolesGoal);
     }
 
-    /*
-     * Do a better estimate of the element potentials.
-     * We have found that the current estimate may not be good
-     * enough to avoid drastic numerical issues associated with
-     * the use of a numerically generated Jacobian.
-     *
-     * The Brinkley algorithm assumes a constant T, P system
-     * and uses a linearized analytical Jacobian that turns out
-     * to be very stable.
-     */
+    // Do a better estimate of the element potentials. We have found that the
+    // current estimate may not be good enough to avoid drastic numerical issues
+    // associated with the use of a numerically generated Jacobian.
+    //
+    // The Brinkley algorithm assumes a constant T, P system and uses a
+    // linearized analytical Jacobian that turns out to be very stable.
     int info = estimateEP_Brinkley(s, x, elMolesGoal);
     if (info == 0) {
         setToEquilState(s, x, s.temperature());
     }
 
-    /*
-     * Install the log(temp) into the last solution unknown
-     * slot.
-     */
+    // Install the log(temp) into the last solution unknown slot.
     x[m_mm] = log(s.temperature());
 
-    /*
-     * Setting the max and min values for x[]. Also, if element
-     * abundance vector is zero, setting x[] to -1000. This
-     * effectively zeroes out all species containing that element.
-     */
+    // Setting the max and min values for x[]. Also, if element abundance vector
+    // is zero, setting x[] to -1000. This effectively zeroes out all species
+    // containing that element.
     vector_fp above(nvar);
     vector_fp below(nvar);
     for (m = 0; m < mm; m++) {
@@ -599,10 +570,9 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
             x[m] = -1000.0;
         }
     }
-    /*
-     * Set the temperature bounds to be 25 degrees different than the max and min
-     * temperatures.
-     */
+
+    // Set the temperature bounds to be 25 degrees different than the max and
+    // min temperatures.
     above[mm] = log(s.maxTemp() + 25.0);
     below[mm] = log(s.minTemp() - 25.0);
 
@@ -625,19 +595,17 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
         for (m = 0; m < nvar; m++) {
             double tval = options.relTolerance;
             if (m < mm) {
-                /*
-                 * Special case convergence requirements for electron element.
-                 * This is a special case because the element coefficients may
-                 * be both positive and negative. And, typically they sum to 0.0.
-                 * Therefore, there is no natural absolute value for this quantity.
-                 * We supply the absolute value tolerance here. Note, this is
-                 * made easier since the element abundances are normalized to one
-                 * within this routine.
-                 *
-                 * Note, the 1.0E-13 value was recently relaxed from 1.0E-15, because
-                 * convergence failures were found to occur for the lower value
-                 * at small pressure (0.01 pascal).
-                 */
+                // Special case convergence requirements for electron element.
+                // This is a special case because the element coefficients may
+                // be both positive and negative. And, typically they sum to
+                // 0.0. Therefore, there is no natural absolute value for this
+                // quantity. We supply the absolute value tolerance here. Note,
+                // this is made easier since the element abundances are
+                // normalized to one within this routine.
+                //
+                // Note, the 1.0E-13 value was recently relaxed from 1.0E-15,
+                // because convergence failures were found to occur for the
+                // lower value at small pressure (0.01 pascal).
                 if (m == m_eloc) {
                     tval = elMolesGoal[m] * options.relTolerance + options.absElemTol
                            + 1.0E-13;
@@ -659,10 +627,9 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
             if (m_eloc != npos) {
                 adjustEloc(s, elMolesGoal);
             }
-            /*
-             * Save the calculated and converged element potentials
-             * to the original ThermoPhase object.
-             */
+
+            // Save the calculated and converged element potentials to the
+            // original ThermoPhase object.
             s.setElementPotentials(m_lambda);
             if (s.temperature() > s.maxTemp() + 1.0 ||
                     s.temperature() < s.minTemp() - 1.0) {
@@ -705,9 +672,7 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
         oldf = f;
         scale(res_trial.begin(), res_trial.end(), res_trial.begin(), -1.0);
 
-        /*
-         * Solve the system
-         */
+        // Solve the system
         try {
             info = solve(jac, res_trial.data());
         } catch (CanteraError& err) {
@@ -776,9 +741,8 @@ int ChemEquil::dampStep(thermo_t& mix, vector_fp& oldx,
 {
     double damp;
 
-    /*
-     * Carry out a delta damping approach on the dimensionless element potentials.
-     */
+    // Carry out a delta damping approach on the dimensionless element
+    // potentials.
     damp = 1.0;
     for (size_t m = 0; m < m_mm; m++) {
         if (m == m_eloc) {
@@ -798,9 +762,7 @@ int ChemEquil::dampStep(thermo_t& mix, vector_fp& oldx,
         }
     }
 
-    /*
-     * Update the solution unknown
-     */
+    // Update the solution unknown
     for (size_t m = 0; m < x.size(); m++) {
         x[m] = oldx[m] + damp * step[m];
     }
@@ -832,11 +794,8 @@ void ChemEquil::equilResidual(thermo_t& s, const vector_fp& x,
         } else if (n >= m_nComponents) {
             resid[m] = x[m];
         } else {
-            /*
-             * Change the calculation for small element number, using
-             * L'Hopital's rule.
-             * The log formulation is unstable.
-             */
+            // Change the calculation for small element number, using
+            // L'Hopital's rule. The log formulation is unstable.
             if (elmFracGoal[m] < 1.0E-10 || elmFrac[m] < 1.0E-10 || m == m_eloc) {
                 resid[m] = elmFracGoal[m] - elmFrac[m];
             } else {
@@ -907,10 +866,9 @@ double ChemEquil::calcEmoles(thermo_t& s, vector_fp& x, const double& n_t,
 {
     double n_t_calc = 0.0;
     double tmp;
-    /*
-     * Calculate the activity coefficients of the solution, at the
-     * previous solution state.
-     */
+
+    // Calculate the activity coefficients of the solution, at the previous
+    // solution state.
     vector_fp actCoeff(m_kk, 1.0);
     s.setMoleFractions(Xmol_i_calc.data());
     s.setPressure(pressureConst);
@@ -941,11 +899,8 @@ double ChemEquil::calcEmoles(thermo_t& s, vector_fp& x, const double& n_t,
 int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
                                    vector_fp& elMoles)
 {
-    /*
-     * Before we do anything, we will save the state of the solution.
-     * Then, if things go drastically wrong, we will restore the
-     * saved state.
-     */
+    // Before we do anything, we will save the state of the solution. Then, if
+    // things go drastically wrong, we will restore the saved state.
     vector_fp state;
     s.saveState(state);
     double tmp, sum;
@@ -969,11 +924,8 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
     vector_fp resid(m_mm+1, 0.0);
     vector_int lumpSum(m_mm+1, 0);
 
-    /*
-     * Get the nondimensional Gibbs functions for the species
-     * at their standard states of solution at the current T and P
-     * of the solution.
-     */
+    // Get the nondimensional Gibbs functions for the species at their standard
+    // states of solution at the current T and P of the solution.
     s.getGibbs_RT(m_muSS_RT.data());
 
     vector_fp eMolesCalc(m_mm, 0.0);
@@ -1042,27 +994,18 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
         }
     }
 
-    /*
-     * -------------------------------------------------------------------
-     * Main Loop.
-     */
+    // Main Loop.
     for (int iter = 0; iter < 20* options.maxIterations; iter++) {
-        /*
-         * Save the old solution
-         */
+        // Save the old solution
         for (m = 0; m < m_mm; m++) {
             x_old[m] = x[m];
         }
         x_old[m_mm] = n_t;
-        /*
-         * Calculate the mole numbers of species
-         */
+        // Calculate the mole numbers of species
         if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) {
             writelogf("START ITERATION %d:\n", iter);
         }
-        /*
-         * Calculate the mole numbers of species and elements.
-         */
+        // Calculate the mole numbers of species and elements.
         double n_t_calc = calcEmoles(s, x, n_t, Xmol_i_calc, eMolesCalc, n_i_calc,
                                      pressureConst);
 
@@ -1086,9 +1029,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
 
         double nCutoff;
         bool normalStep = true;
-        /*
-         * Decide if we are to do a normal step or a modified step
-         */
+        // Decide if we are to do a normal step or a modified step
         size_t iM = npos;
         for (m = 0; m < m_mm; m++) {
             if (elMoles[m] > 0.001 * elMolesTotal) {
@@ -1130,30 +1071,28 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
                 resid[m_mm] = std::min(resid[m_mm], 0.0);
             }
         } else {
-            /*
-             * Determine whether the matrix should be dumbed down because
-             * the coefficient matrix of species (with significant concentrations)
-             * is rank deficient.
-             *
-             * The basic idea is that at any time during the calculation only a
-             * small subset of species with sufficient concentration matters.
-             * If the rank of the element coefficient matrix for that subset of species
-             * is less than the number of elements, then the matrix created by
-             * the Brinkley method below may become singular.
-             *
-             * The logic below looks for obvious cases where the current element
-             * coefficient matrix is rank deficient.
-             *
-             * The way around rank-deficiency is to lump-sum the corresponding row
-             * of the matrix. Note, lump-summing seems to work very well in terms of
-             * its stability properties, i.e., it heads in the right direction,
-             * albeit with lousy convergence rates.
-             *
-             * NOTE: This probably should be extended to a full blown Gauss-Jordan
-             *       factorization scheme in the future. For Example
-             *       the scheme below would fail for the set: HCl  NH4Cl, NH3.
-             *       Hopefully, it's caught by the equal rows logic below.
-             */
+            // Determine whether the matrix should be dumbed down because the
+            // coefficient matrix of species (with significant concentrations)
+            // is rank deficient.
+            //
+            // The basic idea is that at any time during the calculation only a
+            // small subset of species with sufficient concentration matters. If
+            // the rank of the element coefficient matrix for that subset of
+            // species is less than the number of elements, then the matrix
+            // created by the Brinkley method below may become singular.
+            //
+            // The logic below looks for obvious cases where the current element
+            // coefficient matrix is rank deficient.
+            //
+            // The way around rank-deficiency is to lump-sum the corresponding
+            // row of the matrix. Note, lump-summing seems to work very well in
+            // terms of its stability properties, i.e., it heads in the right
+            // direction, albeit with lousy convergence rates.
+            //
+            // NOTE: This probably should be extended to a full blown Gauss-
+            // Jordan factorization scheme in the future. For Example the scheme
+            // below would fail for the set: HCl  NH4Cl, NH3. Hopefully, it's
+            // caught by the equal rows logic below.
             for (m = 0; m < m_mm; m++) {
                 lumpSum[m] = 1;
             }
@@ -1189,9 +1128,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
                 }
             }
 
-            /*
-             * Formulate the matrix.
-             */
+            // Formulate the matrix.
             for (im = 0; im < m_mm; im++) {
                 m = m_orderVectorElements[im];
                 if (im < m_nComponents) {
@@ -1212,9 +1149,8 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
             }
             a1(m_mm, m_mm) = 0.0;
 
-            /*
-             * Formulate the residual, resid, and the estimate for the convergence criteria, sum
-             */
+            // Formulate the residual, resid, and the estimate for the
+            // convergence criteria, sum
             sum = 0.0;
             for (im = 0; im < m_mm; im++) {
                 m = m_orderVectorElements[im];
@@ -1223,12 +1159,12 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
                 } else {
                     resid[m] = 0.0;
                 }
-                /*
-                 * For equations with positive and negative coefficients, (electronic charge),
-                 * we must mitigate the convergence criteria by a condition limited by
-                 * finite precision of inverting a matrix.
-                 * Other equations with just positive coefficients aren't limited by this.
-                 */
+
+                // For equations with positive and negative coefficients,
+                // (electronic charge), we must mitigate the convergence
+                // criteria by a condition limited by finite precision of
+                // inverting a matrix. Other equations with just positive
+                // coefficients aren't limited by this.
                 if (m == m_eloc) {
                     tmp = resid[m] / (elMoles[m] + elMolesTotal*1.0E-6 + options.absElemTol);
                 } else {
@@ -1274,21 +1210,17 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
             if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) {
                 writelogf("(it %d) Convergence = %g\n", iter, sum);
             }
-            /*
-             * Insist on 20x accuracy compared to the top routine.
-             * There are instances, for ill-conditioned or
-             * singular matrices where this is needed to move
-             * the system to a point where the matrices aren't
-             * singular.
-             */
+
+            // Insist on 20x accuracy compared to the top routine. There are
+            // instances, for ill-conditioned or singular matrices where this is
+            // needed to move the system to a point where the matrices aren't
+            // singular.
             if (sum < 0.05 * options.relTolerance) {
                 retn = 0;
                 break;
             }
 
-            /*
-             * Row Sum scaling
-             */
+            // Row Sum scaling
             for (m = 0; m <= m_mm; m++) {
                 tmp = 0.0;
                 for (n = 0; n <= m_mm; n++) {
@@ -1322,27 +1254,26 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
                     writelogf("]  =   %10.5g\n", resid[m]);
                 }
             }
-            /*
-             * Next Step: We have row-summed the equations.
-             * However, there are some degenerate cases where two
-             * rows will be multiplies of each other in terms of
-             * 0 < m, 0 < m part of the matrix. This occurs on a case
-             * by case basis, and depends upon the current state of the
-             * element potential values, which affect the concentrations
-             * of species.
-             * So, the way we have found to eliminate this problem is to
-             * lump-sum one of the rows of the matrix, except for the
-             * last column, and stick it all on the diagonal.
-             * Then, we at least have a non-singular matrix, and the
-             * modified equation moves the corresponding unknown in the
-             * correct direction.
-             * The previous row-sum operation has made the identification
-             * of identical rows much simpler.
-             *
-             * Note at least 6E-4 is necessary for the comparison.
-             * I'm guessing 1.0E-3. If two rows are anywhere close to being
-             * equivalent, the algorithm can get stuck in an oscillatory mode.
-             */
+
+            // Next Step: We have row-summed the equations. However, there are
+            // some degenerate cases where two rows will be multiplies of each
+            // other in terms of 0 < m, 0 < m part of the matrix. This occurs on
+            // a case by case basis, and depends upon the current state of the
+            // element potential values, which affect the concentrations of
+            // species.
+            //
+            // So, the way we have found to eliminate this problem is to lump-
+            // sum one of the rows of the matrix, except for the last column,
+            // and stick it all on the diagonal. Then, we at least have a non-
+            // singular matrix, and the modified equation moves the
+            // corresponding unknown in the correct direction.
+            //
+            // The previous row-sum operation has made the identification of
+            // identical rows much simpler.
+            //
+            // Note at least 6E-4 is necessary for the comparison. I'm guessing
+            // 1.0E-3. If two rows are anywhere close to being equivalent, the
+            // algorithm can get stuck in an oscillatory mode.
             modifiedMatrix = false;
             for (m = 0; m < m_mm; m++) {
                 size_t sameAsRow = npos;
@@ -1401,11 +1332,8 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
                                    "unused elements.");
             }
 
-            /*
-             * Figure out the damping coefficient: Use a delta damping
-             * coefficient formulation: magnitude of change is capped
-             * to exp(1).
-             */
+            // Figure out the damping coefficient: Use a delta damping
+            // coefficient formulation: magnitude of change is capped to exp(1).
             beta = 1.0;
             for (m = 0; m < m_mm; m++) {
                 if (resid[m] > 1.0) {
@@ -1419,9 +1347,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
                 writelogf("(it %d) Beta = %g\n", iter, beta);
             }
         }
-        /*
-         * Update the solution vector
-         */
+        // Update the solution vector
         for (m = 0; m < m_mm; m++) {
             x[m] += beta * resid[m];
         }
diff --git a/src/equil/MultiPhase.cpp b/src/equil/MultiPhase.cpp
index 49e9e6f67..f543e6338 100644
--- a/src/equil/MultiPhase.cpp
+++ b/src/equil/MultiPhase.cpp
@@ -107,9 +107,8 @@ void MultiPhase::addPhase(ThermoPhase* p, doublereal moles)
     m_np = m_phase.size();
     m_nsp += p->nSpecies();
 
-    // determine if this phase has new elements
-    // for each new element, add an entry in the map
-    // from names to index number + 1:
+    // determine if this phase has new elements for each new element, add an
+    // entry in the map from names to index number + 1:
 
     string ename;
     // iterate over the elements in this phase
@@ -117,10 +116,9 @@ void MultiPhase::addPhase(ThermoPhase* p, doublereal moles)
     for (m = 0; m < nel; m++) {
         ename = p->elementName(m);
 
-        // if no entry is found for this element name, then
-        // it is a new element. In this case, add the name
-        // to the list of names, increment the element count,
-        // and add an entry to the name->(index+1) map.
+        // if no entry is found for this element name, then it is a new element.
+        // In this case, add the name to the list of names, increment the
+        // element count, and add an entry to the name->(index+1) map.
         if (m_enamemap.find(ename) == m_enamemap.end()) {
             m_enamemap[ename] = m_nel + 1;
             m_enames.push_back(ename);
@@ -135,22 +133,21 @@ void MultiPhase::addPhase(ThermoPhase* p, doublereal moles)
         }
     }
 
-    // If the mixture temperature hasn't been set, then set the
-    // temperature and pressure to the values for the phase being
-    // added. There is no good way to do this. However, this will be overridden later.
+    // If the mixture temperature hasn't been set, then set the temperature and
+    // pressure to the values for the phase being added. There is no good way to
+    // do this. However, this will be overridden later.
     if (m_temp == 298.15 && p->temperature() > 2.0E-3) {
         m_temp = p->temperature();
         m_press = p->pressure();
     }
 
-    // If this is a solution phase, update the minimum and maximum
-    // mixture temperatures. Stoichiometric phases are excluded,
-    // since a mixture may define multiple stoichiometric phases,
-    // each of which has thermo data valid only over a limited
-    // range. For example, a mixture might be defined to contain a
-    // phase representing water ice and one representing liquid
-    // water, only one of which should be present if the mixture
-    // represents an equilibrium state.
+    // If this is a solution phase, update the minimum and maximum mixture
+    // temperatures. Stoichiometric phases are excluded, since a mixture may
+    // define multiple stoichiometric phases, each of which has thermo data
+    // valid only over a limited range. For example, a mixture might be defined
+    // to contain a phase representing water ice and one representing liquid
+    // water, only one of which should be present if the mixture represents an
+    // equilibrium state.
     if (p->nSpecies() > 1) {
         m_Tmin = std::max(p->minTemp(), m_Tmin);
         m_Tmax = std::min(p->maxTemp(), m_Tmax);
@@ -172,8 +169,7 @@ void MultiPhase::init()
     m_elemAbundances.resize(m_nel, 0.0);
 
     // iterate over the elements
-    //   -> fill in m_atoms(m,k), m_snames(k), m_spphase(k),
-    //              m_sptart(ip)
+    //   -> fill in m_atoms(m,k), m_snames(k), m_spphase(k), m_spstart(ip)
     for (m = 0; m < m_nel; m++) {
         sym = m_enames[m];
         k = 0;
@@ -210,8 +206,8 @@ void MultiPhase::init()
         }
     }
 
-    /// set the initial composition within each phase to the
-    /// mole fractions stored in the phase objects
+    // set the initial composition within each phase to the
+    // mole fractions stored in the phase objects
     m_init = true;
     uploadMoleFractionsFromPhases();
     updatePhases();
@@ -437,9 +433,7 @@ void MultiPhase::setMolesByName(const std::string& x)
 
 void MultiPhase::getMoles(doublereal* molNum) const
 {
-    /*
-     * First copy in the mole fractions
-     */
+    // First copy in the mole fractions
     copy(m_moleFractions.begin(), m_moleFractions.end(), molNum);
     size_t ik;
     doublereal* dtmp = molNum;
@@ -605,8 +599,8 @@ double MultiPhase::equilibrate_MultiPhaseEquil(int XY, doublereal err,
                         Hlow = hnow;
                     }
                 } else {
-                    // the current enthalpy is greater than the target; therefore the
-                    // current temperature is too high.
+                    // the current enthalpy is greater than the target;
+                    // therefore the current temperature is too high.
                     if (m_temp < Thigh) {
                         Thigh = m_temp;
                         Hhigh = hnow;
diff --git a/src/equil/MultiPhaseEquil.cpp b/src/equil/MultiPhaseEquil.cpp
index adf6e8eb3..f29cf9027 100644
--- a/src/equil/MultiPhaseEquil.cpp
+++ b/src/equil/MultiPhaseEquil.cpp
@@ -30,17 +30,15 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_
     m_incl_element.resize(m_nel_mix,1);
     for (m = 0; m < m_nel_mix; m++) {
         string enm = mix->elementName(m);
-        // element 'E' or 'e' represents an electron; this
-        // requires special handling, so save its index
-        // for later use
+        // element 'E' or 'e' represents an electron; this requires special
+        // handling, so save its index for later use
         if (enm == "E" || enm == "e") {
             m_eloc = m;
         }
-        // if an element other than electrons is not present in
-        // the mixture, then exclude it and all species containing
-        // it from the calculation. Electrons are a special case,
-        // since a species can have a negative number of 'atoms'
-        // of electrons (positive ions).
+        // if an element other than electrons is not present in the mixture,
+        // then exclude it and all species containing it from the calculation.
+        // Electrons are a special case, since a species can have a negative
+        // number of 'atoms' of electrons (positive ions).
         if (m_mix->elementMoles(m) <= 0.0 && m != m_eloc) {
             m_incl_element[m] = 0;
             for (k = 0; k < m_nsp_mix; k++) {
@@ -65,17 +63,15 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_
         }
     }
 
-    // include pure single-constituent phases only if their thermo
-    // data are valid for this temperature. This is necessary,
-    // since some thermo polynomial fits are done only for a
-    // limited temperature range. For example, using the NASA
-    // polynomial fits for solid ice and liquid water, if this
-    // were not done the calculation would predict solid ice to be
-    // present far above its melting point, since the thermo
-    // polynomial fits only extend to 273.15 K, and give
-    // unphysical results above this temperature, leading
-    // (incorrectly) to Gibbs free energies at high temperature
-    // lower than for liquid water.
+    // include pure single-constituent phases only if their thermo data are
+    // valid for this temperature. This is necessary, since some thermo
+    // polynomial fits are done only for a limited temperature range. For
+    // example, using the NASA polynomial fits for solid ice and liquid water,
+    // if this were not done the calculation would predict solid ice to be
+    // present far above its melting point, since the thermo polynomial fits
+    // only extend to 273.15 K, and give unphysical results above this
+    // temperature, leading (incorrectly) to Gibbs free energies at high
+    // temperature lower than for liquid water.
     size_t ip;
     for (k = 0; k < m_nsp_mix; k++) {
         ip = m_mix->speciesPhaseIndex(k);
@@ -92,8 +88,7 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_
         }
     }
 
-    // Now build the list of all species to be included in the
-    // calculation.
+    // Now build the list of all species to be included in the calculation.
     for (k = 0; k < m_nsp_mix; k++) {
         if (m_incl_species[k] ==1) {
             m_nsp++;
@@ -131,10 +126,9 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_
         m_order[k] = k;
     }
 
-    // if the 'start' flag is set, estimate the initial mole
-    // numbers by doing a linear Gibbs minimization. In this case,
-    // only the elemental composition of the initial mixture state
-    // matters.
+    // if the 'start' flag is set, estimate the initial mole numbers by doing a
+    // linear Gibbs minimization. In this case, only the elemental composition
+    // of the initial mixture state matters.
     if (start) {
         setInitialMoles(loglevel-1);
     }
@@ -160,10 +154,9 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_
     m_force = false;
     updateMixMoles();
 
-    // At this point, the instance has been created, the species
-    // to be included have been determined, and an initial
-    // composition has been selected that has all non-zero mole
-    // numbers for the included species.
+    // At this point, the instance has been created, the species to be included
+    // have been determined, and an initial composition has been selected that
+    // has all non-zero mole numbers for the included species.
 }
 
 doublereal MultiPhaseEquil::equilibrate(int XY, doublereal err,
@@ -219,8 +212,7 @@ int MultiPhaseEquil::setInitialMoles(int loglevel)
     int iter = 0;
 
     while (redo) {
-        // choose a set of components based on the current
-        // composition
+        // choose a set of components based on the current composition
         computeN();
         redo = false;
         iter++;
@@ -331,18 +323,16 @@ void MultiPhaseEquil::getComponents(const std::vector<size_t>& order)
             }
         }
 
-        // If a pivot is zero, exchange columns.  This occurs when
-        // a species has an elemental composition that is not
-        // linearly independent of the component species that have
-        // already been assigned
+        // If a pivot is zero, exchange columns.  This occurs when a species has
+        // an elemental composition that is not linearly independent of the
+        // component species that have already been assigned
         if (m < nColumns && m_A(m,m) == 0.0) {
-            // First, we need to find a good candidate for a
-            // component species to swap in for the one that has
-            // zero pivot. It must contain element m, be linearly
-            // independent of the components processed so far
-            // (m_A(m,k) != 0), and should be a major species if
-            // possible. We'll choose the species with greatest
-            // mole fraction that satisfies these criteria.
+            // First, we need to find a good candidate for a component species
+            // to swap in for the one that has zero pivot. It must contain
+            // element m, be linearly independent of the components processed so
+            // far (m_A(m,k) != 0), and should be a major species if possible.
+            // We'll choose the species with greatest mole fraction that
+            // satisfies these criteria.
             doublereal maxmoles = -999.0;
             size_t kmax = 0;
             for (k = m+1; k < nColumns; k++) {
@@ -480,10 +470,9 @@ doublereal MultiPhaseEquil::stepComposition(int loglevel)
         } else {
             FCTR = 0.9;
         }
-        // if species k is in a multi-species solution phase, then its
-        // mole number must remain positive, unless the entire phase
-        // goes away. First we'll determine an upper bound on omega,
-        // such that all
+        // if species k is in a multi-species solution phase, then its mole
+        // number must remain positive, unless the entire phase goes away. First
+        // we'll determine an upper bound on omega, such that all
         if (m_dsoln[k] == 1) {
             if ((m_moles[k] > MAJOR_THRESHOLD) || (ik < m_nel)) {
                 if (m_moles[k] < MAJOR_THRESHOLD) {
@@ -516,9 +505,9 @@ doublereal MultiPhaseEquil::stepComposition(int loglevel)
 
     // now take a step with this scaled omega
     step(omegamax, m_work);
-    // compute the gradient of G at this new position in the
-    // current direction. If it is positive, then we have overshot
-    // the minimum. In this case, interpolate back.
+    // compute the gradient of G at this new position in the current direction.
+    // If it is positive, then we have overshot the minimum. In this case,
+    // interpolate back.
     doublereal not_mu = 1.0e12;
     m_mix->getValidChemPotentials(not_mu, m_mu.data());
     doublereal grad1 = 0.0;