diff --git a/include/cantera/thermo/HMWSoln.h b/include/cantera/thermo/HMWSoln.h
index df300e686..8ca3281b8 100644
--- a/include/cantera/thermo/HMWSoln.h
+++ b/include/cantera/thermo/HMWSoln.h
@@ -3272,6 +3272,9 @@ private:
      */
     void calc_lambdas(double is) const;
 
+    //! Track the last ionic strength lambdas were calculated at to avoid unnecessarily recalculating them
+    mutable double last_is;
+
     /**
      *  Calculate etheta and etheta_prime
      *
diff --git a/include/cantera/thermo/MaskellSolidSolnPhase.h b/include/cantera/thermo/MaskellSolidSolnPhase.h
index c59a21ce9..f35219148 100644
--- a/include/cantera/thermo/MaskellSolidSolnPhase.h
+++ b/include/cantera/thermo/MaskellSolidSolnPhase.h
@@ -65,17 +65,17 @@ public:
     //! Return the standard concentration for the kth species
     /*!
      * The standard concentration \f$ C^0_k \f$ used to normalize the
-     * generalized concentration. In many cases, this quantity will be the
-     * same for all species in a phase. However, for this case, we will return
-     * a distinct concentration for each species. This is the inverse of the
-     * species molar volume. Units for the standard concentration are kmol
-     * m<SUP>-3</SUP>.
+     * generalized concentration.
      *
-     * @param k Species number: this is a require parameter,
-     * a change from the ThermoPhase base class, where it was
-     * an optional parameter.
+     * @param k Species number: this is an optional parameter,
      */
-    virtual doublereal standardConcentration(size_t k) const;
+    virtual doublereal standardConcentration(size_t k=0) const { return 1.0; }
+
+    //! Natural logarithm of the standard concentration of the kth species.
+    /*!
+     * @param k    index of the species (defaults to zero)
+     */
+    virtual doublereal logStandardConc(size_t k=0) const { return 0.0; }
 
     //! @name Molar Thermodynamic Properties of the Solution
     //! @{
@@ -294,7 +294,12 @@ protected:
      * Function to call through to m_spthermo->update and fill m_h0_RT,
      * m_cp0_R, m_g0_RT, m_s0_R.
      */
-    void _updateThermo() const;
+    bool _updateThermo() const;
+
+    //! Vector containing the last computed activity coefficients at T = m_tlast and r = last_r
+    mutable std::vector<doublereal> last_ac;
+    //! Last value of r used to update activity coeffs
+    mutable doublereal last_r;
 
     //! Vector containing the species reference enthalpies at T = m_tlast
     mutable vector_fp      m_h0_RT;
diff --git a/src/thermo/HMWSoln.cpp b/src/thermo/HMWSoln.cpp
index 10a826f4b..c3dff31e8 100644
--- a/src/thermo/HMWSoln.cpp
+++ b/src/thermo/HMWSoln.cpp
@@ -73,6 +73,7 @@ HMWSoln::HMWSoln() :
     CROP_ln_gamma_o_max(3.0),
     CROP_ln_gamma_k_min(-5.0),
     CROP_ln_gamma_k_max(15.0),
+    last_is(-1.0),
     m_debugCalc(0)
 {
     for (size_t i = 0; i < 17; i++) {
@@ -126,6 +127,7 @@ HMWSoln::HMWSoln(const std::string& inputFile, const std::string& id_) :
     CROP_ln_gamma_o_max(3.0),
     CROP_ln_gamma_k_min(-5.0),
     CROP_ln_gamma_k_max(15.0),
+    last_is(-1.0),
     m_debugCalc(0)
 {
     for (int i = 0; i < 17; i++) {
@@ -180,6 +182,7 @@ HMWSoln::HMWSoln(XML_Node& phaseRoot, const std::string& id_) :
     CROP_ln_gamma_o_max(3.0),
     CROP_ln_gamma_k_min(-5.0),
     CROP_ln_gamma_k_max(15.0),
+    last_is(-1.0),
     m_debugCalc(0)
 {
     for (int i = 0; i < 17; i++) {
@@ -234,6 +237,7 @@ HMWSoln::HMWSoln(const HMWSoln& b) :
     CROP_ln_gamma_o_max(3.0),
     CROP_ln_gamma_k_min(-5.0),
     CROP_ln_gamma_k_max(15.0),
+    last_is(-1.0),
     m_debugCalc(0)
 {
     /*
@@ -447,6 +451,7 @@ HMWSoln::HMWSoln(int testProb) :
     CROP_ln_gamma_o_max(3.0),
     CROP_ln_gamma_k_min(-5.0),
     CROP_ln_gamma_k_max(15.0),
+    last_is(-1.0),
     m_debugCalc(0)
 {
     if (testProb != 1) {
@@ -1845,7 +1850,10 @@ void HMWSoln::s_updatePitzer_CoeffWRTemp(int doDerivs) const
         }
     }
 
-
+    const double twoT = 2.0 * T;
+    const double invT = 1.0 / T;
+    const double invT2 = invT * invT;
+    const double twoinvT3 = 2.0 * invT * invT2;
     for (i = 1; i < m_kk; i++) {
         for (j = 1; j < m_kk; j++) {
             for (size_t k = 1; k < m_kk; k++) {
@@ -1868,14 +1876,14 @@ void HMWSoln::s_updatePitzer_CoeffWRTemp(int doDerivs) const
                                    + Psi_coeff[4]*tln;
 
                     m_Psi_ijk_L[n] = Psi_coeff[1]
-                                     + Psi_coeff[2]*2.0*T
-                                     - Psi_coeff[3]/(T*T)
-                                     + Psi_coeff[4]/T;
+                                     + Psi_coeff[2]*twoT
+                                     - Psi_coeff[3]*invT2
+                                     + Psi_coeff[4]*invT;
 
                     m_Psi_ijk_LL[n] =
                         Psi_coeff[2]*2.0
-                        + 2.0*Psi_coeff[3]/(T*T*T)
-                        - Psi_coeff[4]/(T*T);
+                        + Psi_coeff[3]*twoinvT3
+                        - Psi_coeff[4]*invT2;
                 }
             }
         }
@@ -5456,6 +5464,9 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
 
 void HMWSoln::calc_lambdas(double is) const
 {
+    const double tol = 1.e-12;
+    if( std::abs(is - last_is) < tol ) return;
+
     double aphi, dj, jfunc, jprime, t, x, zprod;
     int i, ij, j;
     /*
@@ -5513,6 +5524,7 @@ void HMWSoln::calc_lambdas(double is) const
 #endif
         }
     }
+    last_is = is;
 }
 
 void HMWSoln::calc_thetas(int z1, int z2,
diff --git a/src/thermo/MaskellSolidSolnPhase.cpp b/src/thermo/MaskellSolidSolnPhase.cpp
index 3aee9aba0..67b71db2c 100644
--- a/src/thermo/MaskellSolidSolnPhase.cpp
+++ b/src/thermo/MaskellSolidSolnPhase.cpp
@@ -23,6 +23,8 @@ namespace Cantera
 MaskellSolidSolnPhase::MaskellSolidSolnPhase() :
     m_Pref(OneAtm),
     m_Pcurrent(OneAtm),
+    last_ac(2),
+    last_r(-1.0),
     m_h0_RT(2),
     m_cp0_R(2),
     m_g0_RT(2),
@@ -36,6 +38,8 @@ MaskellSolidSolnPhase::MaskellSolidSolnPhase() :
 MaskellSolidSolnPhase::MaskellSolidSolnPhase(const MaskellSolidSolnPhase& b) :
     m_Pref(OneAtm),
     m_Pcurrent(OneAtm),
+    last_ac(2),
+    last_r(-1.0),
     m_h0_RT(2),
     m_cp0_R(2),
     m_g0_RT(2),
@@ -64,23 +68,13 @@ ThermoPhase* MaskellSolidSolnPhase::duplMyselfAsThermoPhase() const
 void MaskellSolidSolnPhase::
 getActivityConcentrations(doublereal* c) const
 {
-    std::vector<doublereal> pmv(m_kk);
-    getPartialMolarVolumes(&pmv[0]);
     getActivityCoefficients(c);
     for(unsigned sp=0; sp < m_kk; ++sp)
     {
-      c[sp] *= moleFraction(sp) / pmv[sp];
+      c[sp] *= moleFraction(sp);
     }
 }
 
-doublereal MaskellSolidSolnPhase::standardConcentration(size_t k) const
-{
-    std::vector<doublereal> pmv(m_kk);
-    getPartialMolarVolumes(&pmv[0]);
-    doublereal result = 1.0 / pmv[k];
-    return result;
-}
-
 /********************************************************************
  *            Molar Thermodynamic Properties of the Solution
  ********************************************************************/
@@ -163,17 +157,23 @@ void MaskellSolidSolnPhase::setMolarDensity(const doublereal n)
 void MaskellSolidSolnPhase::
 getActivityCoefficients(doublereal* ac) const
 {
-    _updateThermo();
+    bool temp_changed = _updateThermo();
     const doublereal r = moleFraction(product_species_index);
-    const doublereal pval = p(r);
-    const doublereal fmval = fm(r);
-    const doublereal rfm = r * fmval;
-    const doublereal RT = GasConstant * temperature();
-    const doublereal A = (std::pow(1 - rfm, pval) * std::pow(rfm, pval) * std::pow(r - rfm, 1 - pval))  /
-                         (std::pow(1 - r - rfm, 1 + pval) * (1 - r));
-    const doublereal B = pval * h_mixing / RT;
-    ac[product_species_index] = A * std::exp(B);
-    ac[reactant_species_index] = 1 / (A * r * (1-r) ) * std::exp(-B);
+    const doublereal tol = 1.e-12;
+    if( temp_changed || std::abs(r - last_r) > tol )
+    {
+      const doublereal pval = p(r);
+      const doublereal fmval = fm(r);
+      const doublereal rfm = r * fmval;
+      const doublereal RT = GasConstant * temperature();
+      const doublereal A = (std::pow(1 - rfm, pval) * std::pow(rfm, pval) * std::pow(r - rfm, 1 - pval))  /
+                           (std::pow(1 - r - rfm, 1 + pval) * (1 - r));
+      const doublereal B = pval * h_mixing / RT;
+      last_ac[product_species_index] = A * std::exp(B);
+      last_ac[reactant_species_index] = 1 / (A * r * (1-r) ) * std::exp(-B);
+      last_r = r;
+    }
+    std::copy(last_ac.begin(), last_ac.end(), ac);
 }
 
 void MaskellSolidSolnPhase::
@@ -324,7 +324,7 @@ void MaskellSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string
     VPStandardStateTP::initThermoXML(phaseNode, id_);
 }
 
-void MaskellSolidSolnPhase::_updateThermo() const
+bool MaskellSolidSolnPhase::_updateThermo() const
 {
     assert(m_kk == 2);
     doublereal tnow = temperature();
@@ -339,7 +339,9 @@ void MaskellSolidSolnPhase::_updateThermo() const
             m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k];
         }
         m_tlast = tnow;
+        return true;
     }
+    return false;
 }
 
 doublereal MaskellSolidSolnPhase::s() const
diff --git a/test/thermo/MaskellSolidSolnPhase_Test.cpp b/test/thermo/MaskellSolidSolnPhase_Test.cpp
index be2cd2684..59bd3dbcc 100644
--- a/test/thermo/MaskellSolidSolnPhase_Test.cpp
+++ b/test/thermo/MaskellSolidSolnPhase_Test.cpp
@@ -123,8 +123,8 @@ TEST_F(MaskellSolidSolnPhase_Test, standardConcentrations)
     initializeTestPhaseWithXML(valid_file);
     ASSERT_TRUE(dynamic_cast<MaskellSolidSolnPhase *>(test_phase) != NULL);
 
-    EXPECT_DOUBLE_EQ(1.0 / 0.005, test_phase->standardConcentration(0));
-    EXPECT_DOUBLE_EQ(1.0 / 0.01, test_phase->standardConcentration(1));
+    EXPECT_DOUBLE_EQ(1.0, test_phase->standardConcentration(0));
+    EXPECT_DOUBLE_EQ(1.0, test_phase->standardConcentration(1));
 }
 
 TEST_F(MaskellSolidSolnPhase_Test, activityConcentrations)