diff --git a/include/teqp/models/pcsaft.hpp b/include/teqp/models/pcsaft.hpp
index 9fd6d3f79f94dd7e96cd637ad5044694e5c8e953..d1d5677adefb0868c5e9e0ed948644523aa33752 100644
--- a/include/teqp/models/pcsaft.hpp
+++ b/include/teqp/models/pcsaft.hpp
@@ -1,6 +1,7 @@
 #pragma once
 
 #include "nlohmann/json.hpp"
+#include "teqp/exceptions.hpp"
 
 namespace teqp {
 namespace PCSAFT {
@@ -177,10 +178,23 @@ private:
         sigma_Angstrom, ///< 
         epsilon_over_k; ///< depth of pair potential divided by Boltzman constant
     std::vector<std::string> names;
-    double k_ij; ///< binary interaction parameter
+    Eigen::ArrayXXd kmat; ///< binary interaction parameter matrix
+
+    void check_kmat(int N) {
+        if (kmat.cols() != kmat.rows()) {
+            throw teqp::InvalidArgument("kmat rows and columns are not identical");
+        }
+        if (kmat.cols() == 0) {
+            kmat.resize(N, N); kmat.setZero();
+        }
+        else if (kmat.cols() != N) {
+            throw teqp::InvalidArgument("kmat needs to be a square matrix the same size as the number of components");
+        }
+    };
 public:
-    PCSAFTMixture(const std::vector<std::string> &names) : names(names)
+    PCSAFTMixture(const std::vector<std::string> &names, const Eigen::ArrayXXd& kmat = {}) : names(names), kmat(kmat)
     {
+        check_kmat(names.size());
         m.resize(names.size());
         mminus1.resize(names.size());
         sigma_Angstrom.resize(names.size());
@@ -195,10 +209,11 @@ public:
             epsilon_over_k[i] = coeff.epsilon_over_k;
             i++;
         }
-        k_ij = 0;
     };
-    PCSAFTMixture(const std::vector<SAFTCoeffs> &coeffs) 
+    PCSAFTMixture(const std::vector<SAFTCoeffs> &coeffs, const Eigen::ArrayXXd &kmat = {}) : kmat(kmat)
     {
+        check_kmat(coeffs.size());
+
         m.resize(coeffs.size());
         mminus1.resize(coeffs.size());
         sigma_Angstrom.resize(coeffs.size());
@@ -213,7 +228,6 @@ public:
             names[i] = coeff.name;
             i++;
         }
-        k_ij = 0;
     };
     auto get_m(){ return m; }
     auto get_sigma_Angstrom() { return sigma_Angstrom; }
@@ -259,7 +273,7 @@ public:
             for (std::size_t j = 0; j < N; ++j) {
                 // Eq. A.5
                 auto sigma_ij = 0.5 * sigma_Angstrom[i] + 0.5 * sigma_Angstrom[j];
-                auto eij_over_k = sqrt(epsilon_over_k[i] * epsilon_over_k[j]) * (1.0 - k_ij);
+                auto eij_over_k = sqrt(epsilon_over_k[i] * epsilon_over_k[j]) * (1.0 - kmat(i,j));
                 c.m2_epsilon_sigma3_bar = c.m2_epsilon_sigma3_bar + mole_fractions[i] * mole_fractions[j] * m[i] * m[j] * eij_over_k / T * pow(sigma_ij, 3);
                 c.m2_epsilon2_sigma3_bar = c.m2_epsilon2_sigma3_bar + mole_fractions[i] * mole_fractions[j] * m[i] * m[j] * pow(eij_over_k / T, 2) * pow(sigma_ij, 3);
             }
diff --git a/interface/PCSAFT.cpp b/interface/PCSAFT.cpp
index 578b830dc11bf6c0e25574abb50c58374d528cc9..5796f5788cd68edb3c4d4ccae4835d80dc073c5b 100644
--- a/interface/PCSAFT.cpp
+++ b/interface/PCSAFT.cpp
@@ -17,7 +17,7 @@ void add_PCSAFT(py::module& m) {
 	;
 
 	auto wPCSAFT = py::class_<PCSAFTMixture>(m, "PCSAFTEOS")
-	.def(py::init<const std::vector<std::string>&>(), py::arg("names"))
+	.def(py::init<const std::vector<std::string>&, const Eigen::ArrayXXd&>(), py::arg("names"), py::arg_v("kmat", Eigen::ArrayXXd(2,2).setZero(), "None"))
 	.def(py::init<const std::vector<SAFTCoeffs>&>(), py::arg("coeffs"))
 	.def("print_info", &PCSAFTMixture::print_info)
 	.def("max_rhoN", &PCSAFTMixture::max_rhoN<Eigen::ArrayXd>)
diff --git a/src/tests/catch_test_PCSAFT.cxx b/src/tests/catch_test_PCSAFT.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..8ee5144bed38aefc9ff95a253b1c5d1b14af7613
--- /dev/null
+++ b/src/tests/catch_test_PCSAFT.cxx
@@ -0,0 +1,25 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/catch_approx.hpp>
+
+using Catch::Approx;
+
+#include "teqp/core.hpp"
+#include "teqp/models/pcsaft.hpp"
+using namespace teqp::PCSAFT;
+
+TEST_CASE("Check PCSAFT with kij", "[PCSAFT]")
+{
+    std::vector<std::string> names = { "Methane", "Ethane" };
+    Eigen::ArrayXXd kij_right(2, 2); kij_right.setZero();
+    Eigen::ArrayXXd kij_bad(2, 20); kij_bad.setZero();
+
+    SECTION("No kij") {
+        CHECK_NOTHROW(PCSAFTMixture(names));
+    }
+    SECTION("Correctly shaped kij matrix") {
+        CHECK_NOTHROW(PCSAFTMixture(names, kij_right));
+    }
+    SECTION("Incorrectly shaped kij matrix") {
+        CHECK_THROWS(PCSAFTMixture(names, kij_bad));
+    }
+}
\ No newline at end of file