From 0bbdb956fb4e86a6afb7d6645407035c559ce022 Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Thu, 5 Jun 2025 12:11:17 -0700
Subject: [PATCH 01/28] add objects and evaluation functions for thermochemical
 coverage dependence

---
 src/Calculators/ThermoCoverageDependence.jl | 34 +++++++++++++++++++++
 1 file changed, 34 insertions(+)
 create mode 100644 src/Calculators/ThermoCoverageDependence.jl

diff --git a/src/Calculators/ThermoCoverageDependence.jl b/src/Calculators/ThermoCoverageDependence.jl
new file mode 100644
index 000000000..29aa25645
--- /dev/null
+++ b/src/Calculators/ThermoCoverageDependence.jl
@@ -0,0 +1,34 @@
+using Parameters
+
+abstract type AbstractThermoCoverageDependence end
+export AbstractThermoCoverageDependence
+
+struct EmptyThermoCoverageDependence <: AbstractThermoCoverageDependence end
+@inline getcovdepenthalpycorrection(covdep::EmptyThermoCoverageDependence,T,coverages) = zero(T)
+@inline getcovdepentropycorrection(covdep::EmptyThermoCoverageDependence,T,coverages) = zero(T)
+@inline getcovdepheatcapacitycorrection(covdep::EmptyThermoCoverageDependence,T,coverages) = zero(T)
+export EmptyThermoCoverageDependence
+
+"""
+Polynomial energy coverage dependence
+correction is applied to enthalpy
+poly is a dictionary mapping species names to polynomial coverage dependence corrections to their energy/enthalpy
+polynomial ordering is a_1 + a_2 * theta + a_3 * theta^2 etc.
+"""
+@with_kw struct PolynomialThermoEnergyCoverageDependence <: AbstractThermoCoverageDependence
+    polys::Dict{String,Vector{Float64}} = Dict{String,Vector{Float64}}()
+    indpolys::Dict{Int64,Vector{Float64}} = Dict{Int64,Vector{Float64}}()
+end
+@inline function PolynomialThermoEnergyCoverageDependence(polys)
+    return PolynomialThermoEnergyCoverageDependence(polys,Dict{Int64,Vector{Float64}}())
+end
+@inline function getcovdepenthalpycorrection(covdep::PolynomialThermoEnergyCoverageDependence,T,coverages)
+    E = zero(coverages[1])
+    for (ind,poly) in covdep.indpolys
+        E += evalpoly(coverages[ind],poly)
+    end
+    return E
+end
+@inline getcovdepentropycorrection(covdep::PolynomialThermoEnergyCoverageDependence,T,coverages) = zero(coverages[1])
+@inline getcovdepheatcapacitycorrection(covdep::PolynomialThermoEnergyCoverageDependence,T,coverages) = zero(coverages[1])
+export PolynomialThermoEnergyCoverageDependence
\ No newline at end of file

From 8c00e9c3be7defab358b0e71db445f2a60b730be Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Thu, 5 Jun 2025 12:11:40 -0700
Subject: [PATCH 02/28] add objects and evaluation functions for rate
 coefficient coverage dependence

---
 src/Calculators/RateCoverageDependence.jl | 55 +++++++++++++++++++++++
 1 file changed, 55 insertions(+)
 create mode 100644 src/Calculators/RateCoverageDependence.jl

diff --git a/src/Calculators/RateCoverageDependence.jl b/src/Calculators/RateCoverageDependence.jl
new file mode 100644
index 000000000..942c6ddaa
--- /dev/null
+++ b/src/Calculators/RateCoverageDependence.jl
@@ -0,0 +1,55 @@
+using Parameters
+
+abstract type AbstractRateCoverageDependence end
+export AbstractRateCoverageDependence
+
+struct EmptyRateCoverageDependence <: AbstractRateCoverageDependence end
+export EmptyRateCoverageDependence
+
+getcovdepactivationbarriercorrection(covdep::EmptyRateCoverageDependence,T,coverages) = 0.0
+getcovdepfactorcorrection(covdep::EmptyRateCoverageDependence,T,coverages) = 1.0
+
+export getcovdepactivationbarriercorrection
+export getcovdepfactorcorrection
+
+"""
+Polynomial rate coefficient coverage dependence
+correction form is: 10^(sum(a_i*theta_k^i,i)) * theta_k^m * exp(sum(E_i*theta_k^i,i)/(RT))
+dictionaries map the name of the species (k) to polynomials with ordering a_1 + a_2 * theta + a_3 * theta^2 etc.
+Epolys corresponds to E_i (array of polynomial coefficients), avals corresponds to a_i (array of polynomial coefficients), 
+ms corresponds to m (float)
+"""
+@with_kw struct PolynomialRateCoverageDependence{N,K,L} <: AbstractRateCoverageDependence
+    avals::Dict{String,L} = Dict()
+    Epolys::Dict{String,K} = Dict()
+    ms::Dict{String,L} = Dict()
+    indavals::Dict{N,L} = Dict()
+    indEpolys::Dict{N,K} = Dict()
+    indms::Dict{N,L} = Dict()
+end
+@inline function getcovdepactivationbarriercorrection(covdep::PolynomialRateCoverageDependence,T,coverages)
+    E = 0.0
+    for (ind,Epoly) in covdep.indEpolys
+        E += evalpoly(coverages[ind],indEpoly)
+    end
+    return E
+end
+@inline function getcovdepfactorcorrection(covdep::PolynomialRateCoverageDependence,T,coverages)
+    av = 0.0
+    for (ind,a) in covdep.indavals
+       av += coverages[ind]*a
+    end
+
+    if av != 0
+        A = 10.0^av
+    else 
+        A = 1.0
+    end
+
+    for (ind,m) in covdep.indms 
+        A *= coverages[ind]^m
+    end
+    
+    return A
+end
+export PolynomialRateCoverageDependence
\ No newline at end of file

From 9b39382e64e0933cc2cdb69f32f3126a79f21c29 Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Thu, 5 Jun 2025 12:12:30 -0700
Subject: [PATCH 03/28] adapt individual thermochemical evaluation to handle
 coverage dependence

---
 src/Calculators/Thermo.jl | 56 +++++++++++++++++++++++----------------
 1 file changed, 33 insertions(+), 23 deletions(-)

diff --git a/src/Calculators/Thermo.jl b/src/Calculators/Thermo.jl
index 534873a56..936c69446 100644
--- a/src/Calculators/Thermo.jl
+++ b/src/Calculators/Thermo.jl
@@ -3,18 +3,19 @@ using Parameters
 abstract type AbstractThermo end
 export AbstractThermo
 
-@inline function getGibbs(th::P,T::N) where {N<:Number,P<:AbstractThermo}
-    return @fastmath getEnthalpy(th,T)-T*getEntropy(th,T)
+@inline function getGibbs(th::P,T::N;coverages=nothing) where {N<:Number,P}
+    return @fastmath (getEnthalpy(th,T)+getcovdepenthalpycorrection(th.covdep,T,coverages))-T*(getEntropy(th,T)+getcovdepentropycorrection(th.covdep,T,coverages))
 end
 
-@with_kw struct NASApolynomial <: AbstractThermo
+@with_kw struct NASApolynomial{H<:AbstractThermoCoverageDependence} <: AbstractThermo
     coefs::Array{Float64,1}
     Tmin::Float64
     Tmax::Float64
+    covdep::H = EmptyThermoCoverageDependence()
 end
 export NASApolynomial
 
-@inline function calcHSCpdless(poly::NASApolynomial,T::X) where {X<:Real}
+@inline function calcHSCpdless(poly::NASApolynomial,T::X;coverages=nothing) where {X<:Real}
     if length(poly.coefs) != 7
         Tpoly0 = T
         Tpoly1 = T*T
@@ -58,43 +59,51 @@ export NASApolynomial
         hdivRT =ct2+0.5*ct3+0.33333333333*ct4+0.25*ct5+0.2*ct6+poly.coefs[6]*Tpoly4
         sdivR = Tpoly6*ct2 + ct3 + 0.5*ct4 + 0.33333333333*ct5 + 0.25*ct6 + poly.coefs[7]
     end
+    hdivRT += getcovdepenthalpycorrection(th.covdep,T,coverages)/(R*T)
+    sdivR += getcovdepentropycorrection(th.covdep,T,coverages)/R
+    cpdivR += getcovdepheatcapacitycorrection(th.covdep,T,coverages)/R
     return (cpdivR,hdivRT,sdivR)::Tuple{Float64,Float64,Float64}
 end
 export calcHSCpdless
 
-@inline function getHeatCapacity(poly::NASApolynomial,T::N) where {N<:Number}
+@inline function getHeatCapacity(poly::NASApolynomial,T::N,coverages=nothing) where {N<:Number}
     if length(poly.coefs) == 9
-        return @views @inbounds @fastmath evalpoly(T,poly.coefs[1:7])/T^2*R
+        return @views @inbounds @fastmath evalpoly(T,poly.coefs[1:7])/T^2*R + getcovdepheatcapacitycorrection(poly.covdep,T,coverages)
     elseif length(poly.coefs) == 7
-        return @views @inbounds @fastmath evalpoly(T,poly.coefs[1:5])*R
+        return @views @inbounds @fastmath evalpoly(T,poly.coefs[1:5])*R + getcovdepheatcapacitycorrection(poly.covdep,T,coverages)
     else
         throw(error("NASA polynomial has a number of coefficients not equal to 9 or 7"))
     end
 end
 
-@inline function getEntropy(poly::NASApolynomial,T::N) where {N<:Number}
+@inline function getEntropy(poly::NASApolynomial,T::N,coverages=nothing) where {N<:Number}
     if length(poly.coefs) == 9
-        return @views @inbounds @fastmath ((-poly.coefs[1]/(2*T)-poly.coefs[2])/T+poly.coefs[3]*log(T)+T*evalpoly(T,poly.coefs[4:end-2]./(1:4))+poly.coefs[end])*R
+        return @views @inbounds @fastmath ((-poly.coefs[1]/(2*T)-poly.coefs[2])/T+poly.coefs[3]*log(T)+T*evalpoly(T,poly.coefs[4:end-2]./(1:4))+poly.coefs[end])*R + getcovdepentropycorrection(poly.covdep,T,coverages)
     elseif length(poly.coefs) == 7
-        return @views @inbounds @fastmath (poly.coefs[1]*log(T)+T*evalpoly(T,poly.coefs[2:end-2]./(1:4))+poly.coefs[end])*R
+        return @views @inbounds @fastmath (poly.coefs[1]*log(T)+T*evalpoly(T,poly.coefs[2:end-2]./(1:4))+poly.coefs[end])*R + getcovdepentropycorrection(poly.covdep,T,coverages)
     else
         throw(error("NASA polynomial has a number of coefficients not equal to 9 or 7"))
     end
 end
 
-@inline function getEnthalpy(poly::NASApolynomial,T::N) where {N<:Number}
+@inline function getEnthalpy(poly::NASApolynomial,T::N,coverages=nothing) where {N<:Number}
     if length(poly.coefs) == 9
-        return @views @inbounds @fastmath ((-poly.coefs[1]/T+poly.coefs[2]*log(T))/T+evalpoly(T,poly.coefs[3:end-2]./(1:5)))*R*T+poly.coefs[end-1]*R
+        return @views @inbounds @fastmath ((-poly.coefs[1]/T+poly.coefs[2]*log(T))/T+evalpoly(T,poly.coefs[3:end-2]./(1:5)))*R*T+poly.coefs[end-1]*R+th.covdep(T,coverages) + getcovdepenthalpycorrection(poly.covdep,T,coverages)
     elseif length(poly.coefs) == 7
-        return @views @inbounds @fastmath evalpoly(T,poly.coefs[1:end-2]./(1:5))*R*T+poly.coefs[end-1]*R
+        return @views @inbounds @fastmath evalpoly(T,poly.coefs[1:end-2]./(1:5))*R*T+poly.coefs[end-1]*R + getcovdepenthalpycorrection(poly.covdep,T,coverages)
     else
         throw(error("NASA polynomial has a number of coefficients not equal to 9 or 7"))
     end
 end
 
-@with_kw struct NASA{T<:AbstractThermoUncertainty} <: AbstractThermo
+@with_kw struct NASA{T<:AbstractThermoUncertainty,H<:AbstractThermoCoverageDependence} <: AbstractThermo
     polys::Array{NASApolynomial,1}
     unc::T = EmptyThermoUncertainty()
+    covdep::H = EmptyThermoCoverageDependence()
+end
+function NASA(polys,unc=EmptyThermoUncertainty(),covdep=EmptyThermoCoverageDependence())
+    ps = [polys.covdep isa EmptyThermoCoverageDependence ? NASApolynomial(p.coefs,p.Tmin,p.Tmax,covdep) : p for p in polys]
+    return NASA(ps,unc,covdep)
 end
 export NASA
 
@@ -111,13 +120,13 @@ export NASA
 end
 export selectPoly
 
-@inline getHeatCapacity(nasa::NASA,T::N) where {N<:Number} = getHeatCapacity(selectPoly(nasa,T),T)
+@inline getHeatCapacity(nasa::NASA,T::N) where {N<:Number} = getHeatCapacity(selectPoly(nasa,T),T) #we do not add covdep corrections in NASA because it is also added in the NASAPolynomial, so it would double count
 @inline getEntropy(nasa::NASA,T::N) where {N<:Number} = getEntropy(selectPoly(nasa,T),T)
 @inline getEnthalpy(nasa::NASA,T::N) where {N<:Number} = getEnthalpy(selectPoly(nasa,T),T)
 @inline getGibbs(nasa::NASA,T::N) where {N<:Number} = getGibbs(selectPoly(nasa,T),T)
 @inline calcHSCpdless(nasa::NASA,T::N) where {N<:Real}= calcHSCpdless(selectPoly(nasa,T),T)
 
-@with_kw struct Wilhoit{N,Q,T,P,U,R<:Number,M<:AbstractThermoUncertainty} <: AbstractThermo
+@with_kw struct Wilhoit{N,Q,T,P,U,R<:Number,M<:AbstractThermoUncertainty,H<:AbstractThermoCoverageDependence} <: AbstractThermo
     Cp0::N
     Cpinf::T
     coefs::Array{Q,1}
@@ -125,15 +134,16 @@ export selectPoly
     S0::U
     B::R
     unc::M = EmptyThermoUncertainty()
+    covdep::H = EmptyThermoCoverageDependence()
 end
 export Wilhoit
 
-@inline function getHeatCapacity(w::Wilhoit,T::N) where {N<:Number}
+@inline function getHeatCapacity(w::Wilhoit,T::N,coverages=nothing) where {N<:Number}
     @fastmath y = T/(T+w.B)
-    return @fastmath w.Cp0 + (w.Cpinf-w.Cp0)*y^2*(1+(y-1)*evalpoly(y,w.coefs))
+    return @fastmath w.Cp0 + (w.Cpinf-w.Cp0)*y^2*(1+(y-1)*evalpoly(y,w.coefs)) + getcovdepheatcapacitycorrection(th.covdep,T,coverages)
 end
 
-@inline function getEnthalpy(w::Wilhoit,T::N) where {N<:Number}
+@inline function getEnthalpy(w::Wilhoit,T::N,coverages=nothing) where {N<:Number}
     @fastmath y = T/(T+w.B)
     return @views @fastmath @inbounds w.H0 + w.Cp0 * T - (w.Cpinf - w.Cp0) * T * (
             y * y * ((3 * w.coefs[1] + sum(w.coefs[2:end])) / 6. +
@@ -141,12 +151,12 @@ end
                      (5 * w.coefs[3] + w.coefs[4]) * y^2 / 20. +
                      w.coefs[4] * y^3 / 5.) +
             (2 + sum(w.coefs)) * (y / 2. - 1 + (1.0 / y - 1.) * log(w.B + T))
-        )
+        ) + getcovdepenthalpycorrection(th.covdep,T,coverages)
 end
 
-@inline function getEntropy(w::Wilhoit,T::N) where {N<:Number}
+@inline function getEntropy(w::Wilhoit,T::N,coverages=nothing) where {N<:Number}
     @fastmath y = T/(T+w.B)
-    return @fastmath w.S0 + w.Cpinf*log(T)-(w.Cpinf-w.Cp0)*(log(y)+y*(1+y*evalpoly(y,w.coefs./(2:5))))
+    return @fastmath w.S0 + w.Cpinf*log(T)-(w.Cpinf-w.Cp0)*(log(y)+y*(1+y*evalpoly(y,w.coefs./(2:5)))) + getcovdepentropycorrection(th.covdep,T,coverages)
 end
 
 @with_kw struct ConstantG{B<:Number,J} <: AbstractThermo
@@ -155,7 +165,7 @@ end
 end
 export ConstantG
 
-@inline function getGibbs(cg::ConstantG,T::B) where {B<:Number}
+@inline function getGibbs(cg::ConstantG,T::B,coverages=nothing) where {B<:Number}
     return cg.G
 end
 

From 4639624f5247dad253af24e6bfc72c68a6dcc7c1 Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Thu, 5 Jun 2025 12:12:53 -0700
Subject: [PATCH 04/28] adapt rate coefficient evaluation to handle coverage
 dependence

---
 src/Calculators/Rate.jl | 47 ++++++++++++++++++++++++-----------------
 1 file changed, 28 insertions(+), 19 deletions(-)

diff --git a/src/Calculators/Rate.jl b/src/Calculators/Rate.jl
index c873c58a0..65d165a30 100644
--- a/src/Calculators/Rate.jl
+++ b/src/Calculators/Rate.jl
@@ -7,35 +7,38 @@ export AbstractRate
 abstract type AbstractFalloffRate <: AbstractRate end
 export AbstractFalloffRate
 
-@with_kw struct Arrhenius{N<:Real,K<:Real,Q<:Real,P<:AbstractRateUncertainty} <: AbstractRate
+@with_kw struct Arrhenius{N<:Real,K<:Real,Q<:Real,P<:AbstractRateUncertainty,H<:AbstractRateCoverageDependence} <: AbstractRate
         A::N
         n::K
         Ea::Q
         unc::P = EmptyRateUncertainty()
+        covdep::H = EmptyRateCoverageDependence()
 end
-@inline (arr::Arrhenius)(;T::Q,P::N=0.0,C::S=0.0,phi=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp(-arr.Ea/(R*T))
-@inline (arr::Arrhenius)(T::Q;P::N=0.0,C::S=0.0,phi=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp(-arr.Ea/(R*T))
+@inline (arr::Arrhenius)(;T::Q,P::N=0.0,C::S=0.0,phi=0.0,coverages=nothing) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp(-(arr.Ea+getcovdepactivationbarriercorrection(arr.covdep,T,coverages))/(R*T))*getcovdepfactorcorrection(arr.covdep,T,coverages)
+@inline (arr::Arrhenius)(T::Q;P::N=0.0,C::S=0.0,phi=0.0,coverages=nothing) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp(-(arr.Ea+getcovdepactivationbarriercorrection(arr.covdep,T,coverages))/(R*T))*getcovdepfactorcorrection(arr.covdep,T,coverages)
 export Arrhenius
 
-@with_kw struct StickingCoefficient{N<:Real,K<:Real,Q<:Real,P<:AbstractRateUncertainty} <: AbstractRate
+@with_kw struct StickingCoefficient{N<:Real,K<:Real,Q<:Real,P<:AbstractRateUncertainty,H<:AbstractRateCoverageDependence} <: AbstractRate
         A::N
         n::K
         Ea::Q
         unc::P = EmptyRateUncertainty()
+        covdep::H = EmptyRateCoverageDependence()
 end
-@inline (arr::StickingCoefficient)(;T::Q,P::N=0.0,C::S=0.0,phi=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath min(arr.A*T^arr.n*exp(-arr.Ea/(R*T)),1.0)
-@inline (arr::StickingCoefficient)(T::Q;P::N=0.0,C::S=0.0,phi=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath min(arr.A*T^arr.n*exp(-arr.Ea/(R*T)),1.0)
+@inline (arr::StickingCoefficient)(;T::Q,P::N=0.0,C::S=0.0,phi=0.0,coverages=nothing) where {Q<:Real,N<:Real,S<:Real} = @fastmath min(arr.A*T^arr.n*exp(-(arr.Ea+getcovdepactivationbarriercorrection(arr.covdep,T,coverages))/(R*T))*getcovdepfactorcorrection(arr.covdep,T,coverages),1.0)
+@inline (arr::StickingCoefficient)(T::Q;P::N=0.0,C::S=0.0,phi=0.0,coverages=nothing) where {Q<:Real,N<:Real,S<:Real} = @fastmath min(arr.A*T^arr.n*exp(-(arr.Ea+getcovdepactivationbarriercorrection(arr.covdep,T,coverages))/(R*T))*getcovdepfactorcorrection(arr.covdep,T,coverages),1.0)
 export StickingCoefficient
 
-@with_kw struct Arrheniusq{N<:Real,K<:Real,Q<:Real,P<:AbstractRateUncertainty,B} <: AbstractRate
+@with_kw struct Arrheniusq{N<:Real,K<:Real,Q<:Real,P<:AbstractRateUncertainty,B,H<:AbstractRateCoverageDependence} <: AbstractRate
         A::N
         n::K
         Ea::Q
         q::B = 0.0
         unc::P = EmptyRateUncertainty()
+        covdep::H = EmptyRateCoverageDependence()
 end
-@inline (arr::Arrheniusq)(;T::Q,P::N=0.0,C::S=0.0,phi=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp((-arr.Ea-arr.q*F*phi)/(R*T))
-@inline (arr::Arrheniusq)(T::Q;P::N=0.0,C::S=0.0,phi=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp((-arr.Ea-arr.q*F*phi)/(R*T))
+@inline (arr::Arrheniusq)(;T::Q,P::N=0.0,C::S=0.0,phi=0.0,coverages=nothing) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp((-(arr.Ea+getcovdepactivationbarriercorrection(arr.covdep,T,coverages))-arr.q*F*phi)/(R*T))*getcovdepfactorcorrection(arr.covdep,T,coverages)
+@inline (arr::Arrheniusq)(T::Q;P::N=0.0,C::S=0.0,phi=0.0,coverages=nothing) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp((-(arr.Ea+getcovdepactivationbarriercorrection(arr.covdep,T,coverages))-arr.q*F*phi)/(R*T))*getcovdepfactorcorrection(arr.covdep,T,coverages)
 export Arrheniusq
 
 @with_kw struct PdepArrhenius{T<:Real,Q<:AbstractRateUncertainty,Z<:AbstractRate} <: AbstractRate
@@ -45,7 +48,7 @@ export Arrheniusq
 end
 PdepArrhenius(Ps::Array{Q,1},arrs::Array{Z,1}) where {Q<:Real,Z<:AbstractRate} = PdepArrhenius(sort(Ps),arrs)
 
-@inline function (parr::PdepArrhenius)(;T::Q=nothing,P::V=nothing,C::S=0.0,phi=0.0) where {Q<:Real,V<:Real,S<:Real}
+@inline function (parr::PdepArrhenius)(;T::Q=nothing,P::V=nothing,C::S=0.0,phi=0.0,coverages=nothing) where {Q<:Real,V<:Real,S<:Real}
     inds = getBoundingIndsSorted(P,parr.Ps)::Tuple{Int64,Int64}
 
     if inds[2] == -1
@@ -60,15 +63,16 @@ PdepArrhenius(Ps::Array{Q,1},arrs::Array{Z,1}) where {Q<:Real,Z<:AbstractRate} =
 end
 export PdepArrhenius
 
-@with_kw struct MultiArrhenius{Q<:AbstractRateUncertainty} <: AbstractRate
+@with_kw struct MultiArrhenius{Q<:AbstractRateUncertainty,H<:AbstractRateCoverageDependence} <: AbstractRate
     arrs::Array{Arrhenius,1}
     unc::Q = EmptyRateUncertainty()
+    covdep::H = EmptyRateCoverageDependence()
 end
 
-@inline function (marr::MultiArrhenius)(;T::Q,P::R=0.0,C::S=0.0,phi=0.0) where {Q<:Real,R<:Real,S<:Real}
+@inline function (marr::MultiArrhenius)(;T::Q,P::R=0.0,C::S=0.0,phi=0.0,coverages=nothing) where {Q<:Real,R<:Real,S<:Real}
     out = 0.0
     for arr in marr.arrs
-        @fastmath out += arr(T)
+        @fastmath out += arr(T;coverages=coverages)
     end
     return out
 end
@@ -79,7 +83,7 @@ export MultiArrhenius
     unc::Q = EmptyRateUncertainty()
 end
 
-@inline function (parr::MultiPdepArrhenius)(;T::Q,P::R=0.0,C::S=0.0,phi=0.0) where {Q<:Real,R<:Real,S<:Real}
+@inline function (parr::MultiPdepArrhenius)(;T::Q,P::R=0.0,C::S=0.0,phi=0.0,coverages=nothing) where {Q<:Real,R<:Real,S<:Real}
     out = 0.0
     for pdar in parr.parrs
         @fastmath out += pdar(T=T,P=P)
@@ -95,7 +99,7 @@ export MultiPdepArrhenius
     unc::Q = EmptyRateUncertainty()
 end
 
-(tbarr::ThirdBody)(;T::Q=nothing,P::R=0.0,C::S=nothing,phi=0.0) where {Q<:Real,R<:Real,S<:Real} = C*(tbarr.arr(T))
+(tbarr::ThirdBody)(;T::Q=nothing,P::R=0.0,C::S=nothing,phi=0.0,coverages=nothing) where {Q<:Real,R<:Real,S<:Real} = C*(tbarr.arr(T))
 export ThirdBody
 
 @with_kw struct Lindemann{N<:Integer,K<:AbstractFloat,Q<:AbstractRateUncertainty} <: AbstractFalloffRate
@@ -106,7 +110,7 @@ export ThirdBody
     unc::Q = EmptyRateUncertainty()
 end
 
-@inline function (lnd::Lindemann)(;T::Q=nothing,P::R=0.0,C::S=nothing,phi=0.0) where {Q<:Real,R<:Real,S<:Real}
+@inline function (lnd::Lindemann)(;T::Q=nothing,P::R=0.0,C::S=nothing,phi=0.0,coverages=nothing) where {Q<:Real,R<:Real,S<:Real}
     k0 = lnd.arrlow(T=T)
     kinf = lnd.arrhigh(T=T)
     @fastmath Pr = k0*C/kinf
@@ -126,7 +130,7 @@ export Lindemann
     unc::R = EmptyRateUncertainty()
 end
 
-@inline function (tr::Troe)(;T::Q,P::R=0.0,C::S=nothing,phi=0.0) where {Q<:Real,R<:Real,S<:Real}
+@inline function (tr::Troe)(;T::Q,P::R=0.0,C::S=nothing,phi=0.0,coverages=nothing) where {Q<:Real,R<:Real,S<:Real}
     k0 = tr.arrlow(T=T)
     kinf = tr.arrhigh(T=T)
     @fastmath Pr = k0*C/kinf
@@ -195,7 +199,7 @@ export getredtemp
 end
 export getredpress
 
-@inline function (ch::Chebyshev)(;T::N,P::Q=0.0,C::B=0.0,phi=0.0) where {N<:Real,B<:Real,Q<:Real}
+@inline function (ch::Chebyshev)(;T::N,P::Q=0.0,C::B=0.0,phi=0.0,coverages=nothing) where {N<:Real,B<:Real,Q<:Real}
     k = 0.0
     Tred = getredtemp(ch,T)
     Pred = getredpress(ch,P)
@@ -210,7 +214,12 @@ end
 export Chebyshev
 
 @inline function getkineticstype(kin::B) where {B<:AbstractRate}
-    return extracttypename(typeof(kin).name)
+    cdep = [string(x) for x in typeof(kin).parameters if occursin("CoverageDependence",string(x))]
+    if length(cdep) == 0 || cdep[1] == "EmptyRateCoverageDependence"
+        return extracttypename(typeof(kin).name)
+    else
+        return extracttypename(typeof(kin).name) * "_" * cdep[1]
+    end
 end
 
 @inline function getkineticstype(kin::PdepArrhenius)

From 4171167095b7dfdc3c424a544db8f4c08713acdb Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Thu, 5 Jun 2025 12:16:41 -0700
Subject: [PATCH 05/28] adapt vectorized thermochemistry object construction to
 handle coverage dependence

in particular this makes it so you can stick the coverage dependence on NASA and it will get passed down to
its NASAPolynomial sub objects
---
 src/Calculators/Thermovec.jl | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

diff --git a/src/Calculators/Thermovec.jl b/src/Calculators/Thermovec.jl
index 1c99db8bf..2d66b7c98 100644
--- a/src/Calculators/Thermovec.jl
+++ b/src/Calculators/Thermovec.jl
@@ -5,19 +5,22 @@ export AbstractThermovec
 
 include("Thermo.jl")
 
-@with_kw struct NASApolynomialvec <: AbstractThermo
+@with_kw struct NASApolynomialvec{Q} <: AbstractThermo
     coefs::Array{Float64,2}
     Tmin::Float64
     Tmax::Float64
+    covdeps::Q = nothing
 end
 
-@with_kw struct NASAvec{T<:AbstractThermoUncertainty} <: AbstractThermo
+@with_kw struct NASAvec{T<:AbstractThermoUncertainty,Q} <: AbstractThermo
     polys::Array{NASApolynomialvec,1}
     unc::T = EmptyThermoUncertainty()
+    covdeps::Q = nothing
 end
 function NASAvec(nasas::B) where {B<:Array}
     Tmin = nasas[1].polys[1].Tmin
     Tmax = nasas[1].polys[end].Tmax
+    covdeps = [nasa.covdep for nasa in nasas]
     Ts = Array{Float64,1}()
     for nasa in nasas
         for (i,poly) in enumerate(nasa.polys)
@@ -53,12 +56,10 @@ function NASAvec(nasas::B) where {B<:Array}
                 polyvec[:,j] = nasapoly.coefs
             end
         end
-        nasapvec = NASApolynomialvec(polyvec,Ts[i],Ts[i+1])
+        nasapvec = NASApolynomialvec(polyvec,Ts[i],Ts[i+1],nothing)
         push!(polyvecs,nasapvec)
     end
-    
-    return NASAvec(polys=polyvecs) 
-        
+    return NASAvec(polys=polyvecs,covdeps=covdeps)
 end
 @inline function selectPoly(nasa::NASAvec,T::N) where {N<:Real}
     """

From 0aa2ad747c60988c33f54177d71e61d8ed3ec59c Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Thu, 5 Jun 2025 12:18:22 -0700
Subject: [PATCH 06/28] handle coverage dependence in vectorized
 thermochemistry evaluation

in particular this adjust the dispatching to help handle the fact that coverage dependence complicates automatic differentiation
---
 src/Calculators/Thermovec.jl | 60 ++++++++++++++++++++++++++++++++++--
 1 file changed, 57 insertions(+), 3 deletions(-)

diff --git a/src/Calculators/Thermovec.jl b/src/Calculators/Thermovec.jl
index 2d66b7c98..6cd0bd7b7 100644
--- a/src/Calculators/Thermovec.jl
+++ b/src/Calculators/Thermovec.jl
@@ -74,7 +74,7 @@ end
 end
 export selectPoly
 
-@inline function calcHSCpdless(poly::NASApolynomialvec,T::X) where {X<:Real}
+@inline function calcHSCpdless(poly::NASApolynomialvec{Array{1,EmptyThermoCoverageDependence}},T::X;coverages=nothing) where {X<:Real}
     if size(poly.coefs)[1] != 7
         Tpoly0 = T
         Tpoly1 = T*T
@@ -119,9 +119,63 @@ export selectPoly
     return (cpdivR,hdivRT,sdivR)
 end
 
-@inline function calcHSCpdless(nasavec::NASAvec,T::X) where {X<:Real}
+@inline function calcHSCpdless(poly::NASApolynomialvec{Q},T::X;coverages=nothing) where {X<:Real,Q}
+    hdivRT = similar(coverages)
+    sdivR = similar(coverages)
+    cpdivR = similar(coverages)
+    if size(poly.coefs)[1] != 7
+        Tpoly0 = T
+        Tpoly1 = T*T
+        Tpoly2 = Tpoly1*T
+        Tpoly3 = Tpoly2*T
+        Tpoly4 = 1.0/T
+        Tpoly5 = Tpoly4*Tpoly4
+        Tpoly6 = log(T)
+
+        ct0 = poly.coefs[1,:].*Tpoly5
+        ct1 = poly.coefs[2,:].*Tpoly4
+        ct2 = poly.coefs[3,:]
+        ct3 = poly.coefs[4,:].*Tpoly0
+        ct4 = poly.coefs[5,:].*Tpoly1
+        ct5 = poly.coefs[6,:].*Tpoly2
+        ct6 = poly.coefs[7,:].*Tpoly3
+
+        cpdivR .= ct0 .+ ct1 .+ ct2 .+ ct3 .+ ct4 .+ ct5 .+ ct6
+        hdivRT .= -ct0.+Tpoly6.*ct1.+ct2.+0.5.*ct3.+0.33333333333.*ct4.+0.25.*ct5+0.2.*ct6+poly.coefs[8,:].*Tpoly4
+        sdivR .= -0.5.*ct0 .- ct1 .+ Tpoly6.*ct2 .+ ct3 .+ 0.5.*ct4 .+ 0.33333333333.*ct5 .+ 0.25*ct6 .+ poly.coefs[9,:]
+    else
+        Tpoly0 = T
+        Tpoly1 = T*T
+        Tpoly2 = Tpoly1*T
+        Tpoly3 = Tpoly2*T
+        Tpoly4 = 1.0/T
+        #Tpoly5 = Tpoly4*Tpoly4
+        Tpoly6 = log(T)
+
+        #ct0 = 0.0
+        #ct1 = 0.0
+        ct2 = poly.coefs[1,:]
+        ct3 = poly.coefs[2,:].*Tpoly0
+        ct4 = poly.coefs[3,:].*Tpoly1
+        ct5 = poly.coefs[4,:].*Tpoly2
+        ct6 = poly.coefs[5,:].*Tpoly3
+
+        cpdivR .= ct2 .+ ct3 .+ ct4 .+ ct5 .+ ct6
+        hdivRT .= ct2.+0.5.*ct3.+0.33333333333.*ct4.+0.25.*ct5.+0.2.*ct6.+poly.coefs[6,:].*Tpoly4
+        sdivR .= Tpoly6.*ct2 .+ ct3 .+ 0.5.*ct4 .+ 0.33333333333.*ct5 .+ 0.25.*ct6 .+ poly.coefs[7,:]
+    end
+    for i in 1:length(hdivRT)
+        hdivRT .+= getcovdepenthalpycorrection(poly.covdeps[i],T,coverages)/(R*T)
+        sdivR .+= getcovdepentropycorrection(poly.covdeps[i],T,coverages)/R 
+        cpdivR .+= getcovdepheatcapacitycorrection(poly.covdeps[i],T,coverages)/R
+    end
+    return (cpdivR,hdivRT,sdivR)
+end
+
+
+@inline function calcHSCpdless(nasavec::NASAvec,T::X;coverages=nothing) where {X<:Real}
     poly = selectPoly(nasavec,T)
-    return calcHSCpdless(poly,T)
+    return calcHSCpdless(poly,T;coverages=coverages)
 end
 
 @with_kw struct ConstantGvec{B,J} <: AbstractThermo

From 7ed25d32ba1906c79f68fcb641b824a167f87f08 Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Thu, 5 Jun 2025 12:19:07 -0700
Subject: [PATCH 07/28] adapt vectorized rate coefficient computation to handle
 coverage inputs

---
 src/Calculators/Ratevec.jl | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/Calculators/Ratevec.jl b/src/Calculators/Ratevec.jl
index 7b5045743..ec65c8af9 100644
--- a/src/Calculators/Ratevec.jl
+++ b/src/Calculators/Ratevec.jl
@@ -10,7 +10,7 @@ export AbstractRatevec
         n::K
         Ea::Q
 end
-function Arrheniusvec(arrs::T) where {T<:AbstractArray}
+function Arrheniusvec(arrs)
     A = zeros(length(arrs)) 
     n = zeros(length(arrs)) 
     Ea = zeros(length(arrs))
@@ -21,8 +21,8 @@ function Arrheniusvec(arrs::T) where {T<:AbstractArray}
     end
     return Arrheniusvec(A=A,n=n,Ea=Ea)
 end
-@inline (arr::Arrheniusvec)(;T::Q,P::N=0.0,C::S=0.0,phi=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath @inbounds arr.A.*exp.(arr.n.*log(T).-arr.Ea.*(1.0/(R*T)))
-@inline (arr::Arrheniusvec)(T::Q;P::N=0.0,C::S=0.0,phi=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath @inbounds arr.A.*exp.(arr.n.*log(T).-arr.Ea.*(1.0/(R*T)))
+@inline (arr::Arrheniusvec)(;T::Q,P::N=0.0,C::S=0.0,phi=0.0,coverages=nothing) where {Q<:Real,N<:Real,S<:Real} = @fastmath @inbounds arr.A.*exp.(arr.n.*log(T).-arr.Ea.*(1.0/(R*T)))
+@inline (arr::Arrheniusvec)(T::Q;P::N=0.0,C::S=0.0,phi=0.0,coverages=nothing) where {Q<:Real,N<:Real,S<:Real} = @fastmath @inbounds arr.A.*exp.(arr.n.*log(T).-arr.Ea.*(1.0/(R*T)))
 export Arrheniusvec
 
 @with_kw struct Chebyshevvec{T<:AbstractArray,Q<:Real,S<:Real,V<:Real,B<:Real} <: AbstractRate
@@ -86,7 +86,7 @@ export getredtemp
 end
 export getredpress
 
-@inline function (ch::Chebyshevvec)(;T::N,P::Q=0.0,C::B=0.0,phi=0.0) where {N<:Real,B<:Real,Q<:Real}
+@inline function (ch::Chebyshevvec)(;T::N,P::Q=0.0,C::B=0.0,phi=0.0,coverages=nothing) where {N<:Real,B<:Real,Q<:Real}
     k = zeros(N,size(ch.coefs)[1])
     Tred = getredtemp(ch,T)
     Pred = getredpress(ch,P)
@@ -170,7 +170,7 @@ function Troevec(troes::T) where {T<:AbstractArray}
 end
 export Troevec
     
-@inline function (tr::Troevec)(;T::Q=nothing,P::R=0.0,C::S=nothing,phi=0.0) where {Q<:Real,R<:Real,S<:Real}
+@inline function (tr::Troevec)(;T::Q=nothing,P::R=0.0,C::S=nothing,phi=0.0,coverages=nothing) where {Q<:Real,R<:Real,S<:Real}
     k0 = tr.arrlow(T=T)
     kinf = tr.arrhigh(T=T)
     @fastmath Pr = k0.*C./kinf
@@ -201,7 +201,7 @@ function PdepArrheniusvec(pdeparrs::T) where {T<:AbstractArray}
 end
 export PdepArrheniusvec
 
-@inline function (parr::PdepArrheniusvec)(;T::Q=nothing,P::V=nothing,C::S=0.0,phi=0.0) where {Q<:Real,V<:Real,S<:Real}
+@inline function (parr::PdepArrheniusvec)(;T::Q=nothing,P::V=nothing,C::S=0.0,phi=0.0,coverages=nothing) where {Q<:Real,V<:Real,S<:Real}
     inds = getBoundingIndsSorted(P,parr.Ps)::Tuple{Int64,Int64}
     if inds[2] == -1
         return @inbounds parr.arrvecs[inds[1]](T=T)

From 40948037d5c411bf52c0d81102e68cabfd7aa579 Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Thu, 5 Jun 2025 12:30:31 -0700
Subject: [PATCH 08/28] assign covdep indices in getreactionindices

When we put information in the mechanism file we don't yet know the index associated with the species coverage dependence is with respect to
this fills in the index (as opposed to species name) dictionaries on the coverage dependence objects so they are ready for efficient evaluation
---
 src/Phase.jl | 38 +++++++++++++++++++++++++++++++++++++-
 1 file changed, 37 insertions(+), 1 deletion(-)

diff --git a/src/Phase.jl b/src/Phase.jl
index 784ed5c80..5165f7456 100644
--- a/src/Phase.jl
+++ b/src/Phase.jl
@@ -317,9 +317,23 @@ export iterate
 Broadcast.broadcastable(p::T) where {T<:AbstractPhase} = Ref(p)
 export broadcastable
 
-function getreactionindices(spcs,rxns) where {Q<:AbstractPhase}
+function getreactionindices(spcs,rxns)
     arr = zeros(Int64,(8,length(rxns)))
     names = [spc.name for spc in spcs]
+    for spc in spcs
+        if hasproperty(spc.thermo,:covdep)
+            if spc.thermo.covdep isa PolynomialThermoEnergyCoverageDependence
+                for (name,v) in spc.thermo.covdep.polys
+                    ind = findfirst(isequal(name),names)
+                    spc.thermo.covdep.indpolys[ind] = v
+                end 
+            elseif spc.thermo.covdep isa EmptyThermoCoverageDependence
+
+            else 
+                throw(TypeError(spc.thermo.covdep,"Thermo Coverage Dependence Type Not Understood"))
+            end
+        end
+    end
     for (i,rxn) in enumerate(rxns)
         inds = [findfirst(isequal(spc),spcs) for spc in rxn.reactants]
         for (j,spc) in enumerate(rxn.reactants)
@@ -343,6 +357,28 @@ function getreactionindices(spcs,rxns) where {Q<:AbstractPhase}
                 end
             end
         end
+        if hasproperty(rxn.kinetics,:covdep) 
+            if rxn.kinetics.covdep isa PolynomialRateCoverageDependence
+                for (name,v) in rxn.kinetics.covdep.Epolys
+                    ind = findfirst(isequal(name),names)
+                    rxn.kinetics.covdep.indEpolys[ind] = v
+                end 
+                for (name,v) in rxn.kinetics.covdep.avals
+                    ind = findfirst(isequal(name),names)
+                    rxn.kinetics.covdep.indavals[ind] = v
+                end 
+                for (name,v) in rxn.kinetics.covdep.ms
+                    ind = findfirst(isequal(name),names)
+                    rxn.kinetics.covdep.indms[ind] = v
+                end 
+            elseif rxn.kinetics.covdep isa EmptyRateCoverageDependence
+
+            else 
+                throw(TypeError(rxn.kinetics.covdep,"Kinetic Coverage Dependence Type Not Understood"))
+            end
+
+                
+        end
     end
     return arr
 end

From 0710f7ede31734314a1b97aebd8f05dcdc39f5b5 Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Thu, 5 Jun 2025 12:31:49 -0700
Subject: [PATCH 09/28] adapt calculations of phase thermochemistry and
 kinetics to handle coverage dependence

---
 src/PhaseState.jl | 52 +++++++++++++++++++++++------------------------
 1 file changed, 26 insertions(+), 26 deletions(-)

diff --git a/src/PhaseState.jl b/src/PhaseState.jl
index ff2683e8f..821136447 100644
--- a/src/PhaseState.jl
+++ b/src/PhaseState.jl
@@ -5,23 +5,23 @@ using Tracker
 using ReverseDiff
 using RecursiveArrayTools
 
-@inline function calcgibbs(ph::U,T::W) where {U<:IdealPhase,W<:Real}
-    return getGibbs.(getfield.(ph.species,:thermo),T)
+@inline function calcgibbs(ph::U,T::W;coverages=nothing) where {U<:IdealPhase,W<:Real}
+    return getGibbs.(getfield.(ph.species,:thermo),T,coverages=coverages)
 end
 export calcgibbs
 
-@inline function calcenthalpyinternalgibbs(ph::U,T::W,P::Z,V::Q) where {U<:IdealPhase,W,Z,Q<:Real}
-    Hs = getEnthalpy.(getfield.(ph.species,:thermo),T)
+@inline function calcenthalpyinternalgibbs(ph::U,T::W,P::Z,V::Q;coverages=nothing) where {U<:IdealPhase,W,Z,Q<:Real}
+    Hs = getEnthalpy.(getfield.(ph.species,:thermo),T,coverages=coverages)
     Us = Hs .- P*V
-    Gs = Hs .- T.*getEntropy.(getfield.(ph.species,:thermo),T)
+    Gs = Hs .- T.*getEntropy.(getfield.(ph.species,:thermo),T,coverages=coverages)
     return Hs,Us,Gs
 end
 
 
-@inline function calcenthalpyinternalgibbs(ph::Union{IdealGas,IdealSurface},T::W,P::Z,V::Q) where {W,Z,Q<:Real}
-    Hs = getEnthalpy.(getfield.(ph.species,:thermo),T)
+@inline function calcenthalpyinternalgibbs(ph::Union{IdealGas,IdealSurface},T::W,P::Z,V::Q;coverages=nothing) where {W,Z,Q<:Real}
+    Hs = getEnthalpy.(getfield.(ph.species,:thermo),T,coverages)
     Us = Hs .- R*T
-    Gs = Hs .- T.*getEntropy.(getfield.(ph.species,:thermo),T)
+    Gs = Hs .- T.*getEntropy.(getfield.(ph.species,:thermo),T,coverages)
     return Hs,Us,Gs
 end
 
@@ -50,27 +50,27 @@ end
 
 export makespcsvector
 
-@inline function getkf(rxn::ElementaryReaction,ph,T,P,C,ns,V,phi)
+@inline function getkf(rxn::ElementaryReaction,ph,T,P,C,ns,V,phi;coverages=nothing)
     if isdefined(rxn.kinetics,:efficiencies) && length(rxn.kinetics.efficiencies) > 0
         @views @inbounds @fastmath C += sum([ns[i]*val for (i,val) in rxn.kinetics.efficiencies])/V
     end
-    return rxn.kinetics(T=T,P=P,C=C,phi=phi)
+    return rxn.kinetics(T=T,P=P,C=C,phi=phi,coverages=coverages)
 end
 export getkf
 
-@inline function getkfs(ph::U,T::W1,P::W2,C::W3,ns::Q,V::W4,phi) where {U,W1,W2,W3,W4<:Real,Q<:AbstractArray}
+@inline function getkfs(ph::U,T::W1,P::W2,C::W3,ns::Q,V::W4,phi;coverages=nothing) where {U,W1,W2,W3,W4<:Real,Q<:AbstractArray}
     kfs = similar(ns,length(ph.reactions))
     i = 1
     oldind = 1
     ind = 0
     while i <= length(ph.veckineticsinds) #vectorized kinetics
         @inbounds ind = ph.veckineticsinds[i]
-        @inbounds kfs[oldind:ind] = ph.veckinetics[i](;T=T,P=P,C=C)
+        @inbounds kfs[oldind:ind] = ph.veckinetics[i](;T=T,P=P,C=C,coverages=coverages)
         oldind = ind+1
         i += 1
     end
     @simd for i in ind+1:length(ph.reactions)
-        @inbounds kfs[i] = getkf(ph.reactions[i],ph,T,P,C,ns,V,phi)
+        @inbounds kfs[i] = getkf(ph.reactions[i],ph,T,P,C,ns,V,phi,coverages=coverages)
     end
     return kfs
 end
@@ -181,12 +181,12 @@ export getKcs
 Calculates the forward and reverse rate coefficients for a given reaction, phase and state
 Maintains diffusion limitations if the phase has diffusionlimited=true
 """
-@inline function getkfkrev(rxn::ElementaryReaction,ph::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Q2,diffs::Q3,V::W5,phi::W8;kf::W6=-1.0,f::W7=-1.0) where {U<:AbstractPhase,W8,W6,W7,W5,W4,W1,W2,W3<:Real,Q1,Q2,Q3<:AbstractArray}
+@inline function getkfkrev(rxn::ElementaryReaction,ph::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Q2,diffs::Q3,V::W5,phi::W8;coverages=nothing,kf::W6=-1.0,f::W7=-1.0) where {U<:AbstractPhase,W8,W6,W7,W5,W4,W1,W2,W3<:Real,Q1,Q2,Q3<:AbstractArray}
     if signbit(kf) 
         if signbit(f)
-            kf = getkf(rxn,ph,T,P,C,ns,V,phi)
+            kf = getkf(rxn,ph,T,P,C,ns,V,phi;coverages=coverages)
         else
-            kf = getkf(rxn,ph,T,P,C,ns,V,phi)*f
+            kf = getkf(rxn,ph,T,P,C,ns,V,phi;coverages=coverages)*f
         end
     end
     Kc = getKc(rxn,ph,T,Gs,phi)
@@ -223,7 +223,7 @@ Maintains diffusion limitations if the phase has diffusionlimited=true
 end
 export getkfkrev
 
-@inline function getkfkrevs(phase::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Q2,diffs::Q3,V::W5,phi::W7;kfs::W6=nothing) where {U,W7,W6,W5<:Real,W1<:Real,W2<:Real,W3,W4,Q1<:AbstractArray,Q2,Q3<:AbstractArray}
+@inline function getkfkrevs(phase::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Q2,diffs::Q3,V::W5,phi::W7;coverages=nothing,kfs::W6=nothing) where {U,W7,W6,W5<:Real,W1<:Real,W2<:Real,W3,W4,Q1<:AbstractArray,Q2,Q3<:AbstractArray}
     if !phase.diffusionlimited && kfs === nothing
         kfs = getkfs(phase,T,P,C,ns,V,phi)
         if phi == 0.0
@@ -241,14 +241,14 @@ export getkfkrev
         len = length(phase.reactions)
         krev = zeros(typeof(N),len)
         @simd for i = 1:len
-           @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi;kf=kfs[i])
+           @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi;coverages=coverages,kf=kfs[i])
         end
     else
         len = length(phase.reactions)
         kfs = zeros(typeof(N),len)
         krev = zeros(typeof(N),len)
         @simd for i = 1:len
-           @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi)
+           @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi;coverages=coverages)
         end
     end
     kfs .*= phase.forwardability
@@ -256,7 +256,7 @@ export getkfkrev
     return kfs,krev
 end
 
-@inline function getkfkrevs(phase::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Q2,diffs::Q3,V::W5,phi::W7;kfs::W6=nothing) where {U,W7,W6,W5<:Real,W1<:Real,W2<:Real,W3,W4,Q1<:AbstractArray,Q2<:Union{ReverseDiff.TrackedArray,Tracker.TrackedArray},Q3<:AbstractArray} #autodiff p
+@inline function getkfkrevs(phase::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Q2,diffs::Q3,V::W5,phi::W7;coverages=nothing,kfs::W6=nothing) where {U,W7,W6,W5<:Real,W1<:Real,W2<:Real,W3,W4,Q1<:AbstractArray,Q2<:Union{ReverseDiff.TrackedArray,Tracker.TrackedArray},Q3<:AbstractArray} #autodiff p
     if !phase.diffusionlimited && kfs === nothing
         kfs = getkfs(phase,T,P,C,ns,V,phi)
         if phi == 0.0
@@ -275,7 +275,7 @@ end
         krev = similar(kfs)
         kfsdiff = similar(kfs)
         @simd for i = 1:len
-           @fastmath @inbounds kfsdiff[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi;kf=kfs[i])
+           @fastmath @inbounds kfsdiff[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi;coverages=coverages,kf=kfs[i])
         end
         return kfsdiff, krev
     else
@@ -283,15 +283,15 @@ end
         kfs = zeros(typeof(Gs[1]),len)
         krev = zeros(typeof(Gs[1]),len)
         @simd for i = 1:len
-           @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi)
+           @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi;coverages=coverages)
         end
     end
     return kfs,krev
 end
 
-@inline function getkfkrevs(phase::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Array{Q2,1},diffs::Q3,V::W5,phi::W7;kfs::W6=nothing) where {U,W7,W6,W5<:Real,W1<:Real,W2<:Real,W3,W4,Q1<:AbstractArray,Q2<:ForwardDiff.Dual,Q3<:AbstractArray} #autodiff p
+@inline function getkfkrevs(phase::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Array{Q2,1},diffs::Q3,V::W5,phi::W7;coverages=nothing,kfs::W6=nothing) where {U,W7,W6,W5<:Real,W1<:Real,W2<:Real,W3,W4,Q1<:AbstractArray,Q2<:ForwardDiff.Dual,Q3<:AbstractArray} #autodiff p
     if !phase.diffusionlimited && kfs === nothing
-        kfs = getkfs(phase,T,P,C,ns,V,phi)
+        kfs = getkfs(phase,T,P,C,ns,V,phi;coverages=coverages)
         if phi == 0.0
             krev = @fastmath kfs./getKcs(phase,T,Gs)
         else 
@@ -307,14 +307,14 @@ end
         len = length(phase.reactions)
         krev = similar(kfs)
         @simd for i = 1:len
-           @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi;kf=kfs[i])
+           @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi;coverages=coverages,kf=kfs[i])
         end
     else
         len = length(phase.reactions)
         kfs = zeros(typeof(Gs[1]),len)
         krev = zeros(typeof(Gs[1]),len)
         @simd for i = 1:len
-           @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi)
+           @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi;coverages=coverages)
         end
     end
     kfs .*= phase.forwardability

From 8840e859bd20ddc1217b8c7e08e78ae65b43ac79 Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Thu, 5 Jun 2025 12:33:45 -0700
Subject: [PATCH 10/28] adapt ConstantTAPhiDomain to handle coverage dependence

---
 src/Domain.jl | 17 +++++++++++++----
 1 file changed, 13 insertions(+), 4 deletions(-)

diff --git a/src/Domain.jl b/src/Domain.jl
index 82d2e05ea..c2477c86d 100644
--- a/src/Domain.jl
+++ b/src/Domain.jl
@@ -656,6 +656,7 @@ mutable struct ConstantTAPhiDomain{N<:AbstractPhase,S<:Integer,W<:Real,W2<:Real,
     t::MArray{Tuple{1},W2,1,1}
     p::Array{W,1}
     thermovariabledict::Dict{String,Int64}
+    iscovdep::Bool
 end
 function ConstantTAPhiDomain(; phase::E2, initialconds::Dict{X,X2}, constantspecies::Array{X3,1}=Array{String,1}(),
     sparse::Bool=false, sensitivity::Bool=false, stationary::Bool=false) where {E<:Real,E2<:AbstractPhase,W<:Real,X,X2,X3}
@@ -683,6 +684,7 @@ function ConstantTAPhiDomain(; phase::E2, initialconds::Dict{X,X2}, constantspec
     @assert T != 0.0
     ns = y0
     N = sum(ns)
+    coverages = ns ./ phase.sitedensity
 
     if length(constantspecies) > 0
         spcnames = getfield.(phase.species, :name)
@@ -691,15 +693,22 @@ function ConstantTAPhiDomain(; phase::E2, initialconds::Dict{X,X2}, constantspec
         constspcinds = Array{Int64,1}()
     end
     efficiencyinds = [rxn.index for rxn in phase.reactions if typeof(rxn.kinetics) <: AbstractFalloffRate && length(rxn.kinetics.efficiencies) > 0]
-    Gs = calcgibbs(phase, T)
+    Gs = calcgibbs(phase, T; coverages=coverages)
     if :solvent in fieldnames(typeof(phase)) && typeof(phase.solvent) != EmptySolvent
         mu = phase.solvent.mu(T)
     else
         mu = 0.0
     end
     C = 0.0 #this currently shouldn't matter here, on a surface you shouldn't have pdep
-    kfs, krevs = getkfkrevs(phase, T, 0.0, C, N, ns, Gs, [], A, phi)
-    p = vcat(Gs, kfs)
+    iscovdep = any([!(spc.thermo.covdep isa EmptyThermoCoverageDependence) for spc in phase.species]) || any([!(rxn.kinetics.covdep isa EmptyRateCoverageDependence) for rxn in phase.reactions])
+    kfs, krevs = getkfkrevs(phase, T, 0.0, C, N, ns, Gs, [], A, phi; coverages=coverages)
+    if !iscovdep
+        p = vcat(Gs, kfs)
+        alternativepformat = false 
+    else 
+        p = vcat(zeros(length(phase.species)), ones(length(phase.reactions)))
+        alternativepformat = true 
+    end 
     if sparse
         jacobian = spzeros(typeof(T), length(phase.species), length(phase.species))
     else
@@ -707,7 +716,7 @@ function ConstantTAPhiDomain(; phase::E2, initialconds::Dict{X,X2}, constantspec
     end
     rxnarray = getreactionindices(phase)
     return ConstantTAPhiDomain(phase, [1, length(phase.species)], [1, length(phase.species) + length(phase.reactions)], constspcinds,
-        T, A, phi, kfs, krevs, efficiencyinds, Gs, rxnarray, mu, Array{Float64,1}(), jacobian, sensitivity, false, MVector(false), MVector(0.0), p, Dict{String,Int64}()), y0, p
+        T, A, phi, kfs, krevs, efficiencyinds, Gs, rxnarray, mu, Array{Float64,1}(), jacobian, sensitivity, alternativepformat, MVector(false), MVector(0.0), p, Dict{String,Int64}(), iscovdep), y0, p
 end
 export ConstantTAPhiDomain
 

From e302296cbef74b50c5b803812f3c9a9ff46a7894 Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Thu, 5 Jun 2025 12:36:53 -0700
Subject: [PATCH 11/28] adapt general domains to outputting coverages from
 calcthermo

also ensures gibbs free energy is output (before an empty vector was returned if it was not expected to be needed, but coverage dependence makes us needed at the interface even for domains where Gs is constant)
---
 src/Domain.jl | 108 +++++++++++++++++++++++++-------------------------
 1 file changed, 54 insertions(+), 54 deletions(-)

diff --git a/src/Domain.jl b/src/Domain.jl
index c2477c86d..f802fbc80 100644
--- a/src/Domain.jl
+++ b/src/Domain.jl
@@ -931,7 +931,7 @@ export ConstantTLiqFilmDomain
     for ind in d.efficiencyinds #efficiency related rates may have changed
         d.kfs[ind], d.krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, d.Gs, d.diffusivity, V, 0.0)
     end
-    return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+    return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
 end
 
 @inline function calcthermo(d::ConstantTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealGas,Y<:Integer,J<:Array{Float64,1},Q<:Float64} #uses parameter input
@@ -952,13 +952,13 @@ end
         for ind in d.efficiencyinds #efficiency related rates may have changed
             d.kfs[ind], d.krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; f=kfps[ind])
         end
-        return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+        return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
     elseif nothermochg
         d.kfs = kfps
         for ind in d.efficiencyinds #efficiency related rates may have changed
             d.kfs[ind], d.krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; f=kfps[ind])
         end
-        return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+        return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
     else #need to handle thermo changes
         d.kfs .= kfps
         if !d.alternativepformat
@@ -970,7 +970,7 @@ end
         for ind in d.efficiencyinds #efficiency related rates may have changed
             d.kfs[ind], d.krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; f=kfps[ind])
         end
-        return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+        return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
     end
 end
 
@@ -991,7 +991,7 @@ end
     for ind in d.efficiencyinds #efficiency related rates may have changed
         kfs[ind], krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0; f=kfs[ind])
     end
-    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
 end
 
 @inline function calcthermo(d::ConstantTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J,Q} #Autodiff p
@@ -1011,7 +1011,7 @@ end
     for ind in d.efficiencyinds #efficiency related rates may have changed
         kfs[ind], krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0; f=kfs[ind])
     end
-    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
 end
 
 @inline function calcthermo(d::ConstantTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J<:Union{ReverseDiff.TrackedArray,Tracker.TrackedArray},Q} #Autodiff p
@@ -1033,7 +1033,7 @@ end
     for ind in d.efficiencyinds #efficiency related rates may have changed
         kfs[ind], krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0; f=kfs[ind])
     end
-    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
 end
 
 @inline function calcthermo(d::ConstantTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Union{ReverseDiff.TrackedArray,Tracker.TrackedArray},W<:IdealGas,Y<:Integer,J,Q} #Tracker/reversediff
@@ -1050,7 +1050,7 @@ end
         kfs = [ind in d.efficiencyinds ? getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0)[1] * d.p[length(d.phase.species)+ind] * p[d.parameterindexes[1]-1+length(d.phase.species)+ind] : p[d.parameterindexes[1]-1+length(d.phase.species)+ind] for ind in 1:length(d.phase.reactions)]
     end
     krevs = getkfkrevs(d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0; kfs=kfs)[2]
-    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
 end
 
 @inline function calcthermo(d::ConstantTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Union{ReverseDiff.TrackedArray,Tracker.TrackedArray},W<:IdealGas,Y<:Integer,J<:Union{ReverseDiff.TrackedArray,Tracker.TrackedArray},Q} #Tracker/reversediff
@@ -1067,7 +1067,7 @@ end
         kfs = [ind in d.efficiencyinds ? getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0)[1] * d.p[length(d.phase.species)+ind] * p[d.parameterindexes[1]-1+length(d.phase.species)+ind] : p[d.parameterindexes[1]-1+length(d.phase.species)+ind] for ind in 1:length(d.phase.reactions)]
     end
     krevs = getkfkrevs(d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0; kfs=kfs)[2]
-    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
 end
 
 @inline function calcthermo(d::ConstantVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:SciMLBase.NullParameters,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1092,7 +1092,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, d.V, 0.0)
-    return ns, cs, T, P, d.V, C, N, 0.0, kfs, krevs, Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0
+    return ns, cs, T, P, d.V, C, N, 0.0, kfs, krevs, Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0, nothing
 end
 
 @inline function calcthermo(d::ConstantVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1118,7 +1118,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, d.V, 0.0)
-    return @views @fastmath ns, cs, T, P, d.V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0
+    return @views @fastmath ns, cs, T, P, d.V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0, nothing
 end
 
 @inline function calcthermo(d::ConstantVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1144,7 +1144,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, d.V, 0.0)
-    return @views @fastmath ns, cs, T, P, d.V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0
+    return @views @fastmath ns, cs, T, P, d.V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0, nothing
 end
 
 @inline function calcthermo(d::ConstantPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:SciMLBase.NullParameters,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1168,7 +1168,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, d.P, C, N, ns, Gs, diffs, V, 0.0)
-    return ns, cs, T, d.P, V, C, N, 0.0, kfs, krevs, Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0
+    return ns, cs, T, d.P, V, C, N, 0.0, kfs, krevs, Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, nothing
 end
 
 @inline function calcthermo(d::ConstantPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1194,9 +1194,9 @@ end
     end
     kfs, krevs = getkfkrevs(d.phase, T, d.P, C, N, ns, Gs, diffs, V, 0.0)
     if p != SciMLBase.NullParameters()
-        return @views @fastmath ns, cs, T, d.P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0
+        return @views @fastmath ns, cs, T, d.P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, nothing
     else
-        return ns, cs, T, d.P, V, C, N, 0.0, kfs, krevs, Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0
+        return ns, cs, T, d.P, V, C, N, 0.0, kfs, krevs, Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, nothing
     end
 end
 
@@ -1222,7 +1222,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, d.P, C, N, ns, Gs, diffs, V, 0.0)
-    return @views @fastmath ns, cs, T, d.P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0
+    return @views @fastmath ns, cs, T, d.P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, nothing
 end
 
 @inline function calcthermo(d::ParametrizedVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:SciMLBase.NullParameters,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1247,7 +1247,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
-    return ns, cs, T, P, V, C, N, 0.0, kfs, krevs, Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0
+    return ns, cs, T, P, V, C, N, 0.0, kfs, krevs, Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0, nothing
 end
 
 @inline function calcthermo(d::ParametrizedVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1273,7 +1273,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
-    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0
+    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0, nothing
 end
 
 @inline function calcthermo(d::ParametrizedVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1299,7 +1299,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
-    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0
+    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0, nothing
 end
 
 @inline function calcthermo(d::ParametrizedPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:SciMLBase.NullParameters,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1324,7 +1324,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
-    return ns, cs, T, P, V, C, N, 0.0, kfs, krevs, Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0
+    return ns, cs, T, P, V, C, N, 0.0, kfs, krevs, Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, nothing
 end
 
 @inline function calcthermo(d::ParametrizedPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1350,7 +1350,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
-    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0
+    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, nothing
 end
 @inline function calcthermo(d::ParametrizedPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
     P = d.P(t)
@@ -1375,7 +1375,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
-    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0
+    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, nothing
 end
 
 @inline function calcthermo(d::ParametrizedTConstantVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:SciMLBase.NullParameters,W<:IdealDiluteSolution,Y<:Integer,J<:AbstractArray,Q}
@@ -1396,7 +1396,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, d.phi)
-    return ns, cs, T, P, V, C, N, mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0
+    return ns, cs, T, P, V, C, N, mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0, nothing
 end
 
 @inline function calcthermo(d::ParametrizedTConstantVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealDiluteSolution,Y<:Integer,J<:AbstractArray,Q}
@@ -1418,7 +1418,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, d.phi)
-    return @views @fastmath ns, cs, T, P, V, C, N, mu, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), d.phi
+    return @views @fastmath ns, cs, T, P, V, C, N, mu, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), d.phi, nothing
 end
 
 @inline function calcthermo(d::ParametrizedTConstantVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealDiluteSolution,Y<:Integer,J<:AbstractArray,Q}
@@ -1440,7 +1440,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, d.phi)
-    return @views @fastmath ns, cs, T, P, V, C, N, mu, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), d.phi
+    return @views @fastmath ns, cs, T, P, V, C, N, mu, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), d.phi, nothing
 end
 
 @inline function calcthermo(d::ParametrizedTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:SciMLBase.NullParameters,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1465,7 +1465,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
-    return ns, cs, T, P, V, C, N, 0.0, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0
+    return ns, cs, T, P, V, C, N, 0.0, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0, nothing
 end
 
 @inline function calcthermo(d::ParametrizedTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1491,7 +1491,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
-    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0
+    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0, nothing
 end
 
 @inline function calcthermo(d::ParametrizedTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1517,7 +1517,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
-    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0
+    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0, nothing
 end
 
 @inline function calcthermo(d::ConstantTVDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2<:SciMLBase.NullParameters,W<:IdealDiluteSolution,Y<:Integer,J<:AbstractArray,Q}
@@ -1526,7 +1526,7 @@ end
     cs = ns ./ d.V
     C = N / d.V
     P = 1.0e8
-    return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+    return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
 end
 
 @inline function calcthermo(d::ConstantTVDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2<:Array{Float64,1},W<:IdealDiluteSolution,Y<:Integer,J<:Array{Float64,1},Q}
@@ -1539,31 +1539,31 @@ end
         @views nothermochg = d.Gs == p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
         @views nokfchg = d.kfsnondiff == p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
         if nothermochg && nokfchg
-            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
         elseif nothermochg
             d.kfsnondiff = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, d.V, d.phi; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
         else
             d.kfsnondiff = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, d.V, d.phi; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
         end
     else
         @views nothermochg = d.Gs == d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
         @views nokfchg = d.kfsnondiff == d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
         if nothermochg && nokfchg
-            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
         elseif nothermochg
             d.kfsnondiff .= d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, d.V, d.phi; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
         else
             d.kfsnondiff .= d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.Gs .= d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, d.V, d.phi; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
         end
     end
 
@@ -1583,7 +1583,7 @@ end
         kfsnondiff = convert(typeof(y), d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)])
     end
     kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, d.V, d.phi; kfs=kfsnondiff)
-    return ns, cs, d.T, P, d.V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+    return ns, cs, d.T, P, d.V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
 end
 
 @inline function calcthermo(d::ConstantTVDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2,W<:IdealDiluteSolution,Y<:Integer,J<:AbstractArray,Q} #autodiff p
@@ -1600,7 +1600,7 @@ end
         kfsnondiff = d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
     end
     kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, d.V, d.phi; kfs=kfsnondiff)
-    return ns, cs, d.T, P, d.V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+    return ns, cs, d.T, P, d.V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
 end
 
 @inline function calcthermo(d::ConstantTAPhiDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2<:SciMLBase.NullParameters,W<:IdealSurface,Y<:Integer,J<:AbstractArray,Q}
@@ -1683,7 +1683,7 @@ function calcthermo(d::FragmentBasedConstantTrhoDomain{W,Y}, y::J, t::Q, p::Q2=S
     cs = ns ./ V
     C = N / V
     P = 1.0e8
-    return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+    return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
 end
 
 function calcthermo(d::FragmentBasedConstantTrhoDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2<:Array{Float64,1},W<:FragmentBasedIdealFilm,Y<:Integer,J<:Array{Float64,1},Q}
@@ -1698,31 +1698,31 @@ function calcthermo(d::FragmentBasedConstantTrhoDomain{W,Y}, y::J, t::Q, p::Q2=S
         @views nothermochg = d.Gs == p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
         @views nokfchg = d.kfsnondiff == p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
         if nothermochg && nokfchg
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
         elseif nothermochg
             d.kfsnondiff = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
         else
             d.kfsnondiff = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
         end
     else
         @views nothermochg = d.Gs == d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
         @views nokfchg = d.kfsnondiff == d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
         if nothermochg && nokfchg
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
         elseif nothermochg
             d.kfsnondiff .= d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
         else
             d.kfsnondiff .= d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.Gs .= d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
         end
     end
 
@@ -1744,7 +1744,7 @@ function calcthermo(d::FragmentBasedConstantTrhoDomain{W,Y}, y::Array{W2,1}, t::
         kfsnondiff = convert(typeof(y), d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)])
     end
     kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, V, 0.0; kfs=kfsnondiff)
-    return ns, cs, d.T, P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+    return ns, cs, d.T, P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
 end
 
 function calcthermo(d::FragmentBasedConstantTrhoDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2,W<:FragmentBasedIdealFilm,Y<:Integer,J<:AbstractArray,Q} #autodiff p
@@ -1763,7 +1763,7 @@ function calcthermo(d::FragmentBasedConstantTrhoDomain{W,Y}, y::J, t::Q, p::Q2=S
         kfsnondiff = d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
     end
     kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, V, 0.0; kfs=kfsnondiff)
-    return ns, cs, d.T, P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+    return ns, cs, d.T, P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
 end
 
 @inline function calcthermo(d::ConstantTLiqFilmDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2<:Array{Float64,1},W<:IdealDiluteSolution,Y<:Integer,J<:Array{Float64,1},Q}
@@ -1777,31 +1777,31 @@ end
         @views nothermochg = d.Gs == p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
         @views nokfchg = d.kfsnondiff == p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
         if nothermochg && nokfchg
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
         elseif nothermochg
             d.kfsnondiff = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, d.phi; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
         else
             d.kfsnondiff = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, d.phi; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
         end
     else
         @views nothermochg = d.Gs == d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
         @views nokfchg = d.kfsnondiff == d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
         if nothermochg && nokfchg
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
         elseif nothermochg
             d.kfsnondiff .= d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, d.phi; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
         else
             d.kfsnondiff .= d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.Gs .= d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, d.phi; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
         end
     end
 
@@ -1822,7 +1822,7 @@ end
         kfsnondiff = convert(typeof(y), d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)])
     end
     kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, V, d.phi; kfs=kfsnondiff)
-    return ns, cs, d.T, P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+    return ns, cs, d.T, P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
 end
 
 @inline function calcthermo(d::ConstantTLiqFilmDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2,W<:IdealDiluteSolution,Y<:Integer,J<:AbstractArray,Q} #autodiff p
@@ -1840,7 +1840,7 @@ end
         kfsnondiff = d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
     end
     kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, V, d.phi; kfs=kfsnondiff)
-    return ns, cs, d.T, P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+    return ns, cs, d.T, P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
 end
 
 export calcthermo

From 25a2cc03eaf358cf92ccbf088331937a2a52045e Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Thu, 5 Jun 2025 12:45:16 -0700
Subject: [PATCH 12/28] adapt calcthermo dispatches for ConstantTAPhiDomain to
 handle coverage dependence

---
 src/Domain.jl | 107 ++++++++++++++++++++++++++++++++++----------------
 1 file changed, 73 insertions(+), 34 deletions(-)

diff --git a/src/Domain.jl b/src/Domain.jl
index f802fbc80..95501483b 100644
--- a/src/Domain.jl
+++ b/src/Domain.jl
@@ -1609,7 +1609,19 @@ end
     cs = ns ./ d.A
     C = N / d.A
     P = 0.0
-    return ns, cs, d.T, P, d.A, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+    if d.iscovdep
+        coverages = cs ./ d.phase.sitedensity
+        cpdivR, hdivRT, sdivR = calcHSCpdless(d.phase.vecthermo, d.T; coverages=coverages)
+        @fastmath Gs = (hdivRT .- sdivR) * (R * d.T)
+        kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, Array{Float64,1}(), d.A, d.phi; coverages=coverages)
+        return ns, cs, d.T, P, d.A, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, coverages
+    else 
+        kfs = d.kfs 
+        krevs = d.krevs
+        Gs = d.Gs
+        return ns, cs, d.T, P, d.A, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, , similar(cs,0,0)
+    end
+    
 end
 
 @inline function calcthermo(d::ConstantTAPhiDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2<:Array{Float64,1},W<:IdealSurface,Y<:Integer,J<:Array{Float64,1},Q}
@@ -1618,26 +1630,35 @@ end
     cs = ns ./ d.A
     C = N / d.A
     P = 0.0
-    if !d.alternativepformat
-        @views nothermochg = d.Gs == p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
-        if nothermochg
-            return ns, cs, d.T, P, d.A, C, N, d.mu, p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)], d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+    if !d.iscovdep
+        if !d.alternativepformat
+            @views nothermochg = d.Gs == p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
+            if nothermochg
+                return ns, cs, d.T, P, d.A, C, N, d.mu, p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)], d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, , similar(cs,0,0)
+            else
+                d.kfs = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
+                d.Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
+                d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, d.V, d.phi; kfs=d.kfs)[2]
+                return ns, cs, d.T, P, d.A, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, similar(cs,0,0)
+            end
         else
-            d.kfs = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
-            d.Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
-            d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, d.V, d.phi; kfs=d.kfs)[2]
-            return ns, cs, d.T, P, d.A, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+            @views nothermochg = d.Gs == d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
+            if nothermochg
+                return ns, cs, d.T, P, d.A, C, N, d.mu, d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)], d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, 0.0, Array{Float64,1}(), d.phi, similar(cs,0,0)
+            else
+                d.kfs = d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
+                d.Gs = d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
+                d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, d.V, d.phi; kfs=d.kfs)[2]
+                return ns, cs, d.T, P, d.A, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, similar(cs,0,0)
+            end
         end
     else
-        @views nothermochg = d.Gs == d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
-        if nothermochg
-            return ns, cs, d.T, P, d.A, C, N, d.mu, d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)], d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, 0.0, Array{Float64,1}(), d.phi
-        else
-            d.kfs = d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
-            d.Gs = d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
-            d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, d.V, d.phi; kfs=d.kfs)[2]
-            return ns, cs, d.T, P, d.A, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
-        end
+        coverages = cs ./ d.phase.sitedensity
+        cpdivR, hdivRT, sdivR = calcHSCpdless(d.phase.vecthermo, d.T; coverages=coverages)
+        @fastmath @views hdivRT .+= p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)] ./ (R * d.T)
+        @fastmath Gs = (hdivRT .- sdivR) * (R *d.T)
+        kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, Array{Float64,1}(), d.A, d.phi; coverages=coverages)
+        return ns, cs, d.T, P, d.A, C, N, d.mu, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, coverages
     end
 end
 
@@ -1647,15 +1668,24 @@ end
     cs = ns ./ d.A
     C = N / d.A
     P = 0.0
-    if !d.alternativepformat
-        Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
-        kfs = convert(typeof(y), p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)])
-    else
-        Gs = d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
-        kfs = convert(typeof(y), d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)])
+    if !d.iscovdep
+        if !d.alternativepformat
+            Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
+            kfs = convert(typeof(y), p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)])
+        else
+            Gs = d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
+            kfs = convert(typeof(y), d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)])
+        end
+        krevs = convert(typeof(y), getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, d.A, d.phi; kfs=kfs)[2])
+        return ns, cs, d.T, P, d.A, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, similar(cs,0,0)
+    else 
+        coverages = cs ./ d.phase.sitedensity
+        cpdivR, hdivRT, sdivR = calcHSCpdless(d.phase.vecthermo, d.T; coverages=coverages)
+        @fastmath @views hdivRT .+= p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)] ./ (R * d.T)
+        @fastmath Gs = (hdivRT .- sdivR) * (R * d.T)
+        kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, Array{Float64,1}(), d.A, d.phi; coverages=coverages)
+        return ns, cs, d.T, P, d.A, C, N, d.mu, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, coverages
     end
-    krevs = convert(typeof(y), getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, d.A, d.phi; kfs=kfs)[2])
-    return ns, cs, d.T, P, d.A, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
 end
 
 @inline function calcthermo(d::ConstantTAPhiDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2,W<:IdealSurface,Y<:Integer,J<:AbstractArray,Q} #autodiff p
@@ -1664,15 +1694,24 @@ end
     cs = ns ./ d.A
     C = N / d.A
     P = 0.0
-    if !d.alternativepformat
-        Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
-        kfs = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
-    else
-        Gs = d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
-        kfs = d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
+    if !d.iscovdep
+        if !d.alternativepformat
+            Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
+            kfs = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
+        else
+            Gs = d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
+            kfs = d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
+        end
+        krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, d.A, d.phi; kfs=kfs)[2]
+        return ns, cs, d.T, P, d.A, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, similar(cs,0,0)
+    else 
+        coverages = cs ./ d.phase.sitedensity
+        cpdivR, hdivRT, sdivR = calcHSCpdless(d.phase.vecthermo, d.T; coverages=coverages)
+        @fastmath @views hdivRT .+= p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)] ./ (R * T)
+        @fastmath Gs = (hdivRT .- sdivR) * (R * d.T)
+        kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, Array{Float64,1}(), d.A, d.phi; coverages=coverages)
+        return ns, cs, d.T, P, d.A, C, N, d.mu, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, coverages
     end
-    krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, d.A, d.phi; kfs=kfs)[2]
-    return ns, cs, d.T, P, d.A, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
 end
 
 function calcthermo(d::FragmentBasedConstantTrhoDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2<:SciMLBase.NullParameters,W<:FragmentBasedIdealFilm,Y<:Integer,J<:AbstractArray,Q}

From 19cfdafd4671aa3c4b79b078cb76dce0655b331a Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Thu, 5 Jun 2025 12:45:47 -0700
Subject: [PATCH 13/28] tweaking of calcthermo outputs (squash)

---
 src/Domain.jl | 108 +++++++++++++++++++++++++-------------------------
 1 file changed, 54 insertions(+), 54 deletions(-)

diff --git a/src/Domain.jl b/src/Domain.jl
index 95501483b..4b63880fe 100644
--- a/src/Domain.jl
+++ b/src/Domain.jl
@@ -931,7 +931,7 @@ export ConstantTLiqFilmDomain
     for ind in d.efficiencyinds #efficiency related rates may have changed
         d.kfs[ind], d.krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, d.Gs, d.diffusivity, V, 0.0)
     end
-    return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
+    return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealGas,Y<:Integer,J<:Array{Float64,1},Q<:Float64} #uses parameter input
@@ -952,13 +952,13 @@ end
         for ind in d.efficiencyinds #efficiency related rates may have changed
             d.kfs[ind], d.krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; f=kfps[ind])
         end
-        return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
+        return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
     elseif nothermochg
         d.kfs = kfps
         for ind in d.efficiencyinds #efficiency related rates may have changed
             d.kfs[ind], d.krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; f=kfps[ind])
         end
-        return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
+        return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
     else #need to handle thermo changes
         d.kfs .= kfps
         if !d.alternativepformat
@@ -970,7 +970,7 @@ end
         for ind in d.efficiencyinds #efficiency related rates may have changed
             d.kfs[ind], d.krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; f=kfps[ind])
         end
-        return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
+        return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
     end
 end
 
@@ -991,7 +991,7 @@ end
     for ind in d.efficiencyinds #efficiency related rates may have changed
         kfs[ind], krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0; f=kfs[ind])
     end
-    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
+    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J,Q} #Autodiff p
@@ -1011,7 +1011,7 @@ end
     for ind in d.efficiencyinds #efficiency related rates may have changed
         kfs[ind], krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0; f=kfs[ind])
     end
-    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
+    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J<:Union{ReverseDiff.TrackedArray,Tracker.TrackedArray},Q} #Autodiff p
@@ -1033,7 +1033,7 @@ end
     for ind in d.efficiencyinds #efficiency related rates may have changed
         kfs[ind], krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0; f=kfs[ind])
     end
-    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
+    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Union{ReverseDiff.TrackedArray,Tracker.TrackedArray},W<:IdealGas,Y<:Integer,J,Q} #Tracker/reversediff
@@ -1050,7 +1050,7 @@ end
         kfs = [ind in d.efficiencyinds ? getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0)[1] * d.p[length(d.phase.species)+ind] * p[d.parameterindexes[1]-1+length(d.phase.species)+ind] : p[d.parameterindexes[1]-1+length(d.phase.species)+ind] for ind in 1:length(d.phase.reactions)]
     end
     krevs = getkfkrevs(d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0; kfs=kfs)[2]
-    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
+    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Union{ReverseDiff.TrackedArray,Tracker.TrackedArray},W<:IdealGas,Y<:Integer,J<:Union{ReverseDiff.TrackedArray,Tracker.TrackedArray},Q} #Tracker/reversediff
@@ -1067,7 +1067,7 @@ end
         kfs = [ind in d.efficiencyinds ? getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0)[1] * d.p[length(d.phase.species)+ind] * p[d.parameterindexes[1]-1+length(d.phase.species)+ind] : p[d.parameterindexes[1]-1+length(d.phase.species)+ind] for ind in 1:length(d.phase.reactions)]
     end
     krevs = getkfkrevs(d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0; kfs=kfs)[2]
-    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
+    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:SciMLBase.NullParameters,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1092,7 +1092,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, d.V, 0.0)
-    return ns, cs, T, P, d.V, C, N, 0.0, kfs, krevs, Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0, nothing
+    return ns, cs, T, P, d.V, C, N, 0.0, kfs, krevs, Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1118,7 +1118,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, d.V, 0.0)
-    return @views @fastmath ns, cs, T, P, d.V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0, nothing
+    return @views @fastmath ns, cs, T, P, d.V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1144,7 +1144,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, d.V, 0.0)
-    return @views @fastmath ns, cs, T, P, d.V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0, nothing
+    return @views @fastmath ns, cs, T, P, d.V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:SciMLBase.NullParameters,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1168,7 +1168,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, d.P, C, N, ns, Gs, diffs, V, 0.0)
-    return ns, cs, T, d.P, V, C, N, 0.0, kfs, krevs, Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, nothing
+    return ns, cs, T, d.P, V, C, N, 0.0, kfs, krevs, Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1194,9 +1194,9 @@ end
     end
     kfs, krevs = getkfkrevs(d.phase, T, d.P, C, N, ns, Gs, diffs, V, 0.0)
     if p != SciMLBase.NullParameters()
-        return @views @fastmath ns, cs, T, d.P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, nothing
+        return @views @fastmath ns, cs, T, d.P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, Array{Float64,1}()
     else
-        return ns, cs, T, d.P, V, C, N, 0.0, kfs, krevs, Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, nothing
+        return ns, cs, T, d.P, V, C, N, 0.0, kfs, krevs, Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, Array{Float64,1}()
     end
 end
 
@@ -1222,7 +1222,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, d.P, C, N, ns, Gs, diffs, V, 0.0)
-    return @views @fastmath ns, cs, T, d.P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, nothing
+    return @views @fastmath ns, cs, T, d.P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ParametrizedVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:SciMLBase.NullParameters,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1247,7 +1247,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
-    return ns, cs, T, P, V, C, N, 0.0, kfs, krevs, Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0, nothing
+    return ns, cs, T, P, V, C, N, 0.0, kfs, krevs, Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ParametrizedVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1273,7 +1273,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
-    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0, nothing
+    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ParametrizedVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1299,7 +1299,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
-    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0, nothing
+    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ParametrizedPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:SciMLBase.NullParameters,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1324,7 +1324,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
-    return ns, cs, T, P, V, C, N, 0.0, kfs, krevs, Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, nothing
+    return ns, cs, T, P, V, C, N, 0.0, kfs, krevs, Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ParametrizedPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1350,7 +1350,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
-    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, nothing
+    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, Array{Float64,1}()
 end
 @inline function calcthermo(d::ParametrizedPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
     P = d.P(t)
@@ -1375,7 +1375,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
-    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, nothing
+    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ParametrizedTConstantVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:SciMLBase.NullParameters,W<:IdealDiluteSolution,Y<:Integer,J<:AbstractArray,Q}
@@ -1396,7 +1396,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, d.phi)
-    return ns, cs, T, P, V, C, N, mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0, nothing
+    return ns, cs, T, P, V, C, N, mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ParametrizedTConstantVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealDiluteSolution,Y<:Integer,J<:AbstractArray,Q}
@@ -1418,7 +1418,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, d.phi)
-    return @views @fastmath ns, cs, T, P, V, C, N, mu, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), d.phi, nothing
+    return @views @fastmath ns, cs, T, P, V, C, N, mu, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ParametrizedTConstantVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealDiluteSolution,Y<:Integer,J<:AbstractArray,Q}
@@ -1440,7 +1440,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, d.phi)
-    return @views @fastmath ns, cs, T, P, V, C, N, mu, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), d.phi, nothing
+    return @views @fastmath ns, cs, T, P, V, C, N, mu, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ParametrizedTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:SciMLBase.NullParameters,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1465,7 +1465,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
-    return ns, cs, T, P, V, C, N, 0.0, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0, nothing
+    return ns, cs, T, P, V, C, N, 0.0, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ParametrizedTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1491,7 +1491,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
-    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0, nothing
+    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ParametrizedTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1517,7 +1517,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
-    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0, nothing
+    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantTVDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2<:SciMLBase.NullParameters,W<:IdealDiluteSolution,Y<:Integer,J<:AbstractArray,Q}
@@ -1526,7 +1526,7 @@ end
     cs = ns ./ d.V
     C = N / d.V
     P = 1.0e8
-    return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
+    return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantTVDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2<:Array{Float64,1},W<:IdealDiluteSolution,Y<:Integer,J<:Array{Float64,1},Q}
@@ -1539,31 +1539,31 @@ end
         @views nothermochg = d.Gs == p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
         @views nokfchg = d.kfsnondiff == p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
         if nothermochg && nokfchg
-            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
+            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}()
         elseif nothermochg
             d.kfsnondiff = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, d.V, d.phi; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
+            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}()
         else
             d.kfsnondiff = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, d.V, d.phi; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
+            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}()
         end
     else
         @views nothermochg = d.Gs == d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
         @views nokfchg = d.kfsnondiff == d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
         if nothermochg && nokfchg
-            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
+            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}()
         elseif nothermochg
             d.kfsnondiff .= d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, d.V, d.phi; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
+            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}()
         else
             d.kfsnondiff .= d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.Gs .= d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, d.V, d.phi; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
+            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}()
         end
     end
 
@@ -1583,7 +1583,7 @@ end
         kfsnondiff = convert(typeof(y), d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)])
     end
     kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, d.V, d.phi; kfs=kfsnondiff)
-    return ns, cs, d.T, P, d.V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
+    return ns, cs, d.T, P, d.V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantTVDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2,W<:IdealDiluteSolution,Y<:Integer,J<:AbstractArray,Q} #autodiff p
@@ -1600,7 +1600,7 @@ end
         kfsnondiff = d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
     end
     kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, d.V, d.phi; kfs=kfsnondiff)
-    return ns, cs, d.T, P, d.V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
+    return ns, cs, d.T, P, d.V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantTAPhiDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2<:SciMLBase.NullParameters,W<:IdealSurface,Y<:Integer,J<:AbstractArray,Q}
@@ -1722,7 +1722,7 @@ function calcthermo(d::FragmentBasedConstantTrhoDomain{W,Y}, y::J, t::Q, p::Q2=S
     cs = ns ./ V
     C = N / V
     P = 1.0e8
-    return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
+    return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
 end
 
 function calcthermo(d::FragmentBasedConstantTrhoDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2<:Array{Float64,1},W<:FragmentBasedIdealFilm,Y<:Integer,J<:Array{Float64,1},Q}
@@ -1737,31 +1737,31 @@ function calcthermo(d::FragmentBasedConstantTrhoDomain{W,Y}, y::J, t::Q, p::Q2=S
         @views nothermochg = d.Gs == p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
         @views nokfchg = d.kfsnondiff == p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
         if nothermochg && nokfchg
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
         elseif nothermochg
             d.kfsnondiff = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
         else
             d.kfsnondiff = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}
         end
     else
         @views nothermochg = d.Gs == d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
         @views nokfchg = d.kfsnondiff == d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
         if nothermochg && nokfchg
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}
         elseif nothermochg
             d.kfsnondiff .= d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}
         else
             d.kfsnondiff .= d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.Gs .= d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}
         end
     end
 
@@ -1783,7 +1783,7 @@ function calcthermo(d::FragmentBasedConstantTrhoDomain{W,Y}, y::Array{W2,1}, t::
         kfsnondiff = convert(typeof(y), d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)])
     end
     kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, V, 0.0; kfs=kfsnondiff)
-    return ns, cs, d.T, P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
+    return ns, cs, d.T, P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}
 end
 
 function calcthermo(d::FragmentBasedConstantTrhoDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2,W<:FragmentBasedIdealFilm,Y<:Integer,J<:AbstractArray,Q} #autodiff p
@@ -1802,7 +1802,7 @@ function calcthermo(d::FragmentBasedConstantTrhoDomain{W,Y}, y::J, t::Q, p::Q2=S
         kfsnondiff = d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
     end
     kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, V, 0.0; kfs=kfsnondiff)
-    return ns, cs, d.T, P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, nothing
+    return ns, cs, d.T, P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}
 end
 
 @inline function calcthermo(d::ConstantTLiqFilmDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2<:Array{Float64,1},W<:IdealDiluteSolution,Y<:Integer,J<:Array{Float64,1},Q}
@@ -1816,31 +1816,31 @@ end
         @views nothermochg = d.Gs == p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
         @views nokfchg = d.kfsnondiff == p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
         if nothermochg && nokfchg
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}
         elseif nothermochg
             d.kfsnondiff = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, d.phi; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}
         else
             d.kfsnondiff = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, d.phi; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}
         end
     else
         @views nothermochg = d.Gs == d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
         @views nokfchg = d.kfsnondiff == d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
         if nothermochg && nokfchg
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}
         elseif nothermochg
             d.kfsnondiff .= d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, d.phi; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}
         else
             d.kfsnondiff .= d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.Gs .= d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, d.phi; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}
         end
     end
 
@@ -1861,7 +1861,7 @@ end
         kfsnondiff = convert(typeof(y), d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)])
     end
     kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, V, d.phi; kfs=kfsnondiff)
-    return ns, cs, d.T, P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
+    return ns, cs, d.T, P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}
 end
 
 @inline function calcthermo(d::ConstantTLiqFilmDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2,W<:IdealDiluteSolution,Y<:Integer,J<:AbstractArray,Q} #autodiff p
@@ -1879,7 +1879,7 @@ end
         kfsnondiff = d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
     end
     kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, V, d.phi; kfs=kfsnondiff)
-    return ns, cs, d.T, P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, nothing
+    return ns, cs, d.T, P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}
 end
 
 export calcthermo

From c3c4d671bf537122891790437b450bda1012ca6f Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Thu, 5 Jun 2025 12:46:24 -0700
Subject: [PATCH 14/28] adapting ReactiveInternalInterface to coverage
 dependence

---
 src/Interface.jl | 20 ++++++++++++++------
 1 file changed, 14 insertions(+), 6 deletions(-)

diff --git a/src/Interface.jl b/src/Interface.jl
index 466022abc..456261e00 100644
--- a/src/Interface.jl
+++ b/src/Interface.jl
@@ -49,20 +49,28 @@ function ReactiveInternalInterface(domain1, domain2, reactions, A)
 end
 export ReactiveInternalInterface
 
-function getkfskrevs(ri::ReactiveInternalInterface, T1, T2, phi1, phi2, Gs1, Gs2, cstot::Array{Q,1}) where {Q}
-    kfs = getkfs(ri, T1, 0.0, 0.0, Array{Q,1}(), ri.A, phi1)
+function getkfskrevs(ri::ReactiveInternalInterface, T1, T2, phi1, phi2, Gs1, Gs2, coverages1, coverages2, cstot::Array{Q,1}) where {Q}
+    if !(coverages1 === nothing) && coverages2 === nothing
+        kfs = getkfs(ri, T1, 0.0, 0.0, Array{Q,1}(), ri.A, phi1; coverages=coverages1)
+    elseif !(coverages2 === nothing) && coverages1 === nothing
+        kfs = getkfs(ri, T1, 0.0, 0.0, Array{Q,1}(), ri.A, phi1; coverages=coverages2)
+    elseif coverages1 === nothing && coverages2 === nothing 
+        kfs = getkfs(ri, T1, 0.0, 0.0, Array{Q,1}(), ri.A, phi1)
+    else 
+        throw(DomainError("Coverage dependence of an interface between two surfaces is indeterminate (which coverage should be used?)"))
+    end
     Kc = getKc.(ri.reactions, ri.domain1.phase, ri.domain2.phase, Ref(Gs1), Ref(Gs2), T1, phi1)
     krevs = kfs ./ Kc
     return kfs, krevs
 end
 
-function evaluate(ri::ReactiveInternalInterface, dydt, domains, T1, T2, phi1, phi2, Gs1, Gs2, cstot, p::W) where {W<:SciMLBase.NullParameters}
-    kfs, krevs = getkfskrevs(ri, T1, T2, phi1, phi2, Gs1, Gs2, cstot)
+function evaluate(ri::ReactiveInternalInterface, dydt, domains, T1, T2, phi1, phi2, Gs1, Gs2, coverages1, coverages2, cstot, p::W) where {W<:SciMLBase.NullParameters}
+    kfs, krevs = getkfskrevs(ri, T1, T2, phi1, phi2, Gs1, Gs2, coverages1, coverages2, cstot)
     addreactionratecontributions!(dydt, ri.rxnarray, cstot, kfs, krevs, ri.A)
 end
 
-function evaluate(ri::ReactiveInternalInterface, dydt, domains, T1, T2, phi1, phi2, Gs1, Gs2, cstot, p)
-    kfs, krevs = getkfskrevs(ri, T1, T2, phi1, phi2, Gs1, Gs2, cstot)
+function evaluate(ri::ReactiveInternalInterface, dydt, domains, T1, T2, phi1, phi2, Gs1, Gs2, coverages1, coverages2, cstot, p)
+    kfs, krevs = getkfskrevs(ri, T1, T2, phi1, phi2, Gs1, Gs2, coverages1, coverages2, cstot)
     addreactionratecontributions!(dydt, ri.rxnarray, cstot, kfs .* p[ri.parameterindexes[1]:ri.parameterindexes[2]], krevs .* p[ri.parameterindexes[1]:ri.parameterindexes[2]], ri.A)
 end
 export evaluate

From 7ee4048aae6357796f47a60d172b327156a32da7 Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Thu, 5 Jun 2025 12:46:47 -0700
Subject: [PATCH 15/28] adapting ReactiveInternalInterfaceConstantTPhi to
 coverage dependence

---
 src/Interface.jl | 95 ++++++++++++++++++++++++++++++++++++++----------
 1 file changed, 76 insertions(+), 19 deletions(-)

diff --git a/src/Interface.jl b/src/Interface.jl
index 456261e00..edb8d781c 100644
--- a/src/Interface.jl
+++ b/src/Interface.jl
@@ -93,12 +93,15 @@ struct ReactiveInternalInterfaceConstantTPhi{J,N,B,B2,B3,C,C2,Q<:AbstractReactio
     p::Array{Float64,1}
     reversibility::Array{Bool,1}
     forwardability::Array{Bool,1}
+    iscovdep::Bool
+    phi::Float64
 end
 function ReactiveInternalInterfaceConstantTPhi(domain1, domain2, reactions, T, A, phi=0.0)
     @assert domain1.T == domain2.T
     reactions = upgradekinetics(reactions, domain1, domain2)
     rxnarray = getinterfacereactioninds(domain1, domain2, reactions)
-    kfs = getkf.(reactions, nothing, T, 0.0, 0.0, Ref([]), A, phi)
+    iscovdep = any([!(spc.thermo.covdep isa EmptyThermoCoverageDependence) for spc in domain1.phase.species]) || any([!(spc.thermo.covdep isa EmptyThermoCoverageDependence) for spc in domain2.phase.species]) || any([!(rxn.kinetics.covdep isa EmptyRateCoverageDependence) for rxn in reactions])
+    kfs = getkf.(reactions, nothing, T, 0.0, 0.0, Ref([]), A, phi; coverages=nothing)
     Kc = getKc.(reactions, domain1.phase, domain2.phase, Ref(domain1.Gs), Ref(domain2.Gs), T, phi)
     krevs = kfs ./ Kc
     M, Nrp1, Nrp2 = getstoichmatrix(domain1, domain2, reactions)
@@ -110,33 +113,87 @@ function ReactiveInternalInterfaceConstantTPhi(domain1, domain2, reactions, T, A
     if isa(kfs, Vector{Any})
         kfs = convert(Vector{Float64}, kfs)
     end
-    return ReactiveInternalInterfaceConstantTPhi(domain1, domain2, reactions,
-        rxnarray, M, Nrp1, Nrp2, kfs, krevs, T, A, [1, length(reactions)],
-        [0, 1], kfs[1:end], reversibility, forwardability), kfs[1:end]
+    if !iscovdep
+        return ReactiveInternalInterfaceConstantTPhi(domain1, domain2, reactions,
+            rxnarray, M, Nrp1, Nrp2, kfs, krevs, T, A, [1, length(reactions)],
+            [0, 1], kfs[1:end], reversibility, forwardability, iscovdep, phi), kfs[1:end]
+    else 
+        return ReactiveInternalInterfaceConstantTPhi(domain1, domain2, reactions,
+            rxnarray, M, Nrp1, Nrp2, kfs, krevs, T, A, [1, length(reactions)],
+            [0, 1], kfs[1:end], reversibility, forwardability, iscovdep, phi), ones(length(kfs))
+    end
 end
 export ReactiveInternalInterfaceConstantTPhi
 
-function getkfskrevs(ri::ReactiveInternalInterfaceConstantTPhi, T1, T2, phi1, phi2, Gs1, Gs2, cstot)
-    return ri.kfs, ri.krevs
-end
-
-function evaluate(ri::ReactiveInternalInterfaceConstantTPhi, dydt, domains, T1, T2, phi1, phi2, Gs1, Gs2, cstot, p::W) where {W<:SciMLBase.NullParameters}
-    addreactionratecontributions!(dydt, ri.rxnarray, cstot, ri.kfs, ri.krevs, ri.A)
+function getkfskrevs(ri::ReactiveInternalInterfaceConstantTPhi, T1, T2, phi1, phi2, Gs1, Gs2, coverages1, coverages2, cstot)
+    if ri.iscovdep
+        if !(coverages1 === nothing) && coverages2 === nothing
+          coverages = coverages1
+        elseif !(coverages2 === nothing) && coverages1 === nothing
+            coverages = coverages2
+        elseif coverages1 === nothing && coverages2 === nothing 
+            coverages = nothing
+        else 
+            throw(DomainError("Coverage dependence of an interface between two surfaces is indeterminate (which coverage should be used?)"))
+        end
+        kfs = getkf.(ri.reactions, nothing, ri.T, 0.0, 0.0, Ref([]), ri.A, ri.phi; coverages=coverages)
+        Kc = getKc.(ri.reactions, ri.domain1.phase, ri.domain2.phase, Ref(Gs1), Ref(Gs2), ri.T, ri.phi)
+        krevs = kfs ./ Kc
+        return kfs,krevs
+    else
+        return ri.kfs, ri.krevs
+    end
 end
 
-function evaluate(ri::ReactiveInternalInterfaceConstantTPhi, dydt, domains, T1, T2, phi1, phi2, Gs1, Gs2, cstot, p)
-    if p[ri.parameterindexes[1]:ri.parameterindexes[2]] == ri.kfs
-        kfs = ri.kfs
+function evaluate(ri::ReactiveInternalInterfaceConstantTPhi, dydt, domains, T1, T2, phi1, phi2, Gs1, Gs2, coverages1, coverages2, cstot, p::W) where {W<:SciMLBase.NullParameters}
+    if ri.iscovdep
+        if !(coverages1 === nothing) && coverages2 === nothing
+          coverages = coverages1
+        elseif !(coverages2 === nothing) && coverages1 === nothing
+            coverages = coverages2
+        elseif coverages1 === nothing && coverages2 === nothing 
+            coverages = nothing
+        else 
+            throw(DomainError("Coverage dependence of an interface between two surfaces is indeterminate (which coverage should be used?)"))
+        end
+        kfs = getkf.(ri.reactions, nothing, ri.T, 0.0, 0.0, Ref([]), ri.A, ri.phi; coverages=coverages)
+        Kc = getKc.(ri.reactions, ri.domain1.phase, ri.domain2.phase, Ref(Gs1), Ref(Gs2), ri.T, phi1)
+        krevs = kfs ./ Kc
+        addreactionratecontributions!(dydt, ri.rxnarray, cstot, kfs, krevs, ri.A)
     else
-        kfs = p[ri.parameterindexes[1]:ri.parameterindexes[2]]
+        addreactionratecontributions!(dydt, ri.rxnarray, cstot, ri.kfs, ri.krevs, ri.A)
     end
-    if length(Gs1) == 0 || length(Gs2) == 0 || (all(Gs1 .== ri.domain1.Gs) && all(Gs2 .== ri.domain2.Gs))
-        krevs = ri.krevs
-    else
-        Kc = getKcs(ri, T1, Gs1, Gs2)
+end
+
+function evaluate(ri::ReactiveInternalInterfaceConstantTPhi, dydt, domains, T1, T2, phi1, phi2, Gs1, Gs2, coverages1, coverages2, cstot, p)
+    if ri.iscovdep
+        if !(coverages1 === nothing) && coverages2 === nothing
+            coverages = coverages1
+        elseif !(coverages2 === nothing) && coverages1 === nothing
+            coverages = coverages2
+        elseif coverages1 === nothing && coverages2 === nothing 
+            coverages = nothing
+        else 
+            throw(DomainError("Coverage dependence of an interface between two surfaces is indeterminate (which coverage should be used?)"))
+        end
+        kfs = getkf.(ri.reactions, nothing, ri.T, 0.0, 0.0, Ref([]), ri.A, ri.phi; coverages=coverages) .* p[ri.parameterindexes[1]:ri.parameterindexes[2]] 
+        Kc = getKc.(ri.reactions, ri.domain1.phase, ri.domain2.phase, Ref(Gs1), Ref(Gs2), ri.T, ri.phi)
         krevs = kfs ./ Kc
+        addreactionratecontributions!(dydt, ri.rxnarray, cstot, kfs, krevs, ri.A)
+    else
+        if p[ri.parameterindexes[1]:ri.parameterindexes[2]] == ri.kfs
+            kfs = ri.kfs
+        else
+            kfs = p[ri.parameterindexes[1]:ri.parameterindexes[2]]
+        end
+        if length(Gs1) == 0 || length(Gs2) == 0 || (all(Gs1 .== ri.domain1.Gs) && all(Gs2 .== ri.domain2.Gs))
+            krevs = ri.krevs
+        else
+            Kc = getKcs(ri, T1, Gs1, Gs2)
+            krevs = kfs ./ Kc
+        end
+        addreactionratecontributions!(dydt, ri.rxnarray, cstot, kfs, krevs, ri.A)
     end
-    addreactionratecontributions!(dydt, ri.rxnarray, cstot, kfs, krevs, ri.A)
 end
 export evaluate
 

From f49167df8327bf3f581e53bdaeb4ca2f610bb90b Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Thu, 5 Jun 2025 12:47:26 -0700
Subject: [PATCH 16/28] adapting multi-domain reactor to handle coverage
 dependence

---
 src/Reactor.jl | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/Reactor.jl b/src/Reactor.jl
index 549e9ead4..53b79aaf8 100644
--- a/src/Reactor.jl
+++ b/src/Reactor.jl
@@ -647,7 +647,7 @@ end
     cstot .= 0.0
     dydt .= 0.0
     domain = domains[1]
-    ns, cs, T, P, V, C, N, mu, kfs, krevs, Hs, Us, Gs, diffs, Cvave, cpdivR, phi = calcthermo(domain, y, t, p)
+    ns, cs, T, P, V, C, N, mu, kfs, krevs, Hs, Us, Gs, diffs, Cvave, cpdivR, phi, coverages = calcthermo(domain, y, t, p)
     vns = Array{Any,1}(undef, length(domains))
     vns[1] = ns
     vcs = Array{Any,1}(undef, length(domains))
@@ -683,11 +683,13 @@ end
     vcpdivR[1] = cpdivR
     vphi = Array{Any,1}(undef, length(domains))
     vphi[1] = phi
+    vcoverages = Array{Any,1}(undef, length(domains))
+    vcoverages[1] = coverages
     addreactionratecontributions!(dydt, domain.rxnarray, cstot, kfs, krevs)
     @views dydt[domain.indexes[1]:domain.indexes[2]] .*= V
     for (i, domain) in enumerate(@views domains[2:end])
         k = i + 1
-        vns[k], vcs[k], vT[k], vP[k], vV[k], vC[k], vN[k], vmu[k], vkfs[k], vkrevs[k], vHs[k], vUs[k], vGs[k], vdiffs[k], vCvave[k], vcpdivR[k], vphi[k] = calcthermo(domain, y, t, p)
+        vns[k], vcs[k], vT[k], vP[k], vV[k], vC[k], vN[k], vmu[k], vkfs[k], vkrevs[k], vHs[k], vUs[k], vGs[k], vdiffs[k], vCvave[k], vcpdivR[k], vphi[k], vcoverages[k] = calcthermo(domain, y, t, p)
         cstot[domain.indexes[1]:domain.indexes[2]] .= vcs[k]
         addreactionratecontributions!(dydt, domain.rxnarray, cstot, vkfs[k], vkrevs[k])
         @views dydt[domain.indexes[1]:domain.indexes[2]] .*= vV[k]
@@ -696,7 +698,7 @@ end
         if isa(inter, FragmentBasedReactiveFilmGrowthInterfaceConstantT)
             evaluate(inter, dydt, vV[inter.domaininds[1]], cstot)
         elseif isa(inter, AbstractReactiveInternalInterface)
-            evaluate(inter, dydt, domains, vT[inter.domaininds[1]], vT[inter.domaininds[2]], vphi[inter.domaininds[1]], vphi[inter.domaininds[2]], vGs[inter.domaininds[1]], vGs[inter.domaininds[2]], cstot, p)
+            evaluate(inter, dydt, domains, vT[inter.domaininds[1]], vT[inter.domaininds[2]], vphi[inter.domaininds[1]], vphi[inter.domaininds[2]], vGs[inter.domaininds[1]], vGs[inter.domaininds[2]], vcoverages[inter.domaininds[1]], vcoverages[inter.domaininds[2]], cstot, p)
         elseif isa(inter, VaporLiquidMassTransferInternalInterfaceConstantT)
             evaluate(inter, dydt, vV[inter.domaininds[1]], vV[inter.domaininds[2]], vT[inter.domaininds[1]], vT[inter.domaininds[2]], vN[inter.domaininds[1]], vN[inter.domaininds[2]], vP[inter.domaininds[1]], vP[inter.domaininds[2]], vCvave[inter.domaininds[1]], vCvave[inter.domaininds[2]], vns[inter.domaininds[1]], vns[inter.domaininds[2]], vUs[inter.domaininds[1]], vUs[inter.domaininds[2]], cstot, p)
         end

From b5c82785931c1482cec577a7ab84b3320ce6384c Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Thu, 5 Jun 2025 12:50:27 -0700
Subject: [PATCH 17/28] adapt simulation reactor evaluations to handle coverage
 dependence

also fixed bug related to handling of cpdivR in output from calcthermo
---
 src/Simulation.jl | 26 +++++++++++++++++---------
 1 file changed, 17 insertions(+), 9 deletions(-)

diff --git a/src/Simulation.jl b/src/Simulation.jl
index f80f614d4..653ad5911 100644
--- a/src/Simulation.jl
+++ b/src/Simulation.jl
@@ -53,7 +53,7 @@ function Simulation(sol::Q, domain::W, reducedmodelmappings::ReducedModelMapping
         end
         @inbounds @views yunlumped[end-length(domain.thermovariabledict)+1:end] .= y[end-length(domain.thermovariabledict)+1:end]
 
-        ns, cs, T, P, V, C, N, mu, kfs, krevs, Hs, Us, Gs, diffs, Cvave, cpdivR, phi = calcthermo(domain, yunlumped, t, p)
+        ns, cs, T, P, V, C, N, mu, kfs, krevs, Hs, Us, Gs, diffs, Cvave, cpdivR, phi, coverages = calcthermo(domain, yunlumped, t, p)
 
         reducedmodelmappings.qssc!(qssc, cs, kfs, krevs)
         @inbounds yunlumped[reducedmodelmappings.qssindexes] .= qssc .* V
@@ -272,10 +272,12 @@ function rops(ssys::SystemSimulation, t)
     vdiffs = Array{Any,1}(undef, length(domains))
     vCvave = Array{Any,1}(undef, length(domains))
     vphi = Array{Any,1}(undef, length(domains))
+    vcpdivR = Array{Any,1}(undef, length(domains))
+    vcoverages = Array{Any,1}(undef, length(domains))
     ropmat = spzeros(Nrxns, Nspcs)
     start = 0
     for (k, sim) in enumerate(ssys.sims)
-        vns[k], vcs[k], vT[k], vP[k], vV[k], vC[k], vN[k], vmu[k], vkfs[k], vkrevs[k], vHs[k], vUs[k], vGs[k], vdiffs[k], vCvave[k], vphi[k] = calcthermo(sim.domain, ssys.sol(t), t)
+        vns[k], vcs[k], vT[k], vP[k], vV[k], vC[k], vN[k], vmu[k], vkfs[k], vkrevs[k], vHs[k], vUs[k], vGs[k], vdiffs[k], vCvave[k], vcpdivR[k], vphi[k], vcoverages[k] = calcthermo(sim.domain, ssys.sol(t), t)
         cstot[sim.domain.indexes[1]:sim.domain.indexes[2]] = vcs[k]
         rops!(ropmat, sim.domain.rxnarray, cstot, vkfs[k], vkrevs[k], vV[k], start)
         start += length(vkfs[k])
@@ -286,7 +288,7 @@ function rops(ssys::SystemSimulation, t)
             rops!(ropmat, nothing, inter.rxnarray, inter.fragmentbasedrxnarray, cstot, kfs, krevs, vV[inter.domaininds[1]], inter.Mws, inter.domainfilm.indexes[1]:inter.domainfilm.indexes[2], start)
             start += length(kfs)
         elseif hasproperty(inter, :reactions)
-            kfs, krevs = getkfskrevs(inter, vT[inter.domaininds[1]], vT[inter.domaininds[2]], vphi[inter.domaininds[1]], vphi[inter.domaininds[2]], vGs[inter.domaininds[1]], vGs[inter.domaininds[2]], cstot)
+            kfs, krevs = getkfskrevs(inter, vT[inter.domaininds[1]], vT[inter.domaininds[2]], vphi[inter.domaininds[1]], vphi[inter.domaininds[2]], vGs[inter.domaininds[1]], vGs[inter.domaininds[2]], vcoverages[inter.domaininds[1]], vcoverages[inter.domaininds[2]], cstot)
             rops!(ropmat, inter.rxnarray, cstot, kfs, krevs, inter.A, start)
             start += length(kfs)
         end
@@ -316,7 +318,9 @@ function massrops(ssys::SystemSimulation, t)
     vGs = Array{Any,1}(undef, length(domains))
     vdiffs = Array{Any,1}(undef, length(domains))
     vCvave = Array{Any,1}(undef, length(domains))
+    vcpdivR = Array{Any,1}(undef, length(domains))
     vphi = Array{Any,1}(undef, length(domains))
+    vcoverages = Array{Any,1}(undef, length(domains))
     if any([domain isa FragmentBasedConstantTrhoDomain for domain in domains])
         ropmat = spzeros(Nrxns, Nspcs + 1)
     else
@@ -326,7 +330,7 @@ function massrops(ssys::SystemSimulation, t)
     ropmat = spzeros(Nrxns, Nspcs)
     start = 0
     for (k, sim) in enumerate(ssys.sims)
-        vns[k], vcs[k], vT[k], vP[k], vV[k], vC[k], vN[k], vmu[k], vkfs[k], vkrevs[k], vHs[k], vUs[k], vGs[k], vdiffs[k], vCvave[k], vphi[k] = calcthermo(sim.domain, ssys.sol(t), t)
+        vns[k], vcs[k], vT[k], vP[k], vV[k], vC[k], vN[k], vmu[k], vkfs[k], vkrevs[k], vHs[k], vUs[k], vGs[k], vdiffs[k], vCvave[k], vcpdivR[k], vphi[k], vcoverages[k] = calcthermo(sim.domain, ssys.sol(t), t)
         cstot[sim.domain.indexes[1]:sim.domain.indexes[2]] = vcs[k]
         rops!(ropmat, sim.domain.rxnarray, cstot, vkfs[k], vkrevs[k], vV[k], start)
         start += length(vkfs[k])
@@ -337,7 +341,7 @@ function massrops(ssys::SystemSimulation, t)
             rops!(ropmat, massropvec, inter.rxnarray, inter.fragmentbasedrxnarray, cstot, kfs, krevs, vV[inter.domaininds[1]], inter.Mws, inter.domainfilm.indexes[1]:inter.domainfilm.indexes[2], start)
             start += length(kfs)
         elseif hasproperty(inter, :reactions)
-            kfs, krevs = getkfskrevs(inter, vT[inter.domaininds[1]], vT[inter.domaininds[2]], vphi[inter.domaininds[1]], vphi[inter.domaininds[2]], vGs[inter.domaininds[1]], vGs[inter.domaininds[2]], cstot)
+            kfs, krevs = getkfskrevs(inter, vT[inter.domaininds[1]], vT[inter.domaininds[2]], vphi[inter.domaininds[1]], vphi[inter.domaininds[2]], vGs[inter.domaininds[1]], vGs[inter.domaininds[2]], vcoverages[inter.domaininds[1]], vcoverages[inter.domaininds[2]], cstot)
             rops!(ropmat, inter.rxnarray, cstot, kfs, krevs, inter.A, start)
             start += length(kfs)
         end
@@ -392,17 +396,19 @@ function rops(ssys::SystemSimulation, name, t)
     vdiffs = Array{Any,1}(undef, length(domains))
     vCvave = Array{Any,1}(undef, length(domains))
     vphi = Array{Any,1}(undef, length(domains))
+    vcpdivR = Array{Any,1}(undef, length(domains))
+    vcoverages = Array{Any,1}(undef, length(domains))
     ropvec = spzeros(Nrxns)
     start = 0
     for (k, sim) in enumerate(ssys.sims)
-        vns[k], vcs[k], vT[k], vP[k], vV[k], vC[k], vN[k], vmu[k], vkfs[k], vkrevs[k], vHs[k], vUs[k], vGs[k], vdiffs[k], vCvave[k], vphi[k] = calcthermo(sim.domain, ssys.sol(t), t)
+        vns[k], vcs[k], vT[k], vP[k], vV[k], vC[k], vN[k], vmu[k], vkfs[k], vkrevs[k], vHs[k], vUs[k], vGs[k], vdiffs[k], vCvave[k], vcpdivR[k], vphi[k], vcoverages[k] = calcthermo(sim.domain, ssys.sol(t), t)
         cstot[sim.domain.indexes[1]:sim.domain.indexes[2]] = vcs[k]
         rops!(ropvec, sim.domain.rxnarray, cstot, vkfs[k], vkrevs[k], vV[k], start, ind)
         start += length(vkfs[k])
     end
     for inter in ssys.interfaces
         if hasproperty(inter, :reactions)
-            kfs, krevs = getkfskrevs(inter, vT[inter.domaininds[1]], vT[inter.domaininds[2]], vphi[inter.domaininds[1]], vphi[inter.domaininds[2]], vGs[inter.domaininds[1]], vGs[inter.domaininds[2]], cstot)
+            kfs, krevs = getkfskrevs(inter, vT[inter.domaininds[1]], vT[inter.domaininds[2]], vphi[inter.domaininds[1]], vphi[inter.domaininds[2]], vGs[inter.domaininds[1]], vGs[inter.domaininds[2]], vcoverages[inter.domaininds[1]], vcoverages[inter.domaininds[2]], cstot)
             rops!(ropvec, inter.rxnarray, cstot, kfs, krevs, inter.A, start, ind)
             start += length(kfs)
         end
@@ -800,9 +806,11 @@ function rates(ssys::Q, t::X) where {Q<:SystemSimulation,X<:Real}
     vdiffs = Array{Any,1}(undef, length(domains))
     vCvave = Array{Any,1}(undef, length(domains))
     vphi = Array{Any,1}(undef, length(domains))
+    vcpdivR = Array{Any,1}(undef, length(domains))
+    vcoverages = Array{Any,1}(undef, length(domains))
     index = 1
     for (k, sim) in enumerate(ssys.sims)
-        vns[k], vcs[k], vT[k], vP[k], vV[k], vC[k], vN[k], vmu[k], vkfs[k], vkrevs[k], vHs[k], vUs[k], vGs[k], vdiffs[k], vCvave[k], vphi[k] = calcthermo(sim.domain, ssys.sol(t), t)
+        vns[k], vcs[k], vT[k], vP[k], vV[k], vC[k], vN[k], vmu[k], vkfs[k], vkrevs[k], vHs[k], vUs[k], vGs[k], vdiffs[k], vCvave[k], vcpdivR[k], vphi[k], vcoverages[k] = calcthermo(sim.domain, ssys.sol(t), t)
         cstot[sim.domain.indexes[1]:sim.domain.indexes[2]] = vcs[k]
         rts[index:index+length(vkfs[k])-1] .= getrates(sim.domain.rxnarray, cstot, vkfs[k], vkrevs[k]) .* getdomainsize(sim, t)
         index += length(vkfs[k])
@@ -813,7 +821,7 @@ function rates(ssys::Q, t::X) where {Q<:SystemSimulation,X<:Real}
             @views rts[index:index+length(kfs)-1] = getrates(inter.rxnarray, cstot, kfs, krevs) .* vV[inter.domaininds[1]]
             index += length(kfs)
         elseif hasproperty(inter, :reactions)
-            kfs, krevs = getkfskrevs(inter, vT[inter.domaininds[1]], vT[inter.domaininds[2]], vphi[inter.domaininds[1]], vphi[inter.domaininds[2]], vGs[inter.domaininds[1]], vGs[inter.domaininds[2]], cstot)
+            kfs, krevs = getkfskrevs(inter, vT[inter.domaininds[1]], vT[inter.domaininds[2]], vphi[inter.domaininds[1]], vphi[inter.domaininds[2]], vGs[inter.domaininds[1]], vGs[inter.domaininds[2]], vcoverages[inter.domaininds[1]], vcoverages[inter.domaininds[2]], cstot)
             rts[index:index+length(kfs)-1] = getrates(inter.rxnarray, cstot, kfs, krevs) .* inter.A
             index += length(kfs)
         end

From d0d2ae8304383667b060129b37206dbc90a10f88 Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Thu, 5 Jun 2025 12:50:44 -0700
Subject: [PATCH 18/28] add coverage dependence objects to parser

---
 src/Parse.jl | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/src/Parse.jl b/src/Parse.jl
index 3698081f3..5251e06bb 100644
--- a/src/Parse.jl
+++ b/src/Parse.jl
@@ -4,7 +4,9 @@ using PyCall
 using StaticArrays
 
 module Calc
+    include("Calculators/RateCoverageDependence.jl")
     include("Calculators/RateUncertainty.jl")
+    include("Calculators/ThermoCoverageDependence.jl")
     include("Calculators/ThermoUncertainty.jl")
     include("Calculators/Thermo.jl")
     include("Calculators/Diffusion.jl")
@@ -15,6 +17,7 @@ module Calc
     include("Calculators/kLAkH.jl")
 end
 module Spc
+    include("Calculators/ThermoCoverageDependence.jl")
     include("Calculators/ThermoUncertainty.jl")
     include("Calculators/Thermo.jl")
     include("Calculators/Diffusion.jl")
@@ -22,7 +25,9 @@ module Spc
     include("Species.jl")
 end
 module Rxn
+    include("Calculators/RateCoverageDependence.jl")
     include("Calculators/RateUncertainty.jl")
+    include("Calculators/ThermoCoverageDependence.jl")
     include("Calculators/ThermoUncertainty.jl")
     include("Calculators/Thermo.jl")
     include("Calculators/Diffusion.jl")

From 32645d311928ae804c234d76d3c4a4d0c560b05e Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Thu, 5 Jun 2025 12:51:09 -0700
Subject: [PATCH 19/28] add coverage dependence objects to the main module and
 testing

---
 src/ReactionMechanismSimulator.jl | 2 ++
 src/rmstest.jl                    | 2 ++
 2 files changed, 4 insertions(+)

diff --git a/src/ReactionMechanismSimulator.jl b/src/ReactionMechanismSimulator.jl
index 6e30c01dd..40e520141 100644
--- a/src/ReactionMechanismSimulator.jl
+++ b/src/ReactionMechanismSimulator.jl
@@ -52,7 +52,9 @@ module ReactionMechanismSimulator
     end
     include("Constants.jl")
     include("Tools.jl")
+    include("Calculators/RateCoverageDependence.jl")
     include("Calculators/RateUncertainty.jl")
+    include("Calculators/ThermoCoverageDependence.jl")
     include("Calculators/ThermoUncertainty.jl")
     include("Calculators/Thermo.jl")
     include("Calculators/Diffusion.jl")
diff --git a/src/rmstest.jl b/src/rmstest.jl
index 1ffa8bcce..f1a939c4f 100644
--- a/src/rmstest.jl
+++ b/src/rmstest.jl
@@ -19,7 +19,9 @@ copy!(pydot,pyimport_conda("pydot","pydot","rmg"))
 
 include("Constants.jl")
 include("Tools.jl")
+include("Calculators/RateCoverageDependence.jl")
 include("Calculators/RateUncertainty.jl")
+include("Calculators/ThermoCoverageDependence.jl")
 include("Calculators/ThermoUncertainty.jl")
 include("Calculators/Thermo.jl")
 include("Calculators/Diffusion.jl")

From 4b114d227f574c6009ac1c9a1a486b6604c35571 Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Thu, 5 Jun 2025 12:53:35 -0700
Subject: [PATCH 20/28] fix, (goes with ConstantTAPhiDomain calcthermo) squash

---
 src/Domain.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/Domain.jl b/src/Domain.jl
index 4b63880fe..f8aefe805 100644
--- a/src/Domain.jl
+++ b/src/Domain.jl
@@ -1619,7 +1619,7 @@ end
         kfs = d.kfs 
         krevs = d.krevs
         Gs = d.Gs
-        return ns, cs, d.T, P, d.A, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, , similar(cs,0,0)
+        return ns, cs, d.T, P, d.A, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, similar(cs,0,0)
     end
     
 end

From c44d2fc397707394b486ac7a92fd002b30352e1c Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Thu, 5 Jun 2025 13:15:27 -0700
Subject: [PATCH 21/28] fix thermovec

---
 src/Calculators/Thermovec.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/Calculators/Thermovec.jl b/src/Calculators/Thermovec.jl
index 6cd0bd7b7..bc2a951b6 100644
--- a/src/Calculators/Thermovec.jl
+++ b/src/Calculators/Thermovec.jl
@@ -56,7 +56,7 @@ function NASAvec(nasas::B) where {B<:Array}
                 polyvec[:,j] = nasapoly.coefs
             end
         end
-        nasapvec = NASApolynomialvec(polyvec,Ts[i],Ts[i+1],nothing)
+        nasapvec = NASApolynomialvec(polyvec,Ts[i],Ts[i+1],covdeps)
         push!(polyvecs,nasapvec)
     end
     return NASAvec(polys=polyvecs,covdeps=covdeps)

From 011a35693c66cc4b221ceb9f5239e9b6abf51b00 Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Thu, 5 Jun 2025 13:15:57 -0700
Subject: [PATCH 22/28] handle no coverage inputs in coverage dependence

---
 src/Calculators/RateCoverageDependence.jl   | 6 ++++++
 src/Calculators/ThermoCoverageDependence.jl | 3 +++
 2 files changed, 9 insertions(+)

diff --git a/src/Calculators/RateCoverageDependence.jl b/src/Calculators/RateCoverageDependence.jl
index 942c6ddaa..ee85dbc4d 100644
--- a/src/Calculators/RateCoverageDependence.jl
+++ b/src/Calculators/RateCoverageDependence.jl
@@ -28,6 +28,9 @@ ms corresponds to m (float)
     indms::Dict{N,L} = Dict()
 end
 @inline function getcovdepactivationbarriercorrection(covdep::PolynomialRateCoverageDependence,T,coverages)
+    if length(coverages) == 0
+        return zero(typeof(coverages).parameters[1])
+    end
     E = 0.0
     for (ind,Epoly) in covdep.indEpolys
         E += evalpoly(coverages[ind],indEpoly)
@@ -35,6 +38,9 @@ end
     return E
 end
 @inline function getcovdepfactorcorrection(covdep::PolynomialRateCoverageDependence,T,coverages)
+    if length(coverages) == 0
+        return zero(typeof(coverages).parameters[1])
+    end
     av = 0.0
     for (ind,a) in covdep.indavals
        av += coverages[ind]*a
diff --git a/src/Calculators/ThermoCoverageDependence.jl b/src/Calculators/ThermoCoverageDependence.jl
index 29aa25645..5d5c10c1c 100644
--- a/src/Calculators/ThermoCoverageDependence.jl
+++ b/src/Calculators/ThermoCoverageDependence.jl
@@ -23,6 +23,9 @@ end
     return PolynomialThermoEnergyCoverageDependence(polys,Dict{Int64,Vector{Float64}}())
 end
 @inline function getcovdepenthalpycorrection(covdep::PolynomialThermoEnergyCoverageDependence,T,coverages)
+    if length(coverages) == 0
+        return zero(typeof(coverages).parameters[1])
+    end
     E = zero(coverages[1])
     for (ind,poly) in covdep.indpolys
         E += evalpoly(coverages[ind],poly)

From 4542b66e8075ac79103b47c57b6058b4cfdb7502 Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Thu, 5 Jun 2025 13:16:26 -0700
Subject: [PATCH 23/28] fix calcthermo for ConstantTAPhiDomain

---
 src/Domain.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/Domain.jl b/src/Domain.jl
index f8aefe805..860a49196 100644
--- a/src/Domain.jl
+++ b/src/Domain.jl
@@ -1634,7 +1634,7 @@ end
         if !d.alternativepformat
             @views nothermochg = d.Gs == p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
             if nothermochg
-                return ns, cs, d.T, P, d.A, C, N, d.mu, p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)], d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, , similar(cs,0,0)
+                return ns, cs, d.T, P, d.A, C, N, d.mu, p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)], d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, similar(cs,0,0)
             else
                 d.kfs = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
                 d.Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]

From f59eb2b122ebe202e617a652797eb29882b6b129 Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Thu, 5 Jun 2025 13:16:48 -0700
Subject: [PATCH 24/28] fix constantTPhi interface

---
 src/Interface.jl | 34 +++++++++++++++++-----------------
 1 file changed, 17 insertions(+), 17 deletions(-)

diff --git a/src/Interface.jl b/src/Interface.jl
index edb8d781c..3938b7dcc 100644
--- a/src/Interface.jl
+++ b/src/Interface.jl
@@ -127,13 +127,13 @@ export ReactiveInternalInterfaceConstantTPhi
 
 function getkfskrevs(ri::ReactiveInternalInterfaceConstantTPhi, T1, T2, phi1, phi2, Gs1, Gs2, coverages1, coverages2, cstot)
     if ri.iscovdep
-        if !(coverages1 === nothing) && coverages2 === nothing
-          coverages = coverages1
-        elseif !(coverages2 === nothing) && coverages1 === nothing
+        if length(coverages1) != 0 && length(coverages2) == 0
+            coverages = coverages1
+        elseif length(coverages2) != 0 && length(coverages1) == 0
             coverages = coverages2
-        elseif coverages1 === nothing && coverages2 === nothing 
-            coverages = nothing
-        else 
+        elseif length(coverages1) == 0 && length(coverages2) == 0 
+            coverages = coverages1
+        else
             throw(DomainError("Coverage dependence of an interface between two surfaces is indeterminate (which coverage should be used?)"))
         end
         kfs = getkf.(ri.reactions, nothing, ri.T, 0.0, 0.0, Ref([]), ri.A, ri.phi; coverages=coverages)
@@ -147,13 +147,13 @@ end
 
 function evaluate(ri::ReactiveInternalInterfaceConstantTPhi, dydt, domains, T1, T2, phi1, phi2, Gs1, Gs2, coverages1, coverages2, cstot, p::W) where {W<:SciMLBase.NullParameters}
     if ri.iscovdep
-        if !(coverages1 === nothing) && coverages2 === nothing
-          coverages = coverages1
-        elseif !(coverages2 === nothing) && coverages1 === nothing
+        if length(coverages1) != 0 && length(coverages2) == 0
+            coverages = coverages1
+        elseif length(coverages2) != 0 && length(coverages1) == 0
             coverages = coverages2
-        elseif coverages1 === nothing && coverages2 === nothing 
-            coverages = nothing
-        else 
+        elseif length(coverages1) == 0 && length(coverages2) == 0 
+            coverages = coverages1
+        else
             throw(DomainError("Coverage dependence of an interface between two surfaces is indeterminate (which coverage should be used?)"))
         end
         kfs = getkf.(ri.reactions, nothing, ri.T, 0.0, 0.0, Ref([]), ri.A, ri.phi; coverages=coverages)
@@ -167,13 +167,13 @@ end
 
 function evaluate(ri::ReactiveInternalInterfaceConstantTPhi, dydt, domains, T1, T2, phi1, phi2, Gs1, Gs2, coverages1, coverages2, cstot, p)
     if ri.iscovdep
-        if !(coverages1 === nothing) && coverages2 === nothing
+        if length(coverages1) != 0 && length(coverages2) == 0
             coverages = coverages1
-        elseif !(coverages2 === nothing) && coverages1 === nothing
+        elseif length(coverages2) != 0 && length(coverages1) == 0
             coverages = coverages2
-        elseif coverages1 === nothing && coverages2 === nothing 
-            coverages = nothing
-        else 
+        elseif length(coverages1) == 0 && length(coverages2) == 0 
+            coverages = coverages1
+        else
             throw(DomainError("Coverage dependence of an interface between two surfaces is indeterminate (which coverage should be used?)"))
         end
         kfs = getkf.(ri.reactions, nothing, ri.T, 0.0, 0.0, Ref([]), ri.A, ri.phi; coverages=coverages) .* p[ri.parameterindexes[1]:ri.parameterindexes[2]] 

From c9775b892f2193f7299c9b4d4142f6988f48ecae Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Fri, 6 Jun 2025 20:25:10 -0700
Subject: [PATCH 25/28] fix use of d.Gs when Gs has and needs to be computed

---
 src/Domain.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/Domain.jl b/src/Domain.jl
index 860a49196..4d897b1c2 100644
--- a/src/Domain.jl
+++ b/src/Domain.jl
@@ -1658,7 +1658,7 @@ end
         @fastmath @views hdivRT .+= p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)] ./ (R * d.T)
         @fastmath Gs = (hdivRT .- sdivR) * (R *d.T)
         kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, Array{Float64,1}(), d.A, d.phi; coverages=coverages)
-        return ns, cs, d.T, P, d.A, C, N, d.mu, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, coverages
+        return ns, cs, d.T, P, d.A, C, N, d.mu, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, coverages
     end
 end
 
@@ -1684,7 +1684,7 @@ end
         @fastmath @views hdivRT .+= p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)] ./ (R * d.T)
         @fastmath Gs = (hdivRT .- sdivR) * (R * d.T)
         kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, Array{Float64,1}(), d.A, d.phi; coverages=coverages)
-        return ns, cs, d.T, P, d.A, C, N, d.mu, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, coverages
+        return ns, cs, d.T, P, d.A, C, N, d.mu, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, coverages
     end
 end
 
@@ -1710,7 +1710,7 @@ end
         @fastmath @views hdivRT .+= p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)] ./ (R * T)
         @fastmath Gs = (hdivRT .- sdivR) * (R * d.T)
         kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, Array{Float64,1}(), d.A, d.phi; coverages=coverages)
-        return ns, cs, d.T, P, d.A, C, N, d.mu, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, coverages
+        return ns, cs, d.T, P, d.A, C, N, d.mu, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, coverages
     end
 end
 

From 920c090a107d3c9eff9a2d6fce6bb933d3b5fc9c Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Fri, 6 Jun 2025 20:25:48 -0700
Subject: [PATCH 26/28] ensure coverage dependent indices get assigned for
 interface reactions

---
 src/Interface.jl | 28 +++++++++++++++++++++++++++-
 1 file changed, 27 insertions(+), 1 deletion(-)

diff --git a/src/Interface.jl b/src/Interface.jl
index 3938b7dcc..42f0da98e 100644
--- a/src/Interface.jl
+++ b/src/Interface.jl
@@ -35,6 +35,7 @@ struct ReactiveInternalInterface{T,B,C,C2,N,Q<:AbstractReaction,X} <: AbstractRe
     forwardability::Array{Bool,1}
 end
 function ReactiveInternalInterface(domain1, domain2, reactions, A)
+    reactions = upgradekinetics(reactions, domain1, domain2)
     vectuple, vecinds, otherrxns, otherrxninds, posinds = getveckinetics(reactions)
     rxns = vcat(reactions[vecinds], reactions[otherrxninds])
     rxns = [ElementaryReaction(index=i, reactants=rxn.reactants, reactantinds=rxn.reactantinds, products=rxn.products,
@@ -269,6 +270,31 @@ function upgradekinetics(rxns, domain1, domain2)
     end
     newrxns = Array{ElementaryReaction,1}(undef, length(rxns))
     for (i, rxn) in enumerate(rxns)
+        if hasproperty(rxn.kinetics,:covdep) && !(rxn.kinetics.covdep isa EmptyRateCoverageDependence)
+            if domain1.phase isa IdealSurface && !(domain2.phase isa IdealSurface)
+                names = getfield.(domain1.phase.species,:name)
+            elseif domain2.phase isa IdealSurface && !(domain1.phase isa IdealSurface)
+                names = getfield.(domain2.phase.species,:name)
+            else 
+                throw(DomainError(rxn.kinetics.covdep,"Unsure which surface domain should be used for coverage dependence"))
+            end
+            if rxn.kinetics.covdep isa PolynomialRateCoverageDependence
+                for (name,v) in rxn.kinetics.covdep.Epolys
+                    ind = findfirst(isequal(name),names)
+                    rxn.kinetics.covdep.indEpolys[ind] = v
+                end 
+                for (name,v) in rxn.kinetics.covdep.avals
+                    ind = findfirst(isequal(name),names)
+                    rxn.kinetics.covdep.indavals[ind] = v
+                end 
+                for (name,v) in rxn.kinetics.covdep.ms
+                    ind = findfirst(isequal(name),names)
+                    rxn.kinetics.covdep.indms[ind] = v
+                end 
+            else
+                throw(TypeError(rxn.kinetics.covdep,"Kinetic Coverage Dependence Type Not Understood"))
+            end
+        end
         if isa(rxn.kinetics, StickingCoefficient)
             spc = [spc for spc in rxn.reactants if !(spc in surfdomain.phase.species)]
             @assert length(spc) == 1
@@ -293,7 +319,7 @@ function stickingcoefficient2arrhenius(sc, sitedensity, N, mw; Tmin=300.0, Tmax=
     @assert fit.converged
     p = fit.param
     p[1] = abs(p[1])
-    return Arrhenius(; A=p[1], n=p[2], Ea=p[3])
+    return Arrhenius(; A=p[1], n=p[2], Ea=p[3], covdep=sc.covdep)
 end
 
 struct Inlet{Q<:Real,S,V<:AbstractArray,U<:Real,X<:Real,FF<:Function} <: AbstractBoundaryInterface

From b8c3be89fbd36c12bf5562a2fc5a080ca20bfa0a Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Fri, 6 Jun 2025 20:26:13 -0700
Subject: [PATCH 27/28] fix rate coverage dependence

---
 src/Calculators/RateCoverageDependence.jl | 43 ++++++++++++++---------
 1 file changed, 27 insertions(+), 16 deletions(-)

diff --git a/src/Calculators/RateCoverageDependence.jl b/src/Calculators/RateCoverageDependence.jl
index ee85dbc4d..f88f446b9 100644
--- a/src/Calculators/RateCoverageDependence.jl
+++ b/src/Calculators/RateCoverageDependence.jl
@@ -9,9 +9,6 @@ export EmptyRateCoverageDependence
 getcovdepactivationbarriercorrection(covdep::EmptyRateCoverageDependence,T,coverages) = 0.0
 getcovdepfactorcorrection(covdep::EmptyRateCoverageDependence,T,coverages) = 1.0
 
-export getcovdepactivationbarriercorrection
-export getcovdepfactorcorrection
-
 """
 Polynomial rate coefficient coverage dependence
 correction form is: 10^(sum(a_i*theta_k^i,i)) * theta_k^m * exp(sum(E_i*theta_k^i,i)/(RT))
@@ -19,27 +16,38 @@ dictionaries map the name of the species (k) to polynomials with ordering a_1 +
 Epolys corresponds to E_i (array of polynomial coefficients), avals corresponds to a_i (array of polynomial coefficients), 
 ms corresponds to m (float)
 """
-@with_kw struct PolynomialRateCoverageDependence{N,K,L} <: AbstractRateCoverageDependence
-    avals::Dict{String,L} = Dict()
-    Epolys::Dict{String,K} = Dict()
-    ms::Dict{String,L} = Dict()
-    indavals::Dict{N,L} = Dict()
-    indEpolys::Dict{N,K} = Dict()
-    indms::Dict{N,L} = Dict()
+@with_kw struct PolynomialRateCoverageDependence <: AbstractRateCoverageDependence
+    avals::Dict{String,Vector{Float64}} = Dict{String,Vector{Float64}}()
+    Epolys::Dict{String,Vector{Float64}} = Dict{String,Vector{Float64}}()
+    ms::Dict{String,Vector{Float64}} = Dict{String,Float64}()
+    indavals::Dict{Int64,Vector{Float64}} = Dict{Int64,Vector{Float64}}()
+    indEpolys::Dict{Int64,Vector{Float64}} = Dict{Int64,Vector{Float64}}()
+    indms::Dict{Int64,Vector{Float64}} = Dict{Int64,Float64}()
+end
+@inline function PolynomialRateCoverageDependence(Epolys,ms=Dict{String,Float64}(),avals=Dict{String,Vector{Float64}}())
+    return PolynomialRateCoverageDependence(;Epolys=Epolys,
+            ms=ms,avals=avals,indavals=Dict{Int64,Vector{Float64}}(),indEpolys=Dict{Int64,Vector{Float64}}(),
+            indms=Dict{Int64,Float64}())
 end
 @inline function getcovdepactivationbarriercorrection(covdep::PolynomialRateCoverageDependence,T,coverages)
-    if length(coverages) == 0
+    if coverages === nothing
+        return 0.0
+    elseif length(coverages) == 0
         return zero(typeof(coverages).parameters[1])
     end
     E = 0.0
+    
     for (ind,Epoly) in covdep.indEpolys
-        E += evalpoly(coverages[ind],indEpoly)
+        E += evalpoly(coverages[ind],Epoly)
     end
     return E
 end
+export getcovdepactivationbarriercorrection
 @inline function getcovdepfactorcorrection(covdep::PolynomialRateCoverageDependence,T,coverages)
-    if length(coverages) == 0
-        return zero(typeof(coverages).parameters[1])
+    if coverages === nothing
+        return 1.0
+    elseif length(coverages) == 0
+        return one(typeof(coverages).parameters[1])
     end
     av = 0.0
     for (ind,a) in covdep.indavals
@@ -55,7 +63,10 @@ end
     for (ind,m) in covdep.indms 
         A *= coverages[ind]^m
     end
-    
+
     return A
 end
-export PolynomialRateCoverageDependence
\ No newline at end of file
+export PolynomialRateCoverageDependence
+
+export getcovdepactivationbarriercorrection
+export getcovdepfactorcorrection
\ No newline at end of file

From 1597d83781892ae14aaf4d47fa15dfde5a8c9ca0 Mon Sep 17 00:00:00 2001
From: Matt Johnson <mjohnson541@gmail.com>
Date: Fri, 6 Jun 2025 20:26:28 -0700
Subject: [PATCH 28/28] ensure getkfs is called with coverages

---
 src/PhaseState.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/PhaseState.jl b/src/PhaseState.jl
index 821136447..f0922d0d0 100644
--- a/src/PhaseState.jl
+++ b/src/PhaseState.jl
@@ -225,7 +225,7 @@ export getkfkrev
 
 @inline function getkfkrevs(phase::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Q2,diffs::Q3,V::W5,phi::W7;coverages=nothing,kfs::W6=nothing) where {U,W7,W6,W5<:Real,W1<:Real,W2<:Real,W3,W4,Q1<:AbstractArray,Q2,Q3<:AbstractArray}
     if !phase.diffusionlimited && kfs === nothing
-        kfs = getkfs(phase,T,P,C,ns,V,phi)
+        kfs = getkfs(phase,T,P,C,ns,V,phi;coverages=coverages)
         if phi == 0.0
             krev = @fastmath kfs./getKcs(phase,T,Gs)
         else 
@@ -258,7 +258,7 @@ end
 
 @inline function getkfkrevs(phase::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Q2,diffs::Q3,V::W5,phi::W7;coverages=nothing,kfs::W6=nothing) where {U,W7,W6,W5<:Real,W1<:Real,W2<:Real,W3,W4,Q1<:AbstractArray,Q2<:Union{ReverseDiff.TrackedArray,Tracker.TrackedArray},Q3<:AbstractArray} #autodiff p
     if !phase.diffusionlimited && kfs === nothing
-        kfs = getkfs(phase,T,P,C,ns,V,phi)
+        kfs = getkfs(phase,T,P,C,ns,V,phi;coverages=coverages)
         if phi == 0.0
             krev = @fastmath kfs./getKcs(phase,T,Gs)
         else