diff --git a/.github/disabled_workflows/build_clangcl.yml b/.github/disabled_workflows/build_clangcl.yml
new file mode 100644
index 0000000000000000000000000000000000000000..59ebce4ee57a0deb7486200b47e565f0b1656096
--- /dev/null
+++ b/.github/disabled_workflows/build_clangcl.yml
@@ -0,0 +1,27 @@
+name: build and run with clang-cl
+on:
+  push:
+    branches: [master]
+  pull_request:
+    branches: [master]
+jobs:
+  Matrix-build:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        os: [windows-latest]
+    steps:
+      - name: checkout
+        uses: actions/checkout@v2
+      - name: checkout submodules
+        run: git submodule update --init --recursive
+      - name: mkdir
+        run: mkdir build && cd build
+      - name: cmake config
+        run: cd build && cmake .. -DTEQP_SNIPPETS=ON -DTEQP_TEQPC=ON -T ClangCL
+      - name: build all Catch tests
+        run: cmake --build build --target catch_tests --config Release
+      - name: build teqpc shared library
+        run: cmake --build build --target teqpc --config Release
+      - name: run all the compiled Catch tests
+        run: cd build && ctest --verbose
\ No newline at end of file
diff --git a/.github/workflows/build_cibuildwheel.yml b/.github/workflows/build_cibuildwheel.yml
index 70ac497172cafc0d71f6074805198e561343e9b9..0eb80423c7f1054a25f80c042e64dca66a1e32f2 100644
--- a/.github/workflows/build_cibuildwheel.yml
+++ b/.github/workflows/build_cibuildwheel.yml
@@ -21,6 +21,10 @@ jobs:
       - name: Install cibuildwheel
         run: python -m pip install cibuildwheel==2.3.1
 
+      - name: Install LLVM (clang-cl) on windows
+        run: choco install -y llvm
+        if: runner.os == 'Windows'
+
       - name: Build wheels
         run: python -m cibuildwheel --output-dir wheelhouse
         # to supply options, put them in 'env', like:
@@ -28,8 +32,9 @@ jobs:
           CIBW_ARCHS: auto64
           CIBW_ARCHS_MACOS: universal2
           CIBW_ENVIRONMENT_MACOS: MACOSX_DEPLOYMENT_TARGET=10.15
+          CIBW_BUILD_VERBOSITY_WINDOWS: 2
           CIBW_SKIP: "*pypy*linux* *cp36*"
 
       - uses: actions/upload-artifact@v2
         with:
-          path: ./wheelhouse/*.whl
+          path: ./wheelhouse/*.whl
\ No newline at end of file
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 83a4c8e65e337c9fadb14a0f6d0b8a8194d51ec4..e22ad41a50965e34f1ead8322a441129e0e17e05 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,3 +1,4 @@
+set(CMAKE_OSX_ARCHITECTURES x86_64;arm64 CACHE INTERNAL "archs for osx")
 cmake_minimum_required(VERSION 3.16)
 project(teqp)
 enable_testing()
@@ -109,6 +110,9 @@ if (NOT TEQP_NO_PYTHON)
     pybind11_add_module(teqp "${pybind11_files}")
     target_include_directories(teqp PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/externals/pybind11_json/include")
     target_link_libraries(teqp PRIVATE autodiff PRIVATE teqpinterface)
+    if (MSVC)
+      target_compile_options(teqp PRIVATE "/Zm1000")
+    endif()
 endif()
 
 file(GLOB catch_tests_files "${CMAKE_CURRENT_SOURCE_DIR}/src/tests/*.cxx")
diff --git a/README.md b/README.md
index e112ec48c8439cfe05620b8099edf40a832a0927..cb9c6751bdf9f5d1eac2dde5b83dbc66f27a488a 100644
--- a/README.md
+++ b/README.md
@@ -18,6 +18,20 @@ Written by Ian Bell, NIST.
 
 [![PyPI version](https://badge.fury.io/py/teqp.svg)](https://badge.fury.io/py/teqp)
 
+* 0.9.2 :
+
+  * Bugfix: ``kmat`` can be set also when specifying ``sigma`` and ``e/kB`` with PC-SAFT
+
+* 0.9.1 :
+
+  * Transcription error in a coefficient of PC-SAFT
+
+* 0.9.0 :
+
+  * Add ability to obtain ancillaries for multifluid model (``see teqp/models/multifluid_ancillaries.hpp``) or the ``build_ancillaries`` method in python
+
+  * Enable ability to use multiprecision with PC-SAFT
+
 * 0.8.1 :
 
   * Replace the ``get_Ar20`` function that was erroneously removed
diff --git a/include/teqp/algorithms/VLE.hpp b/include/teqp/algorithms/VLE.hpp
index cdfa9a976bac8b4ef3129084f0f5a76e789b5113..38b36c105840da90acd96b53e271543cdd3b5c42 100644
--- a/include/teqp/algorithms/VLE.hpp
+++ b/include/teqp/algorithms/VLE.hpp
@@ -57,7 +57,7 @@ public:
         using tdx = TDXDerivatives<Model,TYPE,EigenArray1>;
 
         const TYPE &T = m_T;
-        const TYPE R = m_model.R(molefracs);
+        //const TYPE R = m_model.R(molefracs);
         double R0_over_Rr = R0 / Rr;
         
         auto derL = tdx::template get_Ar0n<2, backend>(m_model, T, rhomolarL, molefracs);
@@ -110,14 +110,14 @@ auto do_pure_VLE_T(Residual &resid, Scalar rhoL, Scalar rhoV, int maxiter) {
         }
         auto v = J.matrix().colPivHouseholderQr().solve(-r0.matrix()).array().eval();
         auto rhovecnew = (rhovec + v).eval();
-        double r00 = static_cast<double>(r0[0]);
-        double r01 = static_cast<double>(r0[1]);
+        //double r00 = static_cast<double>(r0[0]);
+        //double r01 = static_cast<double>(r0[1]);
         
         // If the solution has stopped improving, stop. The change in rhovec is equal to v in infinite precision, but 
         // not when finite precision is involved, use the minimum non-denormal float as the determination of whether
         // the values are done changing
         auto minval = std::numeric_limits<Scalar>::epsilon();
-        double minvaldbl = static_cast<double>(minval);
+        //double minvaldbl = static_cast<double>(minval);
         if (((rhovecnew - rhovec).cwiseAbs() < minval).all()) {
             break;
         }
@@ -302,8 +302,8 @@ Eigen::ArrayXd extrapolate_from_critical(const Model& model, const Scalar& Tc, c
     auto z = (Eigen::ArrayXd(1) << 1.0).finished();
     auto R = model.R(z);
     auto ders = tdx::template get_Ar0n<4>(model, Tc, rhoc, z);
-    auto dpdrho = R*Tc*(1 + 2 * ders[1] + ders[2]); // Should be zero
-    auto d2pdrho2 = R*Tc/rhoc*(2 * ders[1] + 4 * ders[2] + ders[3]); // Should be zero
+    //auto dpdrho = R*Tc*(1 + 2 * ders[1] + ders[2]); // Should be zero
+    //auto d2pdrho2 = R*Tc/rhoc*(2 * ders[1] + 4 * ders[2] + ders[3]); // Should be zero
     auto d3pdrho3 = R*Tc/(rhoc*rhoc)*(6 * ders[2] + 6 * ders[3] + ders[4]);
     auto Ar11 = tdx::template get_Ar11(model, Tc, rhoc, z);
     auto Ar12 = tdx::template get_Ar12(model, Tc, rhoc, z);
@@ -679,7 +679,7 @@ auto trace_VLE_isotherm_binary(const Model &model, Scalar T, VecType rhovecL0, V
             store_point();
         }
 
-        double dtold = dt;
+        //double dtold = dt;
         auto x0_previous = x0;
 
         if (opt.integration_order == 5) {
diff --git a/include/teqp/algorithms/VLLE.hpp b/include/teqp/algorithms/VLLE.hpp
index 01df6601830d551bf8279b059d6cf95979b27f98..b546789cfac6246c6564218bbadef786508b4975 100644
--- a/include/teqp/algorithms/VLLE.hpp
+++ b/include/teqp/algorithms/VLLE.hpp
@@ -174,7 +174,7 @@ namespace teqp {
                 y.push_back(p);
             }
             auto intersections = get_self_intersections(x, y);
-            auto& trace = traces[0];
+            //auto& trace = traces[0];
 
             auto process_intersection = [&](auto& trace, auto& i) {
                 auto rhoL1_j = traces[0][i.j].at("rhoL / mol/m^3").template get<std::valarray<double>>();
diff --git a/include/teqp/models/multifluid_ancillaries.hpp b/include/teqp/models/multifluid_ancillaries.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..5b3eb8f01c5c6db0272a7ffd9246e6421a14afb2
--- /dev/null
+++ b/include/teqp/models/multifluid_ancillaries.hpp
@@ -0,0 +1,49 @@
+#pragma once
+
+#include <valarray>
+#include "nlohmann/json.hpp"
+
+namespace teqp{
+
+	struct VLEAncillary{
+		const double T_r, Tmax, Tmin, reducing_value;
+		const std::string type, description;
+		const std::valarray<double> n, t;
+		const bool using_tau_r, noexp;
+	
+		VLEAncillary(const nlohmann::json &j) :
+			T_r(j.at("T_r")),
+            Tmax(j.at("Tmax")),
+            Tmin(j.at("Tmin")),
+            description(j.at("description")),
+            n(j.at("n").get<std::valarray<double>>()),
+            t(j.at("t").get<std::valarray<double>>()),
+            reducing_value(j.at("reducing_value")),
+            type(j.at("type")),
+            using_tau_r(j.at("using_tau_r")),
+            noexp(type == "rhoLnoexp"){};
+
+		double operator() (double T) const{
+			auto Theta = 1-T/T_r;
+			auto RHS = (pow(Theta, t)*n).sum();
+			if (using_tau_r){
+				RHS *= T_r/T;
+			}
+			if (noexp){
+	            return reducing_value*(1+RHS);
+			}
+	        else{
+	            return exp(RHS)*reducing_value;
+	        }
+		};
+	};
+
+	struct MultiFluidVLEAncillaries {
+		const VLEAncillary rhoL, rhoV, pL, pV;
+		MultiFluidVLEAncillaries(const nlohmann::json& j) :
+			rhoL(VLEAncillary(j.at("rhoL"))),
+			rhoV(VLEAncillary(j.at("rhoV"))),
+			pL((j.contains("pS")) ? VLEAncillary(j.at("pS")) : VLEAncillary(j.at("pL"))),
+			pV((j.contains("pS")) ? VLEAncillary(j.at("pS")) : VLEAncillary(j.at("pV"))){}
+	};
+}
\ No newline at end of file
diff --git a/include/teqp/models/pcsaft.hpp b/include/teqp/models/pcsaft.hpp
index 48b645ffe29b465278629a4d5f8efc30b5cad827..22d5094f1db9eee479dbc69fdaa0fe0a6d87706b 100644
--- a/include/teqp/models/pcsaft.hpp
+++ b/include/teqp/models/pcsaft.hpp
@@ -74,7 +74,7 @@ auto get_a(TYPE mbar) {
 template<typename TYPE>
 auto get_b(TYPE mbar) {
     // See https://stackoverflow.com/a/35170514/1360263
-    static Eigen::ArrayXd b_0 = (Eigen::ArrayXd(7) << 0.724094694, 2.2382791861, -4.0025849485, -21.003576815, 26.855641363, 206.55133841, -355.60235612).finished();
+    static Eigen::ArrayXd b_0 = (Eigen::ArrayXd(7) << 0.7240946941, 2.2382791861, -4.0025849485, -21.003576815, 26.855641363, 206.55133841, -355.60235612).finished();
     static Eigen::ArrayXd b_1 = (Eigen::ArrayXd(7) << -0.5755498075, 0.6995095521, 3.8925673390, -17.215471648, 192.67226447, -161.82646165, -165.20769346).finished();
     static Eigen::ArrayXd b_2 = (Eigen::ArrayXd(7) << 0.0976883116, -0.2557574982, -9.1558561530, 20.642075974, -38.804430052, 93.626774077, -29.666905585).finished();
     return forceeval(b_0.cast<TYPE>() + (mbar - 1.0) / mbar * b_1.cast<TYPE>() + (mbar - 1.0) / mbar * (mbar - 2.0) / mbar * b_2.cast<TYPE>()).eval();
diff --git a/interface/PCSAFT.cpp b/interface/PCSAFT.cpp
index 0d0d8713d3bff816e3b55a71be30a86d4fa526ac..46db39dc708fea5cc327c2dfa5ff36cc0fc74af8 100644
--- a/interface/PCSAFT.cpp
+++ b/interface/PCSAFT.cpp
@@ -18,7 +18,7 @@ void add_PCSAFT(py::module& m) {
 
 	auto wPCSAFT = py::class_<PCSAFTMixture>(m, "PCSAFTEOS")
 	.def(py::init<const std::vector<std::string>&, const Eigen::ArrayXXd&>(), py::arg("names"), py::arg_v("kmat", Eigen::ArrayXXd(0,0), "None"))
-	.def(py::init<const std::vector<SAFTCoeffs>&>(), py::arg("coeffs"))
+	.def(py::init<const std::vector<SAFTCoeffs>&, const Eigen::ArrayXXd&>(), py::arg("coeffs"), py::arg_v("kmat", Eigen::ArrayXXd(0,0), "None"))
 	.def("print_info", &PCSAFTMixture::print_info)
 	.def("max_rhoN", &PCSAFTMixture::max_rhoN<Eigen::ArrayXd>)
 	.def("get_m", &PCSAFTMixture::get_m)
diff --git a/interface/multifluid.cpp b/interface/multifluid.cpp
index 4f09fa49ad56a74514ed0c4f29be0b081c5bd3be..76e1089d5d9926d99c86cff78e3d50538ea001ef 100644
--- a/interface/multifluid.cpp
+++ b/interface/multifluid.cpp
@@ -1,6 +1,7 @@
 #include "pybind11_wrapper.hpp"
 
 #include "teqp/models/multifluid.hpp"
+#include "teqp/models/multifluid_ancillaries.hpp"
 #include "teqp/derivs.hpp"
 #include "teqp/models/mie/lennardjones.hpp"
 
@@ -8,6 +9,24 @@
 
 void add_multifluid(py::module& m) {
 
+    // A single ancillary curve
+    py::class_<VLEAncillary>(m, "VLEAncillary")
+        .def(py::init<const nlohmann::json&>())
+        .def("__call__", &VLEAncillary::operator())
+        .def_readonly("T_r", &VLEAncillary::T_r)
+        .def_readonly("Tmax", &VLEAncillary::Tmax)
+        .def_readonly("Tmin", &VLEAncillary::Tmin)
+        ;
+
+    // The collection of VLE ancillary curves
+    py::class_<MultiFluidVLEAncillaries>(m, "MultiFluidVLEAncillaries")
+        .def(py::init<const nlohmann::json&>())
+        .def_readonly("rhoL", &MultiFluidVLEAncillaries::rhoL)
+        .def_readonly("rhoV", &MultiFluidVLEAncillaries::rhoV)
+        .def_readonly("pL", &MultiFluidVLEAncillaries::pL)
+        .def_readonly("pV", &MultiFluidVLEAncillaries::pV)
+        ;
+
     // Multifluid model
     m.def("build_multifluid_model", &build_multifluid_model, py::arg("components"), py::arg("coolprop_root"), py::arg("BIPcollectionpath") = "", py::arg("flags") = nlohmann::json{}, py::arg("departurepath") = "");
     using MultiFluid = decltype(build_multifluid_model(std::vector<std::string>{"",""},"",""));
diff --git a/interface/multifluid_shared.hpp b/interface/multifluid_shared.hpp
index 2aa39f9be8851819ba3f72ead9a9d3bb7e1b4cb8..6864d23fa704bb4fa2f993cabdadff5857f5743c 100644
--- a/interface/multifluid_shared.hpp
+++ b/interface/multifluid_shared.hpp
@@ -1,5 +1,8 @@
 #pragma once 
 
+#include "teqp/exceptions.hpp"
+#include "teqp/models/multifluid_ancillaries.hpp"
+
 template<typename Model, typename Wrapper>
 void add_multifluid_methods(Wrapper &wMF){
     wMF.def("get_Tcvec", [](const Model& c) { return c.redfunc.Tc; })
@@ -7,5 +10,13 @@ void add_multifluid_methods(Wrapper &wMF){
        .def("get_Tr", [](const Model& c, const Eigen::ArrayXd &molefrac) { return c.redfunc.get_Tr(molefrac); })
        .def("get_rhor", [](const Model& c, const Eigen::ArrayXd &molefrac) { return c.redfunc.get_rhor(molefrac); })
        .def("set_meta", [](Model& c, const std::string &s) { return c.set_meta(s); })
-       .def("get_meta", [](const Model& c) { return c.get_meta(); });
+       .def("get_meta", [](const Model& c) { return c.get_meta(); })
+       .def("build_ancillaries", [](const Model& c) { 
+           if (c.redfunc.Tc.size() != 1) {
+               throw teqp::InvalidArgument("Can only build ancillaries for pure fluids");
+           }
+           auto jancillaries = nlohmann::json::parse(c.get_meta()).at("pures")[0].at("ANCILLARIES");
+           return teqp::MultiFluidVLEAncillaries(jancillaries);
+           })
+       ;
 }
\ No newline at end of file
diff --git a/interface/teqpversion.hpp b/interface/teqpversion.hpp
index 635c7778ec168038b6370dab530e1d9f48debf8d..e831366f97aed4f7a9d75be42d19bdfcd3d5cefe 100644
--- a/interface/teqpversion.hpp
+++ b/interface/teqpversion.hpp
@@ -1,2 +1,2 @@
 #include <string>
-const std::string TEQPVERSION = "0.8.1.dev0";
\ No newline at end of file
+const std::string TEQPVERSION = "0.9.3.dev0";
\ No newline at end of file
diff --git a/setup.py b/setup.py
index 8f980f43444f8e640fa9f726b4f0c01d54e47fbb..1ea918b45251d63d71a866a19910c0aa45be8510 100644
--- a/setup.py
+++ b/setup.py
@@ -55,6 +55,7 @@ class CMakeBuild(build_ext):
 
         if platform.system() == "Windows":
             cmake_args += ['-DCMAKE_LIBRARY_OUTPUT_DIRECTORY_{}={}'.format(cfg.upper(), extdir)]
+            cmake_args += ['-T ClangCL']
             if sys.maxsize > 2**32:
                 cmake_args += ['-A', 'x64']
             build_args += ['--', '/m']
@@ -67,8 +68,17 @@ class CMakeBuild(build_ext):
                                                               self.distribution.get_version())
         if not os.path.exists(self.build_temp):
             os.makedirs(self.build_temp)
-        subprocess.check_call(['cmake', ext.sourcedir] + cmake_args, cwd=self.build_temp, env=env)
-        subprocess.check_call(['cmake', '--build', '.', '--target', 'teqp'] + build_args, cwd=self.build_temp)
+        
+        # Config
+        cmake_elements = ['cmake', ext.sourcedir] + cmake_args
+        print('cmake config command:', ' '.join(cmake_elements))
+        print('running from:', self.build_temp)
+        subprocess.check_call(cmake_elements, cwd=self.build_temp, env=env)
+
+        # Build
+        build_elements = ['cmake', '--build', '.', '--target', 'teqp'] + build_args
+        print('cmake build command:', ' '.join(build_elements))
+        subprocess.check_call(build_elements, cwd=self.build_temp)
 
 init_template  = r'''import os
 
diff --git a/src/boost_pool_test.cpp b/src/boost_pool_test.cpp
index 68b51687e07aff5122d7aaf535665241927faf88..f13164e5e7b8a7a320cd868ecca301c142698bce 100644
--- a/src/boost_pool_test.cpp
+++ b/src/boost_pool_test.cpp
@@ -16,16 +16,33 @@ int simple() {
     // Launch the pool with four threads.
     boost::asio::thread_pool pool(4);
 
-    std::size_t i = 0;
-    for (auto& num : numbers) {
-        auto& o = outputs[i];
-        // Submit a lambda object to the pool.
-        boost::asio::post(pool, [&num, &o]() {o = "as str: " + std::to_string(num); });
-        i++;
-    }
+    auto serial = [&numbers, &outputs]() {
+        std::size_t i = 0;
+        for (auto& num : numbers) {
+            outputs[i] = "as str: " + std::to_string(num);
+            i++;
+        }
+    };
+    
+    auto parallel = [&numbers, &outputs, &pool]() {
+        std::size_t i = 0;
+        for (auto& num : numbers) {
+            auto& o = outputs[i];
+            // Submit a lambda object to the pool.
+            boost::asio::post(pool, [&num, &o]() {o = "as str: " + std::to_string(num); });
+            i++;
+        }
+        // Wait for all tasks in the pool to complete.
+        pool.join();
+    };
 
-    // Wait for all tasks in the pool to complete.
-    pool.join();
+    auto tic = std::chrono::high_resolution_clock::now();
+    serial();
+    auto tic2 = std::chrono::high_resolution_clock::now();
+    parallel();
+    auto tic3 = std::chrono::high_resolution_clock::now();
+    std::cout << std::chrono::duration<double>(tic2 - tic).count() << std::endl;
+    std::cout << std::chrono::duration<double>(tic3 - tic2).count() << std::endl;
 
     for (auto& o : outputs) {
         std::cout << o << std::endl;
diff --git a/src/test_multifluid_ancillaries.cpp b/src/test_multifluid_ancillaries.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7cecb4f0ddce9fa6f9b82f6427bbe978c6a8d6ae
--- /dev/null
+++ b/src/test_multifluid_ancillaries.cpp
@@ -0,0 +1,16 @@
+
+#include "teqp/models/multifluid_ancillaries.hpp"
+#include "teqp/models/multifluid.hpp"
+#include "teqp/algorithms/VLE.hpp"
+
+int main() {
+	auto model = teqp::build_multifluid_model({"n-Propane"}, "../mycp");
+	auto j = nlohmann::json::parse(model.get_meta());
+	auto jancillaries = j.at("pures")[0].at("ANCILLARIES");
+	teqp::MultiFluidVLEAncillaries anc(jancillaries);
+	double T = 340;
+	auto rhoV = anc.rhoV(T), rhoL = anc.rhoL(T);
+	auto rhovec = teqp::pure_VLE_T(model, T, rhoL, rhoV, 10);
+	//auto [rhoLnew, rhoVnew] = rhovec;
+	int rr = 0;
+}
\ No newline at end of file
diff --git a/src/tests/catch_test_PCSAFT.cxx b/src/tests/catch_test_PCSAFT.cxx
index 8ee5144bed38aefc9ff95a253b1c5d1b14af7613..81ba7f20d7dd56ae6f10812a5198644ddce5b834 100644
--- a/src/tests/catch_test_PCSAFT.cxx
+++ b/src/tests/catch_test_PCSAFT.cxx
@@ -7,6 +7,16 @@ using Catch::Approx;
 #include "teqp/models/pcsaft.hpp"
 using namespace teqp::PCSAFT;
 
+TEST_CASE("Single alphar check value", "[PCSAFT]")
+{
+    std::vector<std::string> names = { "Methane" };
+    auto model = PCSAFTMixture(names);
+    double T = 200, Dmolar = 300;
+    auto z = (Eigen::ArrayXd(1) << 1.0).finished();
+    using tdx = teqp::TDXDerivatives<decltype(model), double>;
+    CHECK(tdx::get_Ar00(model, T, Dmolar, z) == Approx(-0.032400020930842724));
+}
+
 TEST_CASE("Check PCSAFT with kij", "[PCSAFT]")
 {
     std::vector<std::string> names = { "Methane", "Ethane" };
@@ -22,4 +32,30 @@ TEST_CASE("Check PCSAFT with kij", "[PCSAFT]")
     SECTION("Incorrectly shaped kij matrix") {
         CHECK_THROWS(PCSAFTMixture(names, kij_bad));
     }
+}
+
+TEST_CASE("Check PCSAFT with kij and coeffs", "[PCSAFT]")
+{
+    std::vector<teqp::PCSAFT::SAFTCoeffs> coeffs;
+    std::vector<double> eoverk = { 120,130 }, m = { 1,2 }, sigma = { 0.9, 1.1 };
+    for (auto i = 0; i < eoverk.size(); ++i) {
+        teqp::PCSAFT::SAFTCoeffs c;
+        c.m = m[i];
+        c.sigma_Angstrom = sigma[i];
+        c.epsilon_over_k = eoverk[i];
+        coeffs.push_back(c);
+    }
+
+    Eigen::ArrayXXd kij_right(2, 2); kij_right.setZero();
+    Eigen::ArrayXXd kij_bad(2, 20); kij_bad.setZero();
+
+    SECTION("No kij") {
+        CHECK_NOTHROW(PCSAFTMixture(coeffs));
+    }
+    SECTION("Correctly shaped kij matrix") {
+        CHECK_NOTHROW(PCSAFTMixture(coeffs, kij_right));
+    }
+    SECTION("Incorrectly shaped kij matrix") {
+        CHECK_THROWS(PCSAFTMixture(coeffs, kij_bad));
+    }
 }
\ No newline at end of file
diff --git a/src/tests/catch_test_multifluid.cxx b/src/tests/catch_test_multifluid.cxx
index 473b09036f923f590a85b3938d93e437ac49e31f..c78ba0205925ef19a8f46d26633e9f5bfdafcdf3 100644
--- a/src/tests/catch_test_multifluid.cxx
+++ b/src/tests/catch_test_multifluid.cxx
@@ -4,6 +4,7 @@
 using Catch::Approx;
 
 #include "teqp/models/multifluid.hpp"
+#include "teqp/models/multifluid_ancillaries.hpp"
 #include "teqp/algorithms/critical_tracing.hpp"
 #include "teqp/algorithms/VLE.hpp"
 #include "teqp/filesystem.hpp"
@@ -139,6 +140,26 @@ TEST_CASE("Check that all pure fluid models can be instantiated", "[multifluid],
     }    
 }
 
+TEST_CASE("Check that all ancillaries can be instantiated and work properly", "[multifluid],[all]") {
+    std::string root = "../mycp";
+    SECTION("With absolute paths to json file") {
+        int counter = 0;
+        for (auto path : get_files_in_folder(root + "/dev/fluids", ".json")) {
+            if (path.filename().stem() == "Methanol") { continue; }
+            CAPTURE(path.string());
+            auto abspath = std::filesystem::absolute(path).string();
+            auto model = build_multifluid_model({ abspath }, root, root + "/dev/mixtures/mixture_binary_pairs.json");
+            auto jancillaries = nlohmann::json::parse(model.get_meta()).at("pures")[0].at("ANCILLARIES");
+            auto anc = teqp::MultiFluidVLEAncillaries(jancillaries);
+            double T = 0.9*anc.rhoL.T_r;
+            auto rhoV = anc.rhoV(T), rhoL = anc.rhoL(T);
+            auto rhovec = teqp::pure_VLE_T(model, T, rhoL, rhoV, 10);
+            counter += 1;
+        }
+        CHECK(counter > 100);
+    }
+}
+
 TEST_CASE("Check that mixtures can also do absolute paths", "[multifluid],[abspath]") {
     std::string root = "../mycp";
     SECTION("With absolute paths to json file") {