diff --git a/include/teqp/derivs.hpp b/include/teqp/derivs.hpp
index 218f66905cb789e4096a93b90327ba05ddd69c4e..9fb166d415070385dca4f6a05c0a19f721158f21 100644
--- a/include/teqp/derivs.hpp
+++ b/include/teqp/derivs.hpp
@@ -97,17 +97,9 @@ typename ContainerType::value_type get_Ar01(const Model& model, const TType T, c
 
 template <typename Model, typename TType, typename ContainerType>
 typename ContainerType::value_type get_B2vir(const Model& model, const TType T, const ContainerType& molefrac) {
-    auto rhovec = 0.0*molefrac;
-    double max_mole_frac = *std::max_element(std::begin(molefrac), std::end(molefrac));
-    if (max_mole_frac != 1.0) {
-        std::cout << "The value for B_2 is incorrect when not pure!" << std::endl;
-    }
-    decltype(molefrac[0] * T) B2 = 0.0;
-    for (auto i = 0; i < rhovec.size(); ++i) {
-        auto dalphar_drhoi__constTrhoj = derivrhoi([&model](const auto& T, const auto& rhovec) { return model.alphar(T, rhovec); }, T, rhovec, i);
-        B2 += molefrac[i] * dalphar_drhoi__constTrhoj;
-    }
-    
+    double h = 1e-100;
+    // B_2 = lim_rho\to 0 dalphar/drho|T,z
+    auto B2 = model.alphar(T, std::complex<double>(0, h), molefrac).imag()/h;
     return B2;
 }
 
diff --git a/include/teqp/models/multifluid.hpp b/include/teqp/models/multifluid.hpp
index eeafd62c9f1bfa427e42944dcebbd5751e1a4130..ee469d550879d2a9c4bfbad7078db84f1a1f2f97 100644
--- a/include/teqp/models/multifluid.hpp
+++ b/include/teqp/models/multifluid.hpp
@@ -84,11 +84,18 @@ public:
     {
         RhoType::value_type rhotot_ = (rhotot.has_value()) ? rhotot.value() : std::accumulate(std::begin(rhovec), std::end(rhovec), (decltype(rhovec[0]))0.0);
         auto molefrac = rhovec / rhotot_;
+        return alphar(T, rhotot_, molefrac);
+    }
+
+    template<typename TType, typename RhoType, typename MoleFracType>
+    auto alphar(TType T,
+        const RhoType rho,
+        const MoleFracType& molefrac) const
+    {
         auto Tred = redfunc.get_Tr(molefrac);
         auto rhored = redfunc.get_rhor(molefrac);
-        auto delta = rhotot_ / rhored;
+        auto delta = rho / rhored;
         auto tau = Tred / T;
-        using resulttype = decltype(T*rhovec[0]);
         return corr.alphar(tau, delta, molefrac) + dep.alphar(tau, delta, molefrac);
     }
 };