diff --git a/include/cantera/thermo/HMWSoln.h b/include/cantera/thermo/HMWSoln.h
index c6303b07f..e376c1f81 100644
--- a/include/cantera/thermo/HMWSoln.h
+++ b/include/cantera/thermo/HMWSoln.h
@@ -2362,6 +2362,9 @@ public:
      */
     int m_form_A_Debye;
 
+protected:
+    virtual void invalidateCachedDataOnStateChange(StateVariable changed_var);
+
 private:
     /**
      * A_Debye -> this expression appears on the top of the
@@ -2394,14 +2397,10 @@ private:
      *                    dw = C_0 * M_0 (density of water) (kg/m3)
      *                       = 1.0E3 at 25C
      */
+    mutable bool m_A_Debye_valid;
     mutable double m_A_Debye;
-    mutable double m_last_dA_DebyedP_TP;
-    mutable double m_last_dA_DebyedP_TP_T;
-    mutable double m_last_dA_DebyedP_TP_P;
-
-    mutable double m_last_P;
-    mutable double m_last_T;
-    mutable double m_last_A_Debye;
+    mutable bool m_dA_DebyedP_TP_valid;
+    mutable double m_dA_DebyedP_TP;
 
     //!  Water standard state calculator
     /*!
@@ -3173,8 +3172,7 @@ private:
      * natural logarithm of the molality activity coefficients
      */
     void s_update_lnMolalityActCoeff() const;
-    mutable std::vector<doublereal> m_last_cropped_molalities;
-    bool cropped_molalities_changed() const;
+    mutable bool s_update_lnMolalityActCoeff_valid;
 
     //! This function calculates the temperature derivative of the
     //! natural logarithm of the molality activity coefficients.
@@ -3183,6 +3181,7 @@ private:
      * coefficient is on the molality scale. It's derivative is too.
      */
     void s_update_dlnMolalityActCoeff_dT() const;
+    mutable bool s_update_dlnMolalityActCoeff_dT_valid;
 
     /**
      * This function calculates the temperature second derivative
@@ -3190,6 +3189,7 @@ private:
      * coefficients.
      */
     void s_update_d2lnMolalityActCoeff_dT2() const;
+    mutable bool s_update_d2lnMolalityActCoeff_dT2_valid;
 
     /**
      * This function calculates the pressure derivative of the
@@ -3198,6 +3198,7 @@ private:
      * Assumes that the activity coefficients are current.
      */
     void s_update_dlnMolalityActCoeff_dP() const;
+    mutable bool s_update_dlnMolalityActCoeff_dP_valid;
 
     //! This function will be called to update the internally stored
     //! natural logarithm of the molality activity coefficients
@@ -3270,9 +3271,7 @@ private:
      * @param is Ionic strength
      */
     void calc_lambdas(double is) const;
-
-    //! Track the last ionic strength lambdas were calculated at to avoid unnecessarily recalculating them
-    mutable double m_last_is;
+    mutable doublereal m_last_is;
 
     /**
      *  Calculate etheta and etheta_prime
diff --git a/include/cantera/thermo/MaskellSolidSolnPhase.h b/include/cantera/thermo/MaskellSolidSolnPhase.h
index 5c03ac4cb..9fe8600ec 100644
--- a/include/cantera/thermo/MaskellSolidSolnPhase.h
+++ b/include/cantera/thermo/MaskellSolidSolnPhase.h
@@ -274,7 +274,11 @@ public:
 
 
     void set_h_mix(const doublereal hmix) { h_mixing = hmix; }
+
 protected:
+    void invalidateCachedDataOnStateChange(StateVariable changed_var);
+
+private:
     /**
      * Value of the reference pressure for all species in this phase.
      * The T dependent polynomials are evaluated at the reference
@@ -294,12 +298,12 @@ protected:
      * Function to call through to m_spthermo->update and fill m_h0_RT,
      * m_cp0_R, m_g0_RT, m_s0_R.
      */
-    bool _updateThermo() const;
+    void _updateThermo() const;
+    mutable bool m_updateThermo_valid;
 
     //! Vector containing the last computed activity coefficients at T = m_tlast and r = last_r
-    mutable std::vector<doublereal> m_last_ac;
-    //! Last value of r used to update activity coeffs
-    mutable doublereal m_last_mole_frac_product;
+    mutable std::vector<doublereal> m_activity_coeffs;
+    mutable bool m_activity_coeffs_valid;
 
     //! Vector containing the species reference enthalpies at T = m_tlast
     mutable vector_fp      m_h0_RT;
@@ -323,7 +327,6 @@ protected:
     int product_species_index;
     int reactant_species_index;
 
-private:
     // Functions to calculate some of the pieces of the mixing terms.
     doublereal s() const;
     doublereal fm(const doublereal r) const;
diff --git a/include/cantera/thermo/Phase.h b/include/cantera/thermo/Phase.h
index 4e86780a3..13918b899 100644
--- a/include/cantera/thermo/Phase.h
+++ b/include/cantera/thermo/Phase.h
@@ -551,7 +551,14 @@ public:
             throw CanteraError("Phase::setDensity()", "density must be positive");
         }
         m_dens = density_;
+        if(density_ != m_last_dens) {
+            invalidateCachedDataOnStateChange(DENSITY);
+            m_last_dens = density_;
+        }
     }
+private:
+    mutable doublereal m_last_dens;
+public:
 
     //! Set the internally stored molar density (kmol/m^3) of the phase.
     //!     @param[in] molarDensity Input molar density (kmol/m^3).
@@ -565,8 +572,15 @@ public:
                                "temperature must be positive");
         }
         m_temp = temp;
+        if(temp != m_last_temp) {
+            invalidateCachedDataOnStateChange(TEMPERATURE);
+            m_last_temp = temp;
+        }
     }
     //@}
+private:
+    mutable doublereal m_last_temp;
+public:
 
     //! @name Mean Properties
     //!@{
@@ -690,6 +704,35 @@ public:
     }
 
 protected:
+    //! Virtual function to simplify caching of extensive computations in child classes.
+    /*!
+     *  The intent is that this virtual function should be called any time a state variable
+     *  of a ThermoPhase class is changed the class can mark any cached data that is expensive to
+     *  compute as invalid. For example:
+     *  class APhase : public Phase
+     *  {
+     *  public:
+     *    void update_activity_coeffs(); // Expensive function that only needs to run on state change
+     *    void update_activity_coeffs_dT(); // Another expensive function
+     *  private:
+     *    // Variables that track whether the state has been changed since the expensive functions
+     *    // were last run to avoid unnecessary work. These will be updated to false by invalidateCachedDataOnStateChange()
+     *    bool m_activity_coeffs_valid;
+     *    bool m_activity_coeffs_dT_valid;
+     *  }
+     *
+     *  This is somewhat easier to manage than having each individual expensive function track the last T, P, X_i, etc
+     *  that it was calculated at to determine whether it is safe to skip the expensive computation.
+     */
+    enum StateVariable {
+      PRESSURE,
+      TEMPERATURE,
+      SPECIES,
+      DENSITY,
+      OTHER
+    };
+    virtual void invalidateCachedDataOnStateChange(StateVariable changed_var) {}
+
     //! Set the molecular weight of a single species to a given value
     //!     @param k       id of the species
     //!     @param mw      Molecular Weight (kg kmol-1)
@@ -741,6 +784,8 @@ private:
     mutable vector_fp m_ym;
 
     mutable vector_fp m_y; //!< species mass fractions
+    mutable vector_fp m_last_y;
+    bool massFractionsChanged();
 
     vector_fp m_molwts; //!< species molecular weights (kg kmol-1)
 
diff --git a/include/cantera/thermo/ThermoPhase.h b/include/cantera/thermo/ThermoPhase.h
index baa56d137..efd248c71 100644
--- a/include/cantera/thermo/ThermoPhase.h
+++ b/include/cantera/thermo/ThermoPhase.h
@@ -1131,6 +1131,7 @@ private:
     //! Helper function used by setState_HPorUV and setState_SPorSV.
     //! Sets the temperature and (if set_p is true) the pressure.
     void setState_conditional_TP(doublereal t, doublereal p, bool set_p);
+
 public:
 
     //@}
diff --git a/src/thermo/HMWSoln.cpp b/src/thermo/HMWSoln.cpp
index ae28ac989..a0c631cdf 100644
--- a/src/thermo/HMWSoln.cpp
+++ b/src/thermo/HMWSoln.cpp
@@ -38,12 +38,10 @@ HMWSoln::HMWSoln() :
     m_TempPitzerRef(298.15),
     m_IionicMolalityStoich(0.0),
     m_form_A_Debye(A_DEBYE_WATER),
+    m_A_Debye_valid(false),
     m_A_Debye(1.172576),   // units = sqrt(kg/gmol)
-    m_last_dA_DebyedP_TP(-1.0),
-    m_last_dA_DebyedP_TP_T(-1.0),
-    m_last_dA_DebyedP_TP_P(-1.0),
-    m_last_P(-1.0),
-    m_last_T(-1.0),
+    m_dA_DebyedP_TP_valid(false),
+    m_dA_DebyedP_TP(-1.0),
     m_waterSS(0),
     m_densWaterSS(1000.),
     m_waterProps(0),
@@ -75,6 +73,10 @@ HMWSoln::HMWSoln() :
     CROP_ln_gamma_o_max(3.0),
     CROP_ln_gamma_k_min(-5.0),
     CROP_ln_gamma_k_max(15.0),
+    s_update_lnMolalityActCoeff_valid(false),
+    s_update_dlnMolalityActCoeff_dT_valid(false),
+    s_update_d2lnMolalityActCoeff_dT2_valid(false),
+    s_update_dlnMolalityActCoeff_dP_valid(false),
     m_last_is(-1.0),
     m_debugCalc(0)
 {
@@ -94,12 +96,10 @@ HMWSoln::HMWSoln(const std::string& inputFile, const std::string& id_) :
     m_TempPitzerRef(298.15),
     m_IionicMolalityStoich(0.0),
     m_form_A_Debye(A_DEBYE_WATER),
+    m_A_Debye_valid(false),
     m_A_Debye(1.172576),   // units = sqrt(kg/gmol)
-    m_last_dA_DebyedP_TP(-1.0),
-    m_last_dA_DebyedP_TP_T(-1.0),
-    m_last_dA_DebyedP_TP_P(-1.0),
-    m_last_P(-1.0),
-    m_last_T(-1.0),
+    m_dA_DebyedP_TP_valid(false),
+    m_dA_DebyedP_TP(-1.0),
     m_waterSS(0),
     m_densWaterSS(1000.),
     m_waterProps(0),
@@ -131,6 +131,10 @@ 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),
+    s_update_lnMolalityActCoeff_valid(false),
+    s_update_dlnMolalityActCoeff_dT_valid(false),
+    s_update_d2lnMolalityActCoeff_dT2_valid(false),
+    s_update_dlnMolalityActCoeff_dP_valid(false),
     m_last_is(-1.0),
     m_debugCalc(0)
 {
@@ -151,12 +155,10 @@ HMWSoln::HMWSoln(XML_Node& phaseRoot, const std::string& id_) :
     m_TempPitzerRef(298.15),
     m_IionicMolalityStoich(0.0),
     m_form_A_Debye(A_DEBYE_WATER),
+    m_A_Debye_valid(false),
     m_A_Debye(1.172576),   // units = sqrt(kg/gmol)
-    m_last_dA_DebyedP_TP(-1.0),
-    m_last_dA_DebyedP_TP_T(-1.0),
-    m_last_dA_DebyedP_TP_P(-1.0),
-    m_last_P(-1.0),
-    m_last_T(-1.0),
+    m_dA_DebyedP_TP_valid(false),
+    m_dA_DebyedP_TP(-1.0),
     m_waterSS(0),
     m_densWaterSS(1000.),
     m_waterProps(0),
@@ -188,6 +190,10 @@ 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),
+    s_update_lnMolalityActCoeff_valid(false),
+    s_update_dlnMolalityActCoeff_dT_valid(false),
+    s_update_d2lnMolalityActCoeff_dT2_valid(false),
+    s_update_dlnMolalityActCoeff_dP_valid(false),
     m_last_is(-1.0),
     m_debugCalc(0)
 {
@@ -208,12 +214,10 @@ HMWSoln::HMWSoln(const HMWSoln& b) :
     m_TempPitzerRef(298.15),
     m_IionicMolalityStoich(0.0),
     m_form_A_Debye(A_DEBYE_WATER),
+    m_A_Debye_valid(false),
     m_A_Debye(1.172576),   // units = sqrt(kg/gmol)
-    m_last_dA_DebyedP_TP(-1.0),
-    m_last_dA_DebyedP_TP_T(-1.0),
-    m_last_dA_DebyedP_TP_P(-1.0),
-    m_last_P(-1.0),
-    m_last_T(-1.0),
+    m_dA_DebyedP_TP_valid(false),
+    m_dA_DebyedP_TP(-1.0),
     m_waterSS(0),
     m_densWaterSS(1000.),
     m_waterProps(0),
@@ -245,6 +249,10 @@ 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),
+    s_update_lnMolalityActCoeff_valid(false),
+    s_update_dlnMolalityActCoeff_dT_valid(false),
+    s_update_d2lnMolalityActCoeff_dT2_valid(false),
+    s_update_dlnMolalityActCoeff_dP_valid(false),
     m_last_is(-1.0),
     m_debugCalc(0)
 {
@@ -424,12 +432,10 @@ HMWSoln::HMWSoln(int testProb) :
     m_TempPitzerRef(298.15),
     m_IionicMolalityStoich(0.0),
     m_form_A_Debye(A_DEBYE_WATER),
+    m_A_Debye_valid(false),
     m_A_Debye(1.172576),   // units = sqrt(kg/gmol)
-    m_last_dA_DebyedP_TP(-1.0),
-    m_last_dA_DebyedP_TP_T(-1.0),
-    m_last_dA_DebyedP_TP_P(-1.0),
-    m_last_P(-1.0),
-    m_last_T(-1.0),
+    m_dA_DebyedP_TP_valid(false),
+    m_dA_DebyedP_TP(-1.0),
     m_waterSS(0),
     m_densWaterSS(1000.),
     m_waterProps(0),
@@ -461,6 +467,10 @@ 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),
+    s_update_lnMolalityActCoeff_valid(false),
+    s_update_dlnMolalityActCoeff_dT_valid(false),
+    s_update_d2lnMolalityActCoeff_dT2_valid(false),
+    s_update_dlnMolalityActCoeff_dP_valid(false),
     m_last_is(-1.0),
     m_debugCalc(0)
 {
@@ -1051,6 +1061,11 @@ doublereal HMWSoln::satPressure(doublereal t) {
 
 double HMWSoln::A_Debye_TP(double tempArg, double presArg) const
 {
+    if(m_A_Debye_valid) {
+      return m_A_Debye;
+    }
+    m_A_Debye_valid = true;
+
     double T = temperature();
     double A;
     if (tempArg != -1.0) {
@@ -1113,22 +1128,17 @@ double HMWSoln::dA_DebyedP_TP(double tempArg, double presArg) const
         P = presArg;
     }
     double dAdP;
-    const double tol = 1.e-6;
     switch (m_form_A_Debye) {
     case A_DEBYE_CONST:
         dAdP = 0.0;
         break;
     case A_DEBYE_WATER:
-        if( std::abs(T - m_last_dA_DebyedP_TP_T) > tol || std::abs(P - m_last_dA_DebyedP_TP_P) > tol )
-        {
-            dAdP = m_waterProps->ADebye(T, P, 3);
-            m_last_dA_DebyedP_TP_T = T;
-            m_last_dA_DebyedP_TP_P = P;
-            m_last_dA_DebyedP_TP  = dAdP;
-        }
-        else
-        {
-            dAdP = m_last_dA_DebyedP_TP;
+        if(!m_dA_DebyedP_TP_valid) {
+          dAdP = m_waterProps->ADebye(T, P, 3);
+          m_dA_DebyedP_TP  = dAdP;
+          m_dA_DebyedP_TP_valid = true;
+        } else {
+          dAdP = m_dA_DebyedP_TP;
         }
         break;
     default:
@@ -1347,30 +1357,13 @@ void HMWSoln::initLengths()
     counterIJ_setup();
 }
 
-bool HMWSoln::cropped_molalities_changed() const
-{
-    if( m_last_cropped_molalities.size() != m_kk ) {
-      m_last_cropped_molalities.resize(m_kk);
-      std::fill(m_last_cropped_molalities.begin(), m_last_cropped_molalities.end(), -1.0);
-    }
-
-    bool result = false;
-    const doublereal tol = 1.e-12;
-    for(size_t k=0; k < m_kk; ++k) {
-        if( std::abs(m_molalitiesCropped[k] - m_last_cropped_molalities[k]) > tol ) {
-          result = true;
-        }
-    }
-
-    if(result) {
-      std::copy( m_molalitiesCropped.begin(), m_molalitiesCropped.end(), m_last_cropped_molalities.begin() );
-    }
-
-    return result;
-}
-
 void HMWSoln::s_update_lnMolalityActCoeff() const
 {
+    if(s_update_lnMolalityActCoeff_valid) {
+      return;
+    }
+    s_update_lnMolalityActCoeff_valid = true;
+
     /*
      * Calculate the molalities. Currently, the molalities
      * may not be current with respect to the contents of the
@@ -1383,19 +1376,6 @@ void HMWSoln::s_update_lnMolalityActCoeff() const
      */
     calcMolalitiesCropped();
 
-    double T = temperature();
-    double P = pressure();
-    const doublereal tol = 1.e-12;
-    if( std::abs(m_A_Debye - m_last_A_Debye) < tol &&
-        std::abs(T - m_last_T) < tol &&
-        std::abs(P - m_last_P) < tol &&
-        !cropped_molalities_changed() ) {
-      return;
-    }
-    m_last_T = T;
-    m_last_P = P;
-    m_last_A_Debye = m_A_Debye;
-
     /*
      * Calculate the stoichiometric ionic charge. This isn't used in the
      * Pitzer formulation.
@@ -2964,6 +2944,11 @@ s_updatePitzer_lnMolalityActCoeff() const
 
 void HMWSoln::s_update_dlnMolalityActCoeff_dT() const
 {
+    if(s_update_dlnMolalityActCoeff_dT_valid) {
+      return;
+    }
+    s_update_dlnMolalityActCoeff_dT_valid = true;
+
     /*
      *  Zero the unscaled 2nd derivatives
      */
@@ -3813,6 +3798,11 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
 
 void HMWSoln::s_update_d2lnMolalityActCoeff_dT2() const
 {
+    if(s_update_d2lnMolalityActCoeff_dT2_valid) {
+      return;
+    }
+    s_update_d2lnMolalityActCoeff_dT2_valid = true;
+
     /*
      *  Zero the unscaled 2nd derivatives
      */
@@ -4667,6 +4657,10 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
 
 void HMWSoln::s_update_dlnMolalityActCoeff_dP() const
 {
+    if(s_update_dlnMolalityActCoeff_dP_valid) {
+      return;
+    }
+    s_update_dlnMolalityActCoeff_dP_valid = true;
     m_dlnActCoeffMolaldP_Unscaled.assign(m_kk, 0.0);
     s_updatePitzer_dlnMolalityActCoeff_dP();
 
@@ -5504,10 +5498,10 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
 
 void HMWSoln::calc_lambdas(double is) const
 {
-    const double tol = 1.e-12;
-    if( std::abs(is - m_last_is) < tol ) {
+    if( m_last_is == is ) {
       return;
     }
+    m_last_is = is;
 
     double aphi, dj, jfunc, jprime, t, x, zprod;
     int i, ij, j;
@@ -5566,7 +5560,6 @@ void HMWSoln::calc_lambdas(double is) const
 #endif
         }
     }
-    m_last_is = is;
 }
 
 void HMWSoln::calc_thetas(int z1, int z2,
@@ -5885,6 +5878,19 @@ doublereal HMWSoln::s_NBS_CLM_dlnMolalityActCoeff_dP() const
     return - dAdP * sqrtIs /(1.0 + 1.5 * sqrtIs);
 }
 
+void HMWSoln::invalidateCachedDataOnStateChange(StateVariable changed_var)
+{
+    if( changed_var == PRESSURE || changed_var == TEMPERATURE ) {
+      m_A_Debye_valid = false;
+      m_dA_DebyedP_TP_valid = false;
+    }
+    s_update_lnMolalityActCoeff_valid = false;
+    s_update_dlnMolalityActCoeff_dT_valid = false;
+    s_update_d2lnMolalityActCoeff_dT2_valid = false;
+    s_update_dlnMolalityActCoeff_dP_valid = false;
+    MolalityVPSSTP::invalidateCachedDataOnStateChange(changed_var);
+}
+
 int HMWSoln::debugPrinting()
 {
 #ifdef DEBUG_MODE
diff --git a/src/thermo/MaskellSolidSolnPhase.cpp b/src/thermo/MaskellSolidSolnPhase.cpp
index 96ef03ef5..32d5fd8f1 100644
--- a/src/thermo/MaskellSolidSolnPhase.cpp
+++ b/src/thermo/MaskellSolidSolnPhase.cpp
@@ -23,8 +23,9 @@ namespace Cantera
 MaskellSolidSolnPhase::MaskellSolidSolnPhase() :
     m_Pref(OneAtm),
     m_Pcurrent(OneAtm),
-    m_last_ac(2),
-    m_last_mole_frac_product(-1.0),
+    m_updateThermo_valid(false),
+    m_activity_coeffs(2),
+    m_activity_coeffs_valid(false),
     m_h0_RT(2),
     m_cp0_R(2),
     m_g0_RT(2),
@@ -38,8 +39,9 @@ MaskellSolidSolnPhase::MaskellSolidSolnPhase() :
 MaskellSolidSolnPhase::MaskellSolidSolnPhase(const MaskellSolidSolnPhase& b) :
     m_Pref(OneAtm),
     m_Pcurrent(OneAtm),
-    m_last_ac(2),
-    m_last_mole_frac_product(-1.0),
+    m_updateThermo_valid(false),
+    m_activity_coeffs(2),
+    m_activity_coeffs_valid(false),
     m_h0_RT(2),
     m_cp0_R(2),
     m_g0_RT(2),
@@ -157,10 +159,10 @@ void MaskellSolidSolnPhase::setMolarDensity(const doublereal n)
 void MaskellSolidSolnPhase::
 getActivityCoefficients(doublereal* ac) const
 {
-    bool temp_changed = _updateThermo();
-    const doublereal r = moleFraction(product_species_index);
-    const doublereal tol = 1.e-12;
-    if( temp_changed || std::abs(r - m_last_mole_frac_product) > tol ) {
+    _updateThermo();
+    if( !m_activity_coeffs_valid ) {
+      m_activity_coeffs_valid = true;
+      const doublereal r = moleFraction(product_species_index);
       const doublereal pval = p(r);
       const doublereal fmval = fm(r);
       const doublereal rfm = r * fmval;
@@ -168,11 +170,10 @@ getActivityCoefficients(doublereal* ac) const
       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;
-      m_last_ac[product_species_index] = A * std::exp(B);
-      m_last_ac[reactant_species_index] = 1 / (A * r * (1-r) ) * std::exp(-B);
-      m_last_mole_frac_product = r;
+      m_activity_coeffs[product_species_index] = A * std::exp(B);
+      m_activity_coeffs[reactant_species_index] = 1 / (A * r * (1-r) ) * std::exp(-B);
     }
-    std::copy(m_last_ac.begin(), m_last_ac.end(), ac);
+    std::copy(m_activity_coeffs.begin(), m_activity_coeffs.end(), ac);
 }
 
 void MaskellSolidSolnPhase::
@@ -323,14 +324,15 @@ void MaskellSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string
     VPStandardStateTP::initThermoXML(phaseNode, id_);
 }
 
-bool MaskellSolidSolnPhase::_updateThermo() const
+void MaskellSolidSolnPhase::_updateThermo() const
 {
     assert(m_kk == 2);
-    doublereal tnow = temperature();
-    if (m_tlast != tnow) {
+    if (!m_updateThermo_valid) {
+        m_updateThermo_valid = true;
         /*
          * Update the thermodynamic functions of the reference state.
          */
+        doublereal tnow = temperature();
         m_spthermo->update(tnow, DATA_PTR(m_cp0_R), DATA_PTR(m_h0_RT),
                            DATA_PTR(m_s0_R));
         m_tlast = tnow;
@@ -338,9 +340,7 @@ bool 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
@@ -360,4 +360,11 @@ doublereal MaskellSolidSolnPhase::p(const doublereal r) const
   return (1 - 2*r) / std::sqrt(sval*sval - 4 * sval * r + 4 * sval * r * r);
 }
 
+void MaskellSolidSolnPhase::invalidateCachedDataOnStateChange(StateVariable changed_var)
+{
+  m_updateThermo_valid = false;
+  m_activity_coeffs_valid = false;
+  VPStandardStateTP::invalidateCachedDataOnStateChange(changed_var);
+}
+
 }  // end namespace Cantera
diff --git a/src/thermo/Phase.cpp b/src/thermo/Phase.cpp
index ae77b5dbc..627f7443c 100644
--- a/src/thermo/Phase.cpp
+++ b/src/thermo/Phase.cpp
@@ -332,6 +332,9 @@ void Phase::setMoleFractions(const doublereal* const x)
      */
     m_mmw = sum/norm;
     m_stateNum++;
+    if(massFractionsChanged()) {
+        invalidateCachedDataOnStateChange(SPECIES);
+    }
 }
 
 void Phase::setMoleFractions_NoNorm(const doublereal* const x)
@@ -342,6 +345,9 @@ void Phase::setMoleFractions_NoNorm(const doublereal* const x)
     transform(m_ym.begin(), m_ym.begin() + m_kk, m_molwts.begin(),
               m_y.begin(), multiplies<double>());
     m_stateNum++;
+    if(massFractionsChanged()) {
+        invalidateCachedDataOnStateChange(SPECIES);
+    }
 }
 
 void Phase::setMoleFractionsByName(compositionMap& xMap)
@@ -376,6 +382,9 @@ void Phase::setMassFractions(const doublereal* const y)
               m_ym.begin(), multiplies<double>());
     m_mmw = 1.0 / accumulate(m_ym.begin(), m_ym.end(), 0.0);
     m_stateNum++;
+    if(massFractionsChanged()) {
+        invalidateCachedDataOnStateChange(SPECIES);
+    }
 }
 
 void Phase::setMassFractions_NoNorm(const doublereal* const y)
@@ -387,6 +396,9 @@ void Phase::setMassFractions_NoNorm(const doublereal* const y)
     sum = accumulate(m_ym.begin(), m_ym.end(), 0.0);
     m_mmw = 1.0/sum;
     m_stateNum++;
+    if( massFractionsChanged() ) {
+        invalidateCachedDataOnStateChange(SPECIES);
+    }
 }
 
 void Phase::setMassFractionsByName(compositionMap& yMap)
@@ -409,6 +421,22 @@ void Phase::setMassFractionsByName(const std::string& y)
     setMassFractionsByName(c);
 }
 
+bool Phase::massFractionsChanged()
+{
+  if(m_last_y.size() != m_kk) {
+    m_last_y.resize(m_kk);
+    std::fill(m_last_y.begin(), m_last_y.end(), -1.0);
+  }
+
+  bool result = false;
+  for(size_t k=0; k < m_kk; ++k) {
+    if(m_y[k] != m_last_y[k]) {
+      result = true;
+    }
+  }
+  return result;
+}
+
 void Phase::setState_TRX(doublereal t, doublereal dens, const doublereal* x)
 {
     setMoleFractions(x);