diff --git a/src/methods/multicomponent.jl b/src/methods/multicomponent.jl
index 52c673e..e19dea3 100644
--- a/src/methods/multicomponent.jl
+++ b/src/methods/multicomponent.jl
@@ -28,8 +28,6 @@ function IASTModels(m_first::I, m_rest::Vararg{I}) where {I <: IsothermModel}
     return _iastmodels((m_first, m_rest...))
 end
 
-
-
 function loading(model::ExtendedLangmuir{_T, I}, p, T, yᵢ) where {_T, I <: Tuple{Vararg{<:LangmuirS1{_T}}}}
 
     _y = yᵢ/sum(yᵢ) 
@@ -39,30 +37,40 @@ end
 
 function loading(model::ExtendedLangmuir{_T, I}, pᵢ, T) where {_T, I <: Tuple{Vararg{<:LangmuirS1{_T}}}}
 
-    one_plus_Kp = one(eltype(model))
-    loadings = similar(pᵢ, eltype(model))
+    _1_∑kP = one(eltype(model))
+    loadings = similar(pᵢ)
     models = model.isotherms
 
-    for i in eachindex(pᵢ)
+    for i ∈ eachindex(pᵢ)
             loadings[i] = begin
                 M, K₀, E = models[i].M, models[i].K₀, models[i].E
                 MKpy = M*K₀*exp(-E/(Rgas(model)*T))*pᵢ[i]
-                one_plus_Kp += MKpy/M
+                _1_∑kP += MKpy/M
                 MKpy
             end
     end
 
-    return loadings./one_plus_Kp
+    return loadings./_1_∑kP
 end
 
 
 function loading(model::ExtendedLangmuir{_T, I}, p, T, y) where {_T, I <: Tuple{Vararg{<:MultiSite{_T}}}}
         
         unpack_multisite = getfield.(model.isotherms, :isotherms)
-        Ks = collect.(map(x -> getfield.(x, :K₀), unpack_multisite) |> collect)
-        sites_to_match = sortperm.(Ks)
-        models = [ExtendedLangmuir(getindex.(collect(unpack_multisite), getindex.(sites_to_match, i))...) for i in eachindex(Ks)]
+
+        Es = collect.(map(x -> getfield.(x, :E), unpack_multisite) |> collect)
+
+        N_sites =  length(Es[1])
+
+        sites_to_match = sortperm.(Es)
+
+        models = [ExtendedLangmuir(
+            getindex.(collect(unpack_multisite), getindex.(sites_to_match, i))...
+            )
+             for i in 1:N_sites]
+
         loadings = mapreduce(m -> loading(m, p, T, y), +, models)
+
         return loadings
 end
 
@@ -83,4 +91,26 @@ function loading(model::IASTModels, p, T, y; method = IASTNestedLoop(), gas_mode
     return loadings
 end
 
+
+function isosteric_heat(model::ExtendedLangmuir{_T, I}, pᵢ, T; Vᵃ = zeros(eltype(pᵢ)), Vᵍ = Rgas(model).*T./pᵢ) where {_T, I <: Tuple{Vararg{<:LangmuirS1{_T}}}}
+
+    pᵢ_T = [pᵢ; T]
+    f(pᵢ_T) = loading(model, first(pᵢ_T, length(pᵢ)), last(pᵢ_T))
+    ȷ_I = ForwardDiff.jacobian(f, pᵢ_T)
+    
+    ∂pᵢ∂T = @views ȷ_I[:, 1:end - 1]\ȷ_I[:, end]
+    return @. T*(Vᵍ - Vᵃ).*∂pᵢ∂T
+
+end
+
+
+function isosteric_heat(model::ExtendedLangmuir{_T, I}, p, T, yᵢ; Vᵃ = zeros(eltype(p.*yᵢ)), Vᵍ = Rgas(model).*T./(p.*yᵢ)) where {_T, I <: Tuple{Vararg{<:LangmuirS1{_T}}}}
+    _y = yᵢ/sum(yᵢ) 
+    pᵢ = p*_y
+    isosteric_heat(model, pᵢ, T, Vᵃ = Vᵃ, Vᵍ = Vᵍ)
+end
+
+
+
+
 export ExtendedLangmuir, IASTModels
\ No newline at end of file
diff --git a/src/models/thermodynamic_langmuir.jl b/src/models/thermodynamic_langmuir.jl
index 3db4d01..fbb25e9 100644
--- a/src/models/thermodynamic_langmuir.jl
+++ b/src/models/thermodynamic_langmuir.jl
@@ -32,7 +32,7 @@ The activity coefficients `γᵢ` and `γᵩ` are determined using the Gibbs exc
     (E::T, (-Inf, 0.0), "energy parameter")
     (Bᵢᵩ::T, (-Inf, 0.0), "adsorbate-adsorbent interaction coefficient")
 end
-
+ 
 
 function gibbs_excess_free_energy(model::ThermodynamicLangmuir, T, θ)
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 99ad8d4..cf7f520 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -89,6 +89,43 @@ end
 end
 
 
+@testset "Multicomponent Extended Langmuir" begin
+    Lang1 = LangmuirS1(1.727, 16.71e-10, -16152.50)
+    Lang2 = LangmuirS1(4.0, 6.71e-8, -14152.50)
+    multilang = ExtendedLangmuir(Lang1, Lang2)
+    y = [0.3, 0.7]
+    p = 101325.0
+    T = 400.0
+    @test loading(multilang, p, T, y)[1] ≈ Lang1.M * Lang1.K₀*exp(-Lang1.E/(8.31446261815324*T)) * p * y[1] / (1.0 + 
+    Lang2.K₀*exp(-Lang2.E/(8.31446261815324*T))*p*y[2] + Lang1.K₀*exp(-Lang1.E/(8.31446261815324*T))*p*y[1])
+    @test loading(multilang, p, T, y)[2] ≈ Lang2.M * Lang2.K₀*exp(-Lang2.E/(8.31446261815324*T)) * p * y[2] / (1.0 +
+    Lang2.K₀*exp(-Lang2.E/(8.31446261815324*T))*p*y[2] + Lang1.K₀*exp(-Lang1.E/(8.31446261815324*T))*p*y[1])
+end
+
+
+@testset "Multicomponent Extended Multisite Langmuir" begin
+    Lang1 = MultiSite(LangmuirS1(1.727, 16.71e-10, -16152.50), LangmuirS1(4.0, 6.71e-8, -14152.50))
+    Lang2 = MultiSite(LangmuirS1(2.727, 16.71e-12, 1.3*-16152.50), LangmuirS1(7.0, 6.71e-5, 2*-14152.50))
+    multilang = ExtendedLangmuir(Lang1, Lang2)
+    y = [0.6, 0.4]
+    p = 2*101325.0
+    T = 400.0
+
+    den_highE = 1.0 + Lang2.isotherms[2].K₀*exp(-Lang2.isotherms[2].E/(8.31446261815324*T))*p*y[2] + 
+              Lang1.isotherms[1].K₀*exp(-Lang1.isotherms[1].E/(8.31446261815324*T))*p*y[1]
+
+    den_lowE = 1.0 + Lang2.isotherms[1].K₀*exp(-Lang2.isotherms[1].E/(8.31446261815324*T))*p*y[2] +
+                Lang1.isotherms[2].K₀*exp(-Lang1.isotherms[2].E/(8.31446261815324*T))*p*y[1]
+
+    @test loading(multilang, p, T, y)[1] ≈ Lang1.isotherms[1].M * Lang1.isotherms[1].K₀*exp(-Lang1.isotherms[1].E/(8.31446261815324*T)) * p * y[1] / den_highE + 
+    Lang1.isotherms[2].M * Lang1.isotherms[2].K₀*exp(-Lang1.isotherms[2].E/(8.31446261815324*T)) * p * y[1] / den_lowE
+
+    @test loading(multilang, p, T, y)[2] ≈ Lang2.isotherms[1].M * Lang2.isotherms[1].K₀*exp(-Lang2.isotherms[1].E/(8.31446261815324*T)) * p * y[2] / den_lowE + 
+    Lang2.isotherms[2].M * Lang2.isotherms[2].K₀*exp(-Lang2.isotherms[2].E/(8.31446261815324*T)) * p * y[2] / den_highE
+            
+end
+
+
 @testset "IAST" begin
     #10.1002/aic.14684, S.I
     v = @MultiSite{LangmuirS1,LangmuirS1}