diff --git a/CMakeLists.txt b/CMakeLists.txt
index bec518810391bc689e10fde8cae62c75b102658f..72594c0e71089eae8dfcc52f12efb7b3c18531ce 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -26,6 +26,7 @@ add_executable(trace "${CMAKE_CURRENT_SOURCE_DIR}/src/trace_critical.cpp")
 if (MSVC)
     target_sources(trace PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/externals/Eigen/debug/msvc/eigen.natvis")
 endif()
+target_link_libraries(trace autodiff)
 
 add_executable(multifluid "${CMAKE_CURRENT_SOURCE_DIR}/src/multifluid.cpp")
 if (MSVC)
diff --git a/include/teqp/derivs.hpp b/include/teqp/derivs.hpp
index 3e3417107389ec9d2ef11d43612518e7edaba003..8d106aadcf9d1e8dbce4aa8fb5995ed9fe7ca750 100644
--- a/include/teqp/derivs.hpp
+++ b/include/teqp/derivs.hpp
@@ -6,6 +6,11 @@
 
 #include "MultiComplex/MultiComplex.hpp"
 
+// autodiff include
+#include <autodiff/forward/dual.hpp>
+#include <autodiff/forward/dual/eigen.hpp>
+using namespace autodiff;
+
 template <typename TType, typename ContainerType, typename FuncType>
 typename std::enable_if<is_container<ContainerType>::value, typename ContainerType::value_type>::type
 caller(const FuncType& f, TType T, const ContainerType& rho) {
@@ -100,7 +105,7 @@ template <> void setval<std::valarray<std::valarray<double>>, std::size_t, doubl
 
 /***
 * \brief Calculate the Hessian of Psir = ar*rho w.r.t. the molar concentrations
-* 
+*
 * Requires the use of multicomplex derivatives to calculate second partial derivatives
 */
 template<typename Model, typename TType, typename RhoType>
@@ -108,6 +113,26 @@ auto build_Psir_Hessian(const Model& model, const TType T, const RhoType& rho) {
     // Double derivatives in each component's concentration
     // N^N matrix (symmetric)
 
+    dual2nd u; // the output scalar u = f(x), evaluated together with Hessian below
+    VectorXdual2nd g;
+    VectorXdual2nd rhovecc(2); rhovecc << rho[0], rho[1];
+    auto hfunc = [&model, &T](const VectorXdual2nd& rho_) {
+        auto rhotot_ = std::accumulate(std::begin(rho_), std::end(rho_), (decltype(rho_[0]))0.0);
+        return eval(model.alphar(T, rho_) * model.R * T * rhotot_);
+    };
+    return autodiff::hessian(hfunc, wrt(rhovecc), at(rhovecc), u, g); // evaluate the function value u, its gradient, and its Hessian matrix H
+}
+
+/***
+* \brief Calculate the Hessian of Psir = ar*rho w.r.t. the molar concentrations
+* 
+* Requires the use of multicomplex derivatives to calculate second partial derivatives
+*/
+template<typename Model, typename TType, typename RhoType>
+auto build_Psir_Hessian_mcx(const Model& model, const TType T, const RhoType& rho) {
+    // Double derivatives in each component's concentration
+    // N^N matrix (symmetric)
+
     // Lambda function for getting Psir with multicomplex concentrations
     auto func = [&model, &T](const std::vector<MultiComplex<double>>& rhovec) {
         auto N = rhovec.size();
diff --git a/src/test_autodiff.cpp b/src/test_autodiff.cpp
index de0718e9386ea684ee2675f13d168d2a78cf1e45..7bf7057a68500e976013bf33ad82435e76bfa452 100644
--- a/src/test_autodiff.cpp
+++ b/src/test_autodiff.cpp
@@ -7,7 +7,6 @@
 
 #include "MultiComplex/MultiComplex.hpp"
 
-
 // autodiff include
 #include <autodiff/forward/dual.hpp>
 #include <autodiff/forward/dual/eigen.hpp>
diff --git a/src/trace_critical.cpp b/src/trace_critical.cpp
index 268c0222ad99c54c1521f27852f9062039047c28..284cb1d9166bce19185f2bc388a48c0af7a64760 100644
--- a/src/trace_critical.cpp
+++ b/src/trace_critical.cpp
@@ -14,17 +14,20 @@ void trace() {
     auto T = Tc_K[0];
     const auto dT = 1;
     std::valarray<double> rhovec = { rhoc0, 0.0 };
+    auto tic0 = std::chrono::steady_clock::now();
     for (auto iter = 0; iter < 1000; ++iter){
         auto drhovecdT = get_drhovec_dT_crit(vdW, T, rhovec);
         rhovec += drhovecdT*dT;
         T += dT;
         int rr = 0;
         auto z0 = rhovec[0] / rhovec.sum();
-        std::cout << z0 <<  "," << rhovec[0] << "," << T << "," << get_splus(vdW, T, rhovec) << std::endl;
+        //std::cout << z0 <<  "," << rhovec[0] << "," << T << "," << get_splus(vdW, T, rhovec) << std::endl;
         if (z0 < 0) {
             break;
         }
     }
+    auto tic1 = std::chrono::steady_clock::now();
+    std::cout << std::chrono::duration<double>(tic1 - tic0).count() << " s to trace" << std::endl;
     int rr = 0;
 }
 int main() {