diff --git a/include/teqp/models/squarewell.hpp b/include/teqp/models/squarewell.hpp
index 76d53bc42a73656678a24213a9694c773df56997..c6c1c9fd072e0f57d28914ae1646b219004ba9e9 100644
--- a/include/teqp/models/squarewell.hpp
+++ b/include/teqp/models/squarewell.hpp
@@ -7,6 +7,8 @@
 namespace teqp{
 namespace squarewell{
 
+#include "teqp/types.hpp"
+
 /**
  Rodolfo Espíndola-Heredia, Fernando del Río and Anatol Malijevsky
  Optimized equation of the state of the
@@ -122,7 +124,7 @@ private:
             den += thetai[n]*pow(lambda_, n-4);
         }
         den = 1.0 + rhostar*den;
-        return num/den;
+        return forceeval(num/den);
     }
     
     template<typename RhoType>
diff --git a/interface/pybind11_wrapper.cpp b/interface/pybind11_wrapper.cpp
index a1d84da06dedd39c13125d8292cead1b4b975dea..992308a5dcfdf509cc2c1c0e6649cb59f388f446 100644
--- a/interface/pybind11_wrapper.cpp
+++ b/interface/pybind11_wrapper.cpp
@@ -16,6 +16,7 @@ void add_CPA(py::module& m);
 void add_multifluid(py::module& m);
 void add_multifluid_mutant(py::module& m);
 void add_cubics(py::module& m);
+void add_model_potentials(py::module& m);
 
 template<typename Model, int iT, int iD, typename Class>
 void add_ig_deriv_impl(Class& cls) {
@@ -169,6 +170,7 @@ void init_teqp(py::module& m) {
     add_multifluid(m);
     add_multifluid_mutant(m);
     add_cubics(m);
+    add_model_potentials(m);
 
     call_method_factory(m, "get_Ar00iso");
     call_method_factory(m, "get_Ar10iso");
diff --git a/src/tests/catch_test_SW.cxx b/src/tests/catch_test_SW.cxx
index 459c63db59cca145ed755dd22c4900461ca1d88f..c087b4df42b08b02eb2ce439780b2459e96b629b 100644
--- a/src/tests/catch_test_SW.cxx
+++ b/src/tests/catch_test_SW.cxx
@@ -3,6 +3,8 @@
 
 using Catch::Approx;
 
+#include "teqp/algorithms/critical_pure.hpp"
+
 #include "teqp/models/squarewell.hpp"
 #include "teqp/derivs.hpp"
 #include <Eigen/Dense>
@@ -13,13 +15,16 @@ TEST_CASE("simple evaluation of s^+ for square well EOS", "[squarewell]")
 {
     auto model = squarewell::EspindolaHeredia2009(1.5);
     std::valarray<double> z = {1.0};
-    auto ar = model.alphar(1.3144366466267958, 0.2862336473147125, z);
+    auto Tstar=1.3144366466267958, rhostar = 0.2862336473147125;
+    auto ar = model.alphar(Tstar, rhostar, z);
     
     using id = IsochoricDerivatives<decltype(model)>;
-    auto rhovec = (Eigen::ArrayXd(1) << 0.2862336473147125).finished();
-    auto splus = id::get_splus(model, 1.3144366466267958, rhovec);
+    auto rhovec = (Eigen::ArrayXd(1) << rhostar).finished();
+    auto splus = id::get_splus(model, Tstar, rhovec);
     auto alphar_target = -0.8061758248466638;
     auto splus_target = 1.0031288550954747;
     CHECK(splus_target == Approx(splus));
-
+    
+    auto cr = teqp::get_pure_critical_conditions_Jacobian(model, Tstar, rhostar);
+    CHECK(std::isfinite(std::get<1>(cr)(0,0)));
 }