diff --git a/include/AbstractState.h b/include/AbstractState.h
index 3e7ed1bc..c4220af0 100644
--- a/include/AbstractState.h
+++ b/include/AbstractState.h
@@ -234,7 +234,7 @@ protected:
     
     /// Using this backend, calculate a phase given by the state string
     /// @param state A string that describes the state desired, one of "hs_anchor", "critical"/"crit", "reducing"
-    virtual const CoolProp::SimpleState calc_state(const std::string &state){throw NotImplementedError("calc_state is not implemented for this backend");};
+    virtual const CoolProp::SimpleState & calc_state(const std::string &state){throw NotImplementedError("calc_state is not implemented for this backend");};
 
 public:
 
@@ -397,7 +397,7 @@ public:
     /// @param Q The quality for the parameter that is given (0 = saturated liquid, 1 = saturated vapor)
     /// @param given The key for the parameter that is given
     /// @param value The value for the parameter that is given
-    double saturation_ancillary(int param, int Q, int given, double value);
+    //double saturation_ancillary(int param, int Q, int given, double value);
 
     // ----------------------------------------
     // Transport properties
diff --git a/src/Backends/Helmholtz/FlashRoutines.cpp b/src/Backends/Helmholtz/FlashRoutines.cpp
index e1fa32ae..60f37698 100644
--- a/src/Backends/Helmholtz/FlashRoutines.cpp
+++ b/src/Backends/Helmholtz/FlashRoutines.cpp
@@ -140,7 +140,7 @@ void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS)
             pmin_sat = std::max(pmin_satL, pmin_satV);
             
             // Check limits
-            if (!is_in_closed_range(pmin_sat, pmax_sat, static_cast<long double>(HEOS._p))){
+            if (!is_in_closed_range(pmin_sat*0.999999, pmax_sat*1.000001, static_cast<long double>(HEOS._p))){
                 throw ValueError(format("Pressure to PQ_flash [%6g Pa] must be in range [%8g Pa, %8g Pa]",HEOS._p, pmin_sat, pmax_sat));
             }
             // ------------------
diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp
index 94978749..c3e51aa0 100644
--- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp
+++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp
@@ -131,7 +131,7 @@ void HelmholtzEOSMixtureBackend::update_states(void)
     // Clear again just to be sure
     clear();
 }
-const CoolProp::SimpleState HelmholtzEOSMixtureBackend::calc_state(const std::string &state)
+const CoolProp::SimpleState & HelmholtzEOSMixtureBackend::calc_state(const std::string &state)
 {
     if (is_pure_or_pseudopure)
     {
@@ -738,7 +738,7 @@ void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int ot
         }
     }
     // Check between triple point pressure and psat_max
-    else if (_p > components[0]->pEOS->ptriple && _p < _crit.p)
+    else if (_p >= components[0]->pEOS->ptriple*0.9999 && _p <= _crit.p)
     {
         // First try the ancillaries, use them to determine the state if you can
         
@@ -911,10 +911,13 @@ void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int ot
         _rhomolar = 1/(_Q/HEOS.SatV->rhomolar() + (1-_Q)/HEOS.SatL->rhomolar());
         return;
     }
-    else if (_p < components[0]->pEOS->ptriple)
+    else if (_p < components[0]->pEOS->ptriple*0.9999)
     {
         throw NotImplementedError(format("for now, we don't support p [%g Pa] below ptriple [%g Pa]",_p, components[0]->pEOS->ptriple));
     }
+    else{
+        throw ValueError(format("The pressure [%g Pa] cannot be used in p_phase_determination",_p));
+    }
 }
 void HelmholtzEOSMixtureBackend::T_phase_determination_pure_or_pseudopure(int other, long double value)
 {
diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h
index 2f3a9538..a78fa65e 100644
--- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h
+++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h
@@ -49,7 +49,7 @@ public:
 	long double calc_melting_line(int param, int given, long double value);
 	int calc_phase(void){return _phase;};
     
-    const CoolProp::SimpleState calc_state(const std::string &state);
+    const CoolProp::SimpleState &calc_state(const std::string &state);
 
     const std::vector<CoolPropFluid*> &get_components(){return components;};
     std::vector<long double> &get_K(){return K;};