Fix SurfPhase::setCoveragesNoNorm to retain negative coverages

This prevents Jacobian singularities due to (trivial) negative coverages,
e.g. in ReactingSurf1D.

include/cantera/thermo/Phase.h

@@ -514,6 +514,9 @@ public:
     //!                     is the number of species in the phase.
     virtual void setConcentrations(const doublereal* const conc);
 
+    //! Set the concentrations without ignoring negative concentrations
+    virtual void setConcentrationsNoNorm(const double* const conc);
+
     //! Elemental mass fraction of element m
     /*!
      *  The elemental mass fraction \f$Z_{\mathrm{mass},m}\f$ of element \f$m\f$

src/thermo/Phase.cpp

@@ -614,6 +614,23 @@ void Phase::setConcentrations(const doublereal* const conc)
     m_stateNum++;
 }
 
+void Phase::setConcentrationsNoNorm(const double* const conc)
+{
+    doublereal sum = 0.0, norm = 0.0;
+    for (size_t k = 0; k != m_kk; ++k) {
+        sum += conc[k] * m_molwts[k];
+        norm += conc[k];
+    }
+    m_mmw = sum/norm;
+    setDensity(sum);
+    doublereal rsum = 1.0/sum;
+    for (size_t k = 0; k != m_kk; ++k) {
+        m_ym[k] = conc[k] * rsum;
+        m_y[k] = m_ym[k] * m_molwts[k];
+    }
+    m_stateNum++;
+}
+
 doublereal Phase::elementalMassFraction(const size_t m) const
 {
     checkElementIndex(m);

src/thermo/SurfPhase.cpp

@@ -281,8 +281,7 @@ void SurfPhase::setCoveragesNoNorm(const doublereal* theta)
     for (size_t k = 0; k < m_kk; k++) {
         m_work[k] = m_n0*theta[k]/size(k);
     }
-    // Call the Phase:: class function setConcentrations.
-    setConcentrations(m_work.data());
+    setConcentrationsNoNorm(m_work.data());
 }
 
 void SurfPhase::getCoverages(doublereal* theta) const