diff --git a/include/teqp/algorithms/VLE.hpp b/include/teqp/algorithms/VLE.hpp
index 477d7f946ec7ec72b1933e8a42158ea9876125cb..8b8db6c3471e185f93757f71952702205815ca62 100644
--- a/include/teqp/algorithms/VLE.hpp
+++ b/include/teqp/algorithms/VLE.hpp
@@ -265,7 +265,7 @@ auto mix_VLE_Tx(const Model& model, Scalar T, const Vector& rhovecL0, const Vect
     return std::make_tuple(return_code, rhovecLfinal, rhovecVfinal);
 }
 
-struct MixVLETPFlags {
+struct MixVLETpFlags {
     double atol = 1e-10,
         reltol = 1e-10,
         axtol = 1e-10,
@@ -286,7 +286,7 @@ struct MixVLETPFlags {
 * \param flags Flags controlling the iteration and stopping conditions
 */
 template<typename Model, typename Scalar, typename Vector>
-auto mix_VLE_TP(const Model& model, Scalar T, Scalar pgiven, const Vector& rhovecL0, const Vector& rhovecV0, const MixVLETPFlags& flags = {}) {
+auto mix_VLE_Tp(const Model& model, Scalar T, Scalar pgiven, const Vector& rhovecL0, const Vector& rhovecV0, const MixVLETpFlags& flags = {}) {
 
     const Eigen::Index N = rhovecL0.size();
     auto lengths = (Eigen::ArrayXi(2) << rhovecL0.size(), rhovecV0.size()).finished();
@@ -396,7 +396,7 @@ auto mix_VLE_TP(const Model& model, Scalar T, Scalar pgiven, const Vector& rhove
     return std::make_tuple(return_code, message, rhovecLfinal, rhovecVfinal);
 }
 
-struct MixVLEPxFlags {
+struct MixVLEpxFlags {
     double atol = 1e-10,
         reltol = 1e-10,
         axtol = 1e-10,
@@ -416,7 +416,7 @@ struct MixVLEPxFlags {
 * \param flags Additional flags
 */
 template<typename Model, typename Scalar, typename Vector>
-auto mixture_VLE_px(const Model& model, Scalar p_spec, const Vector& xmolar_spec, Scalar T0, const Vector& rhovecL0, const Vector& rhovecV0, const MixVLEPxFlags& flags = {}) {
+auto mixture_VLE_px(const Model& model, Scalar p_spec, const Vector& xmolar_spec, Scalar T0, const Vector& rhovecL0, const Vector& rhovecV0, const MixVLEpxFlags& flags = {}) {
 
     const Eigen::Index N = rhovecL0.size();
     auto lengths = (Eigen::ArrayXi(3) << rhovecL0.size(), rhovecV0.size(), xmolar_spec.size()).finished();
diff --git a/interface/pybind11_wrapper.cpp b/interface/pybind11_wrapper.cpp
index 3b093221abbfc5e18dc5e26e1efc44bd29b75658..ec35d4c5daaef2e83f2e5167e390c066f4626e34 100644
--- a/interface/pybind11_wrapper.cpp
+++ b/interface/pybind11_wrapper.cpp
@@ -57,53 +57,6 @@ inline auto call_method_factory(py::module &m, const std::string& attribute) {
 
 /// Instantiate "instances" of models (really wrapped Python versions of the models), and then attach all derivative methods
 void init_teqp(py::module& m) {
-    add_vdW(m);
-    add_PCSAFT(m);
-    add_CPA(m);
-    add_multifluid(m);
-    add_multifluid_mutant(m);
-    add_cubics(m);
-
-    call_method_factory(m, "get_Ar00iso");
-    call_method_factory(m, "get_Ar10iso");
-    call_method_factory(m, "get_Psiriso"),
-
-    call_method_factory(m, "get_splus");
-    call_method_factory(m, "get_pr");
-    call_method_factory(m, "get_B2vir");
-    call_method_factory(m, "get_B12vir");
-    
-    call_method_factory(m, "pure_VLE_T");
-    call_method_factory(m, "extrapolate_from_critical");
-
-    call_method_factory(m, "build_Psir_Hessian_autodiff");
-    call_method_factory(m, "build_Psi_Hessian_autodiff");
-    call_method_factory(m, "build_Psir_gradient_autodiff");
-    call_method_factory(m, "build_d2PsirdTdrhoi_autodiff");
-    call_method_factory(m, "get_chempotVLE_autodiff");
-    call_method_factory(m, "get_dchempotdT_autodiff");
-    call_method_factory(m, "get_fugacity_coefficients");
-    call_method_factory(m, "get_partial_molar_volumes");
-
-    call_method_factory(m, "trace_critical_arclength_binary");
-    call_method_factory(m, "get_criticality_conditions");
-    call_method_factory(m, "eigen_problem");
-    call_method_factory(m, "get_minimum_eigenvalue_Psi_Hessian");
-    call_method_factory(m, "get_drhovec_dT_crit");
-
-    call_method_factory(m, "get_pure_critical_conditions_Jacobian");
-    call_method_factory(m, "solve_pure_critical");
-    call_method_factory(m, "mix_VLE_Tx");
-    call_method_factory(m, "mixture_VLE_px");
-
-    call_method_factory(m, "get_drhovecdp_Tsat");
-    call_method_factory(m, "trace_VLE_isotherm_binary");
-    call_method_factory(m, "get_drhovecdT_psat");
-    call_method_factory(m, "trace_VLE_isobar_binary");
-    call_method_factory(m, "get_dpsat_dTsat_isopleth");
-
-    call_method_factory(m, "mix_VLLE_T");
-    call_method_factory(m, "find_VLLE_T_binary");
 
     // The options class for critical tracer, not tied to a particular model
     py::class_<TCABOptions>(m, "TCABOptions")
@@ -165,6 +118,24 @@ void init_teqp(py::module& m) {
         .def_readwrite("rho_trivial_threshold", &VLLEFinderOptions::rho_trivial_threshold)
         ;
 
+    py::class_<MixVLETpFlags>(m, "MixVLETpFlags")
+        .def(py::init<>())
+        .def_readwrite("atol", &MixVLETpFlags::atol)
+        .def_readwrite("reltol", &MixVLETpFlags::reltol)
+        .def_readwrite("relaxtoltol", &MixVLETpFlags::axtol)
+        .def_readwrite("relxtol", &MixVLETpFlags::relxtol)
+        .def_readwrite("maxiter", &MixVLETpFlags::maxiter)
+        ;
+
+    py::class_<MixVLEpxFlags>(m, "MixVLEpxFlags")
+        .def(py::init<>())
+        .def_readwrite("atol", &MixVLEpxFlags::atol)
+        .def_readwrite("reltol", &MixVLEpxFlags::reltol)
+        .def_readwrite("relaxtoltol", &MixVLEpxFlags::axtol)
+        .def_readwrite("relxtol", &MixVLEpxFlags::relxtol)
+        .def_readwrite("maxiter", &MixVLEpxFlags::maxiter)
+        ;
+
     py::enum_<VLE_return_code>(m, "VLE_return_code")
         .value("unset", VLE_return_code::unset)
         .value("xtol_satisfied", VLE_return_code::xtol_satisfied)
@@ -177,6 +148,54 @@ void init_teqp(py::module& m) {
     auto alphaig = py::class_<IdealHelmholtz>(m, "IdealHelmholtz").def(py::init<const nlohmann::json&>());
     add_ig_derivatives<IdealHelmholtz>(m, alphaig);
 
+    add_vdW(m);
+    add_PCSAFT(m);
+    add_CPA(m);
+    add_multifluid(m);
+    add_multifluid_mutant(m);
+    add_cubics(m);
+
+    call_method_factory(m, "get_Ar00iso");
+    call_method_factory(m, "get_Ar10iso");
+    call_method_factory(m, "get_Psiriso"),
+
+    call_method_factory(m, "get_splus");
+    call_method_factory(m, "get_pr");
+    call_method_factory(m, "get_B2vir");
+    call_method_factory(m, "get_B12vir");
+    
+    call_method_factory(m, "pure_VLE_T");
+    call_method_factory(m, "extrapolate_from_critical");
+
+    call_method_factory(m, "build_Psir_Hessian_autodiff");
+    call_method_factory(m, "build_Psi_Hessian_autodiff");
+    call_method_factory(m, "build_Psir_gradient_autodiff");
+    call_method_factory(m, "build_d2PsirdTdrhoi_autodiff");
+    call_method_factory(m, "get_chempotVLE_autodiff");
+    call_method_factory(m, "get_dchempotdT_autodiff");
+    call_method_factory(m, "get_fugacity_coefficients");
+    call_method_factory(m, "get_partial_molar_volumes");
+
+    call_method_factory(m, "trace_critical_arclength_binary");
+    call_method_factory(m, "get_criticality_conditions");
+    call_method_factory(m, "eigen_problem");
+    call_method_factory(m, "get_minimum_eigenvalue_Psi_Hessian");
+    call_method_factory(m, "get_drhovec_dT_crit");
+
+    call_method_factory(m, "get_pure_critical_conditions_Jacobian");
+    call_method_factory(m, "solve_pure_critical");
+    call_method_factory(m, "mix_VLE_Tx");
+    call_method_factory(m, "mixture_VLE_px");
+
+    call_method_factory(m, "get_drhovecdp_Tsat");
+    call_method_factory(m, "trace_VLE_isotherm_binary");
+    call_method_factory(m, "get_drhovecdT_psat");
+    call_method_factory(m, "trace_VLE_isobar_binary");
+    call_method_factory(m, "get_dpsat_dTsat_isopleth");
+
+    call_method_factory(m, "mix_VLLE_T");
+    call_method_factory(m, "find_VLLE_T_binary");
+
     // Some functions for timing overhead of interface
     m.def("___mysummer", [](const double &c, const Eigen::ArrayXd &x) { return c*x.sum(); });
     using RAX = Eigen::Ref<Eigen::ArrayXd>;
diff --git a/interface/pybind11_wrapper.hpp b/interface/pybind11_wrapper.hpp
index d0ba4f9a639701fdcb74414695611a63265df924..bbdc155037cb460197f3f2bdcd04eebbf532b371 100644
--- a/interface/pybind11_wrapper.hpp
+++ b/interface/pybind11_wrapper.hpp
@@ -78,9 +78,9 @@ void add_derivatives(py::module &m, Wrapper &cls) {
     
     cls.def("get_pure_critical_conditions_Jacobian", &get_pure_critical_conditions_Jacobian<Model, double, ADBackends::autodiff>, py::arg("T"), py::arg("rho"), py::arg_v("alternative_pure_index", -1), py::arg_v("alternative_length", 2));
     cls.def("solve_pure_critical", &solve_pure_critical<Model, double, ADBackends::autodiff>, py::arg("T"), py::arg("rho"), py::arg_v("flags", std::nullopt, "None"));
-    cls.def("mix_VLE_Tx", &mix_VLE_Tx<Model, double, Eigen::ArrayXd>);
-    cls.def("mix_VLE_TP", &mix_VLE_TP<Model, double, Eigen::ArrayXd>);
-    cls.def("mixture_VLE_px", &mixture_VLE_px<Model, double, Eigen::ArrayXd>, py::arg("p_spec"), py::arg("xmolar_spec").noconvert(), py::arg("T0"), py::arg("rhovecL0").noconvert(), py::arg("rhovecV0").noconvert(), py::arg_v("flags", std::nullopt, "None"));
+    cls.def("mix_VLE_Tx", &mix_VLE_Tx<Model, double, RAX>);
+    cls.def("mix_VLE_Tp", &mix_VLE_Tp<Model, double, RAX>, py::arg("T"), py::arg("p_given"), py::arg("rhovecL0").noconvert(), py::arg("rhovecV0").noconvert(), py::arg_v("flags", MixVLETpFlags{}, "None"));
+    cls.def("mixture_VLE_px", &mixture_VLE_px<Model, double, Eigen::ArrayXd>, py::arg("p_spec"), py::arg("xmolar_spec").noconvert(), py::arg("T0"), py::arg("rhovecL0").noconvert(), py::arg("rhovecV0").noconvert(), py::arg_v("flags", MixVLEpxFlags{}, "None"));
 
     cls.def("get_drhovecdp_Tsat", &get_drhovecdp_Tsat<Model, double, RAX>, py::arg("T"), py::arg("rhovecL").noconvert(), py::arg("rhovecV").noconvert());
     cls.def("trace_VLE_isotherm_binary", &trace_VLE_isotherm_binary<Model, double, Eigen::ArrayXd>, py::arg("T"), py::arg("rhovecL0").noconvert(), py::arg("rhovecV0").noconvert(), py::arg_v("options", std::nullopt, "None"));
diff --git a/src/tests/catch_test_cubics.cxx b/src/tests/catch_test_cubics.cxx
index d4e616e99a4a0e201275e23e8b41190cf7dad3c7..a2148d9dc8eca74e008d690970d3d02bb5774e04 100644
--- a/src/tests/catch_test_cubics.cxx
+++ b/src/tests/catch_test_cubics.cxx
@@ -157,7 +157,7 @@ TEST_CASE("Check manual integration of subcritical VLE isotherm for binary mixtu
                     auto [return_code, rhoL, rhoV] = mix_VLE_Tx(model, T, rhovecL, rhovecV, x, 1e-10, 1e-8, 1e-10, 1e-8, 10);
 
                     // And the other way around just to test the routine for TP solving
-                    auto [return_code2, msg, rhoL_, rhoV_] = mix_VLE_TP(model, T, p, rhovecL, rhovecV);
+                    auto [return_code2, msg, rhoL_, rhoV_] = mix_VLE_Tp(model, T, p, rhovecL, rhovecV);
                     int rr = 0;
 
                 }