Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Modify the Marcus branch to allow PCET reactions #263

Closed
2 changes: 1 addition & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ RMS has been used in many applications:
## How to cite
Please include the following citations for ReactionMechanismSimulator.jl in general and for transitory sensitivities and the automatic mechanism analysis toolkit respectively.

- <div class="csl-entry">Johnson, M. S., Pang, H.-W., Payne, A. M., &#38; Green, W. H. (2023). <i>ReactionMechanismSimulator.jl: A Modern Approach to Chemical Kinetic Mechanism Simulation and Analysis</i>. https://doi.org/10.26434/CHEMRXIV-2023-TJ34T</div>
- <div class="csl-entry">Johnson, M. S., Pang, H.-W., Payne, A. M., &#38; Green, W. H. (2023). <i>ReactionMechanismSimulator.jl: A Modern Approach to Chemical Kinetic Mechanism Simulation and Analysis</i>. Int. J. Chem. Kinet. 2024;1-16 https://doi.org/10.1002/kin.21753</div>
- <div class="csl-entry">Johnson, M. S., McGill, C. J., &#38; Green, W. H. (2022). <i>Transitory Sensitivity in Automatic Chemical Kinetic Mechanism Analysis</i>. https://doi.org/10.26434/CHEMRXIV-2022-ZSFJC</div>

## Installation
Expand Down
51 changes: 36 additions & 15 deletions src/Calculators/Rate.jl
Original file line number Diff line number Diff line change
Expand Up @@ -13,8 +13,8 @@ export AbstractFalloffRate
Ea::Q
unc::P = EmptyRateUncertainty()
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,dGrxn=0.0,d=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,d=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp(-arr.Ea/(R*T))
export Arrhenius

@with_kw struct StickingCoefficient{N<:Real,K<:Real,Q<:Real,P<:AbstractRateUncertainty} <: AbstractRate
Expand All @@ -23,29 +23,50 @@ export Arrhenius
Ea::Q
unc::P = EmptyRateUncertainty()
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,dGrxn=0.0,d=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,dGrxn=0.0,d=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath min(arr.A*T^arr.n*exp(-arr.Ea/(R*T)),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,R<:Real,P<:AbstractRateUncertainty,B} <: AbstractRate
A::N
n::K
Ea::Q
q::B = 0.0
V0::R = 0.0
unc::P = EmptyRateUncertainty()
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,dGrxn=0.0,d=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp((-arr.Ea-arr.q*F*(phi-arr.V0))/(R*T))
@inline (arr::Arrheniusq)(T::Q;P::N=0.0,C::S=0.0,phi=0.0,dGrxn=0.0,d=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp((-arr.Ea-arr.q*F*(phi-arr.V0))/(R*T))
export Arrheniusq

@with_kw struct Marcus{N<:Real,K<:Real,Q,P<:AbstractRateUncertainty,B} <: AbstractRate
A::N
n::K
lmbd_i_coefs::Q
lmbd_o::K
wr::K
wp::K
beta::B
unc::P = EmptyRateUncertainty()
end
@inline function (arr::Marcus)(;T::Q,P::N=0.0,C::S=0.0,phi=0.0,dGrxn=0.0,d=0.0) where {Q<:Real,N<:Real,S<:Real}
@fastmath lmbd = arr.lmbd_o + evalpoly(T,arr.lmbd_i_coefs)
@fastmath arr.A*T^arr.n*exp(-lmbd/4.0*(1.0+dGrxn/lmbd)^2/(R*T)-arr.beta*d)
end
@inline function (arr::Marcus)(T::Q;P::N=0.0,C::S=0.0,phi=0.0,dGrxn=0.0,d=0.0) where {Q<:Real,N<:Real,S<:Real}
@fastmath lmbd = arr.lmbd_o + evalpoly(T,arr.lmbd_i_coefs)
@fastmath return arr.A*T^arr.n*exp(-lmbd/4.0*(1.0+dGrxn/lmbd)^2/(R*T)-arr.beta*d)
end
export Marcus

@with_kw struct PdepArrhenius{T<:Real,Q<:AbstractRateUncertainty,Z<:AbstractRate} <: AbstractRate
Ps::Array{T,1}
arrs::Array{Z,1}
unc::Q = EmptyRateUncertainty()
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,dGrxn=0.0,d=0.0) where {Q<:Real,V<:Real,S<:Real}
inds = getBoundingIndsSorted(P,parr.Ps)::Tuple{Int64,Int64}

if inds[2] == -1
Expand All @@ -65,7 +86,7 @@ export PdepArrhenius
unc::Q = EmptyRateUncertainty()
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,dGrxn=0.0,d=0.0) where {Q<:Real,R<:Real,S<:Real}
out = 0.0
for arr in marr.arrs
@fastmath out += arr(T)
Expand All @@ -79,7 +100,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,dGrxn=0.0,d=0.0) where {Q<:Real,R<:Real,S<:Real}
out = 0.0
for pdar in parr.parrs
@fastmath out += pdar(T=T,P=P)
Expand All @@ -95,7 +116,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,dGrxn=0.0,d=0.0) where {Q<:Real,R<:Real,S<:Real} = C*(tbarr.arr(T))
export ThirdBody

@with_kw struct Lindemann{N<:Integer,K<:AbstractFloat,Q<:AbstractRateUncertainty} <: AbstractFalloffRate
Expand All @@ -106,7 +127,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,dGrxn=0.0,d=0.0) where {Q<:Real,R<:Real,S<:Real}
k0 = lnd.arrlow(T=T)
kinf = lnd.arrhigh(T=T)
@fastmath Pr = k0*C/kinf
Expand All @@ -126,7 +147,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,dGrxn=0.0,d=0.0) where {Q<:Real,R<:Real,S<:Real}
k0 = tr.arrlow(T=T)
kinf = tr.arrhigh(T=T)
@fastmath Pr = k0*C/kinf
Expand Down Expand Up @@ -195,7 +216,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,dGrxn=0.0,d=0.0) where {N<:Real,B<:Real,Q<:Real}
k = 0.0
Tred = getredtemp(ch,T)
Pred = getredpress(ch,P)
Expand Down Expand Up @@ -228,7 +249,7 @@ export getkineticstype

@inline extracttypename(typ::Symbol) = string(typ)
@inline extracttypename(typ) = string(typ.name)

@inline function _calcdkdCeff(tbarr::ThirdBody,T::Float64,Ceff::Float64)
return @fastmath tbarr.arr(T)
end
Expand Down
8 changes: 4 additions & 4 deletions src/Calculators/Ratevec.jl
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ 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,dGrxns=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)))
export Arrheniusvec

Expand Down Expand Up @@ -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,dGrxns=0.0) where {N<:Real,B<:Real,Q<:Real}
k = zeros(N,size(ch.coefs)[1])
Tred = getredtemp(ch,T)
Pred = getredpress(ch,P)
Expand Down Expand Up @@ -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,dGrxns=0.0) where {Q<:Real,R<:Real,S<:Real}
k0 = tr.arrlow(T=T)
kinf = tr.arrhigh(T=T)
@fastmath Pr = k0.*C./kinf
Expand Down Expand Up @@ -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,dGrxns=0.0) 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)
Expand Down
Loading