diff --git a/Project.toml b/Project.toml
index b8e01d0..b44d7f8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "OrnsteinZernike"
 uuid = "51399321-c693-48c0-994d-05e3df0250f2"
 authors = ["Ilian Pihlajamaa <ilian.pihlajamaa@gmail.com>"]
-version = "0.1.2"
+version = "0.1.3"
 
 [deps]
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
@@ -9,6 +9,7 @@ FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Hankel = "74863788-d124-456e-a676-9b76578dd39e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
diff --git a/src/Thermodynamics.jl b/src/Thermodynamics.jl
index 2ff98a9..29b9175 100644
--- a/src/Thermodynamics.jl
+++ b/src/Thermodynamics.jl
@@ -217,9 +217,14 @@ function compute_compressibility(sol::OZSolution, system::SimpleLiquid{dims, spe
 end
 
 function get_ĉ0(sol::OZSolution, ::SimpleLiquid{dims, species, T1, T2, P}) where {dims, species, T1, T2, P}
-    spl = Spline1D(sol.r, _eachslice(sol.cr, dims=1).*sol.r.^(dims-1))
-    ĉ0 = surface_N_sphere(dims)*integrate(spl, zero(eltype(sol.r)), maximum(sol.r))
+    ĉ0 = zeros(eltype(eltype(sol.ck)), species, species)
+    for i = 1:species
+        for j = 1:species
+            spl = Spline1D(sol.r, sol.cr[:, i, j].*sol.r.^(dims-1))
+            ĉ0[i,j] = surface_N_sphere(dims)*integrate(spl, zero(eltype(sol.r)), maximum(sol.r))
+        end
+    end
     return ĉ0
 end
 
-# function pressure_inconsystency(system, closure)
+