diff --git a/include/teqp/algorithms/VLE.hpp b/include/teqp/algorithms/VLE.hpp
index de62203304207df2c941006b6f21dc3eb458ec06..3a17d69c6e90fc2ba4f71db5ddabd146c651442f 100644
--- a/include/teqp/algorithms/VLE.hpp
+++ b/include/teqp/algorithms/VLE.hpp
@@ -168,6 +168,36 @@ auto pure_VLE_T(const Model& model, Scalar T, Scalar rhoL, Scalar rhoV, int maxi
     return do_pure_VLE_T(res, rhoL, rhoV, maxiter);
 }
 
+/***
+ * \brief Calculate the derivative of vapor pressure with respect to temperature
+ * \param model The model to operate on
+ * \param T Temperature
+ * \param rhoL Liquid density
+ * \param rhoV Vapor density
+ *
+ *  Based upon
+ *  \f[
+ * \frac{dp_{\sigma}}{dT} = \frac{h''-h'}{T(v''-v')} = \frac{s''-s'}{v''-v'}
+ *  \f]
+ *  where the \f$h''-h'\f$ is given by the difference in residual enthalpy \f$h''-h' = h^r''-h^r'\f$ because the ideal-gas parts cancel
+ */
+template<typename Model, typename Scalar, ADBackends backend = ADBackends::autodiff>
+auto dpsatdT_pure(const Model& model, Scalar T, Scalar rhoL, Scalar rhoV) {
+    
+    auto molefrac = (Eigen::ArrayXd(1) << 1.0).finished();
+    
+    using tdx = TDXDerivatives<decltype(model), double, decltype(molefrac)>;
+    using iso = IsochoricDerivatives<decltype(model), double, decltype(molefrac)>;
+    
+    auto R = model.R(molefrac);
+    
+    auto hrVLERTV = tdx::get_Ar01(model, T, rhoV, molefrac) + tdx::get_Ar10(model, T, rhoV, molefrac);
+    auto hrVLERTL = tdx::get_Ar01(model, T, rhoL, molefrac) + tdx::get_Ar10(model, T, rhoL, molefrac);
+    auto deltahr_over_T = R*(hrVLERTV-hrVLERTL);
+    auto dpsatdT = deltahr_over_T/(1/rhoV-1/rhoL); // From Clausius-Clapeyron; dp/dT = Deltas/Deltav = Deltah/(T*Deltav); Delta=V-L
+    return dpsatdT;
+}
+
 /***
 * \brief Do a vapor-liquid phase equilibrium problem for a mixture (binary only for now) with mole fractions specified in the liquid phase
 * \param model The model to operate on
diff --git a/src/tests/catch_test_cubics.cxx b/src/tests/catch_test_cubics.cxx
index 6141d2745b1b309d6f5d53e663b5bda0b136d9a4..531cca9dd3c9e1299db15f72fd595fb6c78a4b4a 100644
--- a/src/tests/catch_test_cubics.cxx
+++ b/src/tests/catch_test_cubics.cxx
@@ -80,6 +80,46 @@ TEST_CASE("Check calling superancillary curves", "[cubic][superanc]")
     }
 }
 
+TEST_CASE("Check orthobaric density derivatives for pure fluid", "[cubic][superanc]")
+{
+    std::valarray<double> Tc_K = { 150.687 };
+    std::valarray<double> pc_Pa = { 4863000.0};
+    std::valarray<double> acentric = { 0.0};
+    
+    double T = 130.0, dT = 0.001;
+    auto molefrac = (Eigen::ArrayXd(1) << 1.0).finished();
+    
+    auto model = canonical_PR(Tc_K, pc_Pa, acentric);
+    using tdx = TDXDerivatives<decltype(model)>;
+    using iso = IsochoricDerivatives<decltype(model)>;
+    
+    auto R = model.R(molefrac);
+    auto [rhoL, rhoV] = model.superanc_rhoLV(T);
+    CHECK(rhoL > rhoV);
+    
+    // Finite difference test
+    auto [rhoLp, rhoVp] = model.superanc_rhoLV(T+dT);
+    auto [rhoLm, rhoVm] = model.superanc_rhoLV(T-dT);
+    auto pLp = rhoLp*R*(T+dT) + iso::get_pr(model, T+dT, rhoLp*molefrac);
+    auto pLm = rhoLm*R*(T-dT) + iso::get_pr(model, T-dT, rhoLm*molefrac);
+    
+    // Exact solution for density derivative
+    // Change in enthalpy (Deltah) is equal to change in residual enthalpy (Deltahr) because ideal parts cancel
+    auto hrVLERTV = tdx::get_Ar01(model, T, rhoV, molefrac) + tdx::get_Ar10(model, T, rhoV, molefrac);
+    auto hrVLERTL = tdx::get_Ar01(model, T, rhoL, molefrac) + tdx::get_Ar10(model, T, rhoL, molefrac);
+    auto deltahr_over_T = R*(hrVLERTV-hrVLERTL);
+    auto dpsatdT = deltahr_over_T/(1/rhoV-1/rhoL); // From Clausius-Clapeyron; dp/dT = Deltas/Deltav = Deltah/(T*Deltav); Delta=V-L
+    
+    auto dpsatdT_routine = dpsatdT_pure(model, T, rhoL, rhoV);
+    
+    CHECK(dpsatdT == Approx((pLp - pLm)/(2*dT)));
+    CHECK(dpsatdT_routine == Approx((pLp - pLm)/(2*dT)));
+    
+//    CHECK(drhovecdTL(0) == Approx((rhoLp-rhoLm)/(2*dT)));
+//    auto drhoLdT =
+    
+}
+
 TEST_CASE("Check manual integration of subcritical VLE isotherm for binary mixture", "[cubic][isochoric][traceisotherm]")
 {
     using namespace boost::numeric::odeint;
@@ -389,4 +429,4 @@ TEST_CASE("Check manual integration of subcritical VLE isobar for binary mixture
 
         std::ofstream file("isoP.json"); file << J;
     }
-}
\ No newline at end of file
+}