diff --git a/src/tests/catch_test_PCSAFT.cxx b/src/tests/catch_test_PCSAFT.cxx
index 81ba7f20..8223d96c 100644
--- a/src/tests/catch_test_PCSAFT.cxx
+++ b/src/tests/catch_test_PCSAFT.cxx
@@ -4,8 +4,11 @@
 using Catch::Approx;
 
 #include "teqp/core.hpp"
+#include "teqp/derivs.hpp"
 #include "teqp/models/pcsaft.hpp"
+#include "teqp/finite_derivs.hpp"
 using namespace teqp::PCSAFT;
+using namespace teqp;
 
 TEST_CASE("Single alphar check value", "[PCSAFT]")
 {
@@ -17,6 +20,95 @@ TEST_CASE("Single alphar check value", "[PCSAFT]")
     CHECK(tdx::get_Ar00(model, T, Dmolar, z) == Approx(-0.032400020930842724));
 }
 
+TEST_CASE("Check 0n derivatives", "[PCSAFT]")
+{
+    std::vector<std::string> names = { "Methane", "Ethane" };
+    auto model = PCSAFTMixture(names);
+
+    const double T = 100.0;
+    const double rho = 126.1856883066021;
+    const auto rhovec = (Eigen::ArrayXd(2) << rho, 0).finished();
+    const auto molefrac = rhovec / rhovec.sum();
+
+    using my_float_type = boost::multiprecision::number<boost::multiprecision::cpp_bin_float<100U>>;
+    my_float_type D = rho, h = pow(my_float_type(10.0), -10);
+    auto fD = [&](const auto& x) { return model.alphar(T, x, molefrac); };
+
+    using tdx = TDXDerivatives<decltype(model)>;
+    auto Ar02 = tdx::get_Ar02(model, T, rho, molefrac);
+    auto Ar02n = tdx::get_Ar0n<2>(model, T, rho, molefrac)[2];
+    auto Ar02mp = static_cast<double>((D * D) * centered_diff<2, 4>(fD, D, h));
+    auto Ar02mcx = tdx::get_Ar0n<2, ADBackends::multicomplex>(model, T, rho, molefrac)[2];
+    CAPTURE(Ar02);
+    CAPTURE(Ar02n);
+    CAPTURE(Ar02mp);
+    CAPTURE(Ar02mcx);
+    CHECK(std::abs(Ar02 - Ar02n) < 1e-13);
+    CHECK(std::abs(Ar02 - Ar02mp) < 1e-13);
+    CHECK(std::abs(Ar02 - Ar02mcx) < 1e-13);
+
+    auto Ar01 = tdx::get_Ar01(model, T, rho, molefrac);
+    auto Ar01n = tdx::get_Ar0n<1>(model, T, rho, molefrac)[1];
+    auto Ar01mcx = tdx::get_Ar0n<1, ADBackends::multicomplex>(model, T, rho, molefrac)[1];
+    auto Ar01csd = tdx::get_Ar01<ADBackends::complex_step>(model, T, rho, molefrac);
+    auto Ar01mp = static_cast<double>(D * centered_diff<1, 4>(fD, D, h));
+    CAPTURE(Ar01);
+    CAPTURE(Ar01n);
+    CAPTURE(Ar01mp);
+    CAPTURE(Ar01mcx);
+    CAPTURE(Ar01csd);
+    CHECK(std::abs(Ar01 - Ar01n) < 1e-13);
+    CHECK(std::abs(Ar01 - Ar01mp) < 1e-13);
+    CHECK(std::abs(Ar01 - Ar01mcx) < 1e-13);
+    CHECK(std::abs(Ar01 - Ar01csd) < 1e-13);
+
+    auto Ar03 = tdx::get_Arxy<0, 3, ADBackends::autodiff>(model, T, rho, molefrac);
+    auto Ar03n = tdx::get_Ar0n<3>(model, T, rho, molefrac)[3];
+    auto Ar03mp = static_cast<double>((D * D * D) * centered_diff<3, 4>(fD, D, h));
+    auto Ar03mcx = tdx::get_Ar0n<3, ADBackends::multicomplex>(model, T, rho, molefrac)[3];
+    CAPTURE(Ar03);
+    CAPTURE(Ar03n);
+    CAPTURE(Ar03mp);
+    CAPTURE(Ar03mcx);
+    CHECK(std::abs(Ar03 - Ar03n) < 1e-13);
+    CHECK(std::abs(Ar03 - Ar03mp) < 1e-13);
+    CHECK(std::abs(Ar03 - Ar03mcx) < 1e-13);
+
+    auto Ar04 = tdx::get_Arxy<0, 4, ADBackends::autodiff>(model, T, rho, molefrac);
+    auto Ar04n = tdx::get_Ar0n<4>(model, T, rho, molefrac)[4];
+    auto Ar04mp = static_cast<double>((D * D * D * D) * centered_diff<4, 4>(fD, D, h));
+    auto Ar04mcx = tdx::get_Ar0n<4, ADBackends::multicomplex>(model, T, rho, molefrac)[4];
+    CAPTURE(Ar04);
+    CAPTURE(Ar04n);
+    CAPTURE(Ar04mp);
+    CAPTURE(Ar04mcx);
+    CHECK(std::abs(Ar04 - Ar04n) < 1e-13);
+    CHECK(std::abs(Ar04 - Ar04mp) < 1e-13);
+    CHECK(std::abs(Ar04 - Ar04mcx) < 1e-13);
+}
+
+TEST_CASE("Check neff", "[virial]")
+{
+    double T = 298.15;
+    double rho = 3.0;
+    const Eigen::Array2d molefrac = { 0.5, 0.5 };
+    auto f = [&T, &rho, &molefrac](const auto& model) {
+        auto neff = TDXDerivatives<decltype(model)>::get_neff(model, T, rho, molefrac);
+        CAPTURE(neff);
+        CHECK(neff > 0);
+        CHECK(neff < 100);
+    };
+    // This quantity is undefined for the van der Waals EOS because Ar20 is always 0
+    //SECTION("vdW") {
+    //    f(build_simple());
+    //}
+    SECTION("PCSAFT") {
+        std::vector<std::string> names = { "Methane", "Ethane" };
+        f(PCSAFTMixture(names));
+    }
+}
+
+
 TEST_CASE("Check PCSAFT with kij", "[PCSAFT]")
 {
     std::vector<std::string> names = { "Methane", "Ethane" };
diff --git a/src/tests/catch_tests.cxx b/src/tests/catch_tests.cxx
index 3c09739c..34fb1185 100644
--- a/src/tests/catch_tests.cxx
+++ b/src/tests/catch_tests.cxx
@@ -4,8 +4,6 @@
 using Catch::Approx;
 
 #include "teqp/core.hpp"
-#include "teqp/models/pcsaft.hpp"
-using namespace teqp::PCSAFT;
 #include "teqp/models/cubicsuperancillary.hpp"
 #include "teqp/models/CPA.hpp"
 #include "teqp/models/vdW.hpp"
@@ -100,26 +98,6 @@ TEST_CASE("Check virial coefficients for vdW", "[virial]")
     }
 }
 
-TEST_CASE("Check neff", "[virial]")
-{
-    double T = 298.15;
-    double rho = 3.0;
-    const Eigen::Array2d molefrac = { 0.5, 0.5 };
-    auto f = [&T, &rho, &molefrac](const auto& model) {
-        auto neff = TDXDerivatives<decltype(model)>::get_neff(model, T, rho, molefrac);
-        CAPTURE(neff);
-        CHECK(neff > 0);
-        CHECK(neff < 100);
-    };
-    // This quantity is undefined for the van der Waals EOS because Ar20 is always 0
-    //SECTION("vdW") {
-    //    f(build_simple());
-    //}
-    SECTION("PCSAFT") {
-        std::vector<std::string> names = { "Methane", "Ethane" };
-        f(PCSAFTMixture(names));
-    }
-}
 
 TEST_CASE("Check Hessian of Psir", "[virial]")
 {
@@ -177,73 +155,6 @@ TEST_CASE("Check p four ways for vdW", "[virial][p]")
     CHECK(std::abs(p_ar0n - pexact) / pexact < 1e-8);
 }
 
-TEST_CASE("Check 0n derivatives", "[PCSAFT]")
-{
-    std::vector<std::string> names = { "Methane", "Ethane" };
-    auto model = PCSAFTMixture(names);
-
-    const double T = 100.0;
-    const double rho = 126.1856883066021; 
-    const auto rhovec = (Eigen::ArrayXd(2) << rho, 0).finished();
-    const auto molefrac = rhovec / rhovec.sum();
-
-    using my_float_type = boost::multiprecision::number<boost::multiprecision::cpp_bin_float<100U>>;
-    my_float_type D = rho, h = pow(my_float_type(10.0), -10);
-    auto fD = [&](const auto& x) { return model.alphar(T, x, molefrac); };
-    
-    using tdx = TDXDerivatives<decltype(model)>;
-    auto Ar02 = tdx::get_Ar02(model, T, rho, molefrac);
-    auto Ar02n = tdx::get_Ar0n<2>(model, T, rho, molefrac)[2];
-    auto Ar02mp = static_cast<double>((D*D)*centered_diff<2, 4>(fD, D, h));
-    auto Ar02mcx = tdx::get_Ar0n<2,ADBackends::multicomplex>(model, T, rho, molefrac)[2];
-    CAPTURE(Ar02);
-    CAPTURE(Ar02n);
-    CAPTURE(Ar02mp);
-    CAPTURE(Ar02mcx);
-    CHECK(std::abs(Ar02 - Ar02n) < 1e-13);
-    CHECK(std::abs(Ar02 - Ar02mp) < 1e-13);
-    CHECK(std::abs(Ar02 - Ar02mcx) < 1e-13);
-
-    auto Ar01 = tdx::get_Ar01(model, T, rho, molefrac);
-    auto Ar01n = tdx::get_Ar0n<1>(model, T, rho, molefrac)[1];
-    auto Ar01mcx = tdx::get_Ar0n<1,ADBackends::multicomplex>(model, T, rho, molefrac)[1];
-    auto Ar01csd = tdx::get_Ar01<ADBackends::complex_step>(model, T, rho, molefrac);
-    auto Ar01mp = static_cast<double>(D*centered_diff<1, 4>(fD, D, h));
-    CAPTURE(Ar01);
-    CAPTURE(Ar01n);
-    CAPTURE(Ar01mp);
-    CAPTURE(Ar01mcx);
-    CAPTURE(Ar01csd);
-    CHECK(std::abs(Ar01 - Ar01n) < 1e-13);
-    CHECK(std::abs(Ar01 - Ar01mp) < 1e-13);
-    CHECK(std::abs(Ar01 - Ar01mcx) < 1e-13);
-    CHECK(std::abs(Ar01 - Ar01csd) < 1e-13);
-
-    auto Ar03 = tdx::get_Arxy<0,3,ADBackends::autodiff>(model, T, rho, molefrac);
-    auto Ar03n = tdx::get_Ar0n<3>(model, T, rho, molefrac)[3];
-    auto Ar03mp = static_cast<double>((D*D*D)*centered_diff<3, 4>(fD, D, h));
-    auto Ar03mcx = tdx::get_Ar0n<3,ADBackends::multicomplex>(model, T, rho, molefrac)[3];
-    CAPTURE(Ar03);
-    CAPTURE(Ar03n);
-    CAPTURE(Ar03mp);
-    CAPTURE(Ar03mcx);
-    CHECK(std::abs(Ar03 - Ar03n) < 1e-13);
-    CHECK(std::abs(Ar03 - Ar03mp) < 1e-13);
-    CHECK(std::abs(Ar03 - Ar03mcx) < 1e-13);
-
-    auto Ar04 = tdx::get_Arxy<0,4,ADBackends::autodiff>(model, T, rho, molefrac);
-    auto Ar04n = tdx::get_Ar0n<4>(model, T, rho, molefrac)[4];
-    auto Ar04mp = static_cast<double>((D*D*D*D)*centered_diff<4, 4>(fD, D, h));
-    auto Ar04mcx = tdx::get_Ar0n<4,ADBackends::multicomplex>(model, T, rho, molefrac)[4];
-    CAPTURE(Ar04);
-    CAPTURE(Ar04n);
-    CAPTURE(Ar04mp);
-    CAPTURE(Ar04mcx);
-    CHECK(std::abs(Ar04 - Ar04n) < 1e-13);
-    CHECK(std::abs(Ar04 - Ar04mp) < 1e-13);
-    CHECK(std::abs(Ar04 - Ar04mcx) < 1e-13);
-}
-
 TEST_CASE("Test all vdW derivatives good to numerical precision", "[vdW]") 
 {
     using my_float_type = boost::multiprecision::number<boost::multiprecision::cpp_bin_float<100U>>;