diff --git a/include/teqp/cpp/teqpcpp.hpp b/include/teqp/cpp/teqpcpp.hpp
index 334bb86c6ce07f04c15d2ecc1a050861da93d0af..8c3d4353cc3deb0c8a69c35052c845311eafc9cd 100644
--- a/include/teqp/cpp/teqpcpp.hpp
+++ b/include/teqp/cpp/teqpcpp.hpp
@@ -120,6 +120,7 @@ namespace teqp {
             
             virtual std::tuple<double, double> solve_pure_critical(const double T, const double rho, const std::optional<nlohmann::json>&) const = 0;
             virtual std::tuple<double, double> extrapolate_from_critical(const double Tc, const double rhoc, const double Tgiven) const = 0;
+            virtual std::tuple<EArrayd, EMatrixd> get_pure_critical_conditions_Jacobian(const double T, const double rho, int alternative_pure_index=-1, int alternative_length=2) const  = 0;
             virtual EArray2 pure_VLE_T(const double T, const double rhoL, const double rhoV, int maxiter) const = 0;
             
             virtual std::tuple<EArrayd, EArrayd> get_drhovecdp_Tsat(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const = 0;
diff --git a/interface/CPP/teqpcpp.cpp b/interface/CPP/teqpcpp.cpp
index aeb5b343c0920e23378a92172409bcb93e7e7757..9be1abaf39f37905c3fec2274f11aef376ba6199 100644
--- a/interface/CPP/teqpcpp.cpp
+++ b/interface/CPP/teqpcpp.cpp
@@ -85,6 +85,11 @@ namespace teqp {
                     return teqp::solve_pure_critical(model, T, rho, flags.value());
                 }, m_model);
             }
+            std::tuple<EArrayd, EMatrixd> get_pure_critical_conditions_Jacobian(const double T, const double rho, int alternative_pure_index, int alternative_length) const override {
+                return std::visit([&](const auto& model) {
+                    return teqp::get_pure_critical_conditions_Jacobian(model, T, rho, alternative_pure_index, alternative_length);
+                }, m_model);
+            }
             std::tuple<double, double> extrapolate_from_critical(const double Tc, const double rhoc, const double Tnew) const override {
                 return std::visit([&](const auto& model) {
                     auto mat = teqp::extrapolate_from_critical(model, Tc, rhoc, Tnew);
diff --git a/interface/pybind11_wrapper.cpp b/interface/pybind11_wrapper.cpp
index 99a637ac754a784a9083ac302cc41e3f7545a18c..017cc23fdb5811cc0c601168da43ce99a0506107 100644
--- a/interface/pybind11_wrapper.cpp
+++ b/interface/pybind11_wrapper.cpp
@@ -221,7 +221,7 @@ void init_teqp(py::module& m) {
         .def("get_deriv_mat2", &am::get_deriv_mat2, "T"_a, "rho"_a, "molefrac"_a.noconvert())
     
         // Routines related to pure fluid critical point calculation
-//    .def("get_pure_critical_conditions_Jacobian", &am::get_pure_critical_conditions_Jacobian, "T"_a, "rho"_a, py::arg_v("alternative_pure_index", -1), py::arg_v("alternative_length", 2))
+        .def("get_pure_critical_conditions_Jacobian", &am::get_pure_critical_conditions_Jacobian, "T"_a, "rho"_a, py::arg_v("alternative_pure_index", -1), py::arg_v("alternative_length", 2))
         .def("solve_pure_critical", &am::solve_pure_critical, "T"_a, "rho"_a, py::arg_v("flags", std::nullopt, "None"))
         .def("extrapolate_from_critical", &am::extrapolate_from_critical, "Tc"_a, "rhoc"_a, "T"_a)