Skip to content

Commit c88e830

Browse files
mjohnson541mjohnson541
committed
handle dGrxn and d dependence of reactions in interfaces and phase calculations
1 parent 3a6e483 commit c88e830

File tree

2 files changed

+88
-130
lines changed

2 files changed

+88
-130
lines changed

‎src/Interface.jl

Lines changed: 33 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -34,26 +34,28 @@ struct ReactiveInternalInterface{T,B,C,C2,N,Q<:AbstractReaction,X} <: AbstractRe
3434
reversibililty::Array{Bool,1}
3535
forwardability::Array{Bool,1}
3636
end
37-
function ReactiveInternalInterface(domain1, domain2, reactions, A)
38-
vectuple, vecinds, otherrxns, otherrxninds, posinds = getveckinetics(reactions)
39-
rxns = vcat(reactions[vecinds], reactions[otherrxninds])
40-
rxns = [ElementaryReaction(index=i, reactants=rxn.reactants, reactantinds=rxn.reactantinds, products=rxn.products,
41-
productinds=rxn.productinds, kinetics=rxn.kinetics, radicalchange=rxn.radicalchange, reversible=rxn.reversible, forwardable=rxn.forwardable, pairs=rxn.pairs) for (i, rxn) in enumerate(rxns)]
42-
rxnarray = getinterfacereactioninds(domain1, domain2, rxns)
43-
M, Nrp1, Nrp2 = getstoichmatrix(domain1, domain2, reactions)
44-
reversibility = Array{Bool,1}(getfield.(rxns, :reversible))
45-
forwardability = Array{Bool,1}(getfield.(rxns, :forwardable))
46-
return ReactiveInternalInterface(domain1, domain2,
47-
rxns, vectuple, posinds, rxnarray, M, Nrp1, Nrp2, A, [1, length(reactions)],
48-
[0, 1], ones(length(rxns)), reversibility, forwardability), ones(length(rxns))
37+
function ReactiveInternalInterface(domain1,domain2,reactions,A)
38+
vectuple,vecinds,otherrxns,otherrxninds,posinds = getveckinetics(reactions)
39+
rxns = vcat(reactions[vecinds],reactions[otherrxninds])
40+
rxns = [ElementaryReaction(index=i,reactants=rxn.reactants,reactantinds=rxn.reactantinds,products=rxn.products,
41+
productinds=rxn.productinds,kinetics=rxn.kinetics,electronchange=rxn.electronchange,radicalchange=rxn.radicalchange,reversible=rxn.reversible,forwardable=rxn.forwardable,pairs=rxn.pairs) for (i,rxn) in enumerate(rxns)]
42+
rxnarray = getinterfacereactioninds(domain1,domain2,rxns)
43+
M,Nrp1,Nrp2 = getstoichmatrix(domain1,domain2,reactions)
44+
reversibility = Array{Bool,1}(getfield.(rxns,:reversible))
45+
forwardability = Array{Bool,1}(getfield.(rxns,:forwardable))
46+
return ReactiveInternalInterface(domain1,domain2,
47+
rxns,vectuple,posinds,rxnarray,M,Nrp1,Nrp2,A,[1,length(reactions)],
48+
[0,1],ones(length(rxns)),reversibility,forwardability),ones(length(rxns))
4949
end
5050
export ReactiveInternalInterface
5151

52-
function getkfskrevs(ri::ReactiveInternalInterface, T1, T2, phi1, phi2, Gs1, Gs2, cstot::Array{Q,1}) where {Q}
53-
kfs = getkfs(ri, T1, 0.0, 0.0, Array{Q,1}(), ri.A, phi1)
54-
Kc = getKc.(ri.reactions, ri.domain1.phase, ri.domain2.phase, Ref(Gs1), Ref(Gs2), T1, phi1)
55-
krevs = kfs ./ Kc
56-
return kfs, krevs
52+
function getkfskrevs(ri::ReactiveInternalInterface,T1,T2,phi1,phi2,Gs1,Gs2,cstot::Array{Q,1}) where {Q}
53+
Gpart = ArrayPartition(Gs1,Gs2)
54+
dGrxn = -ri.stoichmatrix*Gpart
55+
kfs = getkfs(ri,T1,0.0,0.0,Array{Q,1}(),ri.A,phi1,dGrxns,0.0)
56+
Kc = getKc.(ri.reactions,ri.domain1.phase,ri.domain2.phase,Ref(Gs1),Ref(Gs2),dGrxns,T1,phi1)
57+
krevs = kfs./Kc
58+
return kfs,krevs
5759
end
5860

5961
function evaluate(ri::ReactiveInternalInterface, dydt, domains, T1, T2, phi1, phi2, Gs1, Gs2, cstot, p::W) where {W<:SciMLBase.NullParameters}
@@ -86,18 +88,20 @@ struct ReactiveInternalInterfaceConstantTPhi{J,N,B,B2,B3,C,C2,Q<:AbstractReactio
8688
reversibility::Array{Bool,1}
8789
forwardability::Array{Bool,1}
8890
end
89-
function ReactiveInternalInterfaceConstantTPhi(domain1, domain2, reactions, T, A, phi=0.0)
90-
@assert domain1.T == domain2.T
91-
reactions = upgradekinetics(reactions, domain1, domain2)
92-
rxnarray = getinterfacereactioninds(domain1, domain2, reactions)
93-
kfs = getkf.(reactions, nothing, T, 0.0, 0.0, Ref([]), A, phi)
94-
Kc = getKc.(reactions, domain1.phase, domain2.phase, Ref(domain1.Gs), Ref(domain2.Gs), T, phi)
95-
krevs = kfs ./ Kc
96-
M, Nrp1, Nrp2 = getstoichmatrix(domain1, domain2, reactions)
97-
reversibility = Array{Bool,1}(getfield.(reactions, :reversible))
98-
forwardability = Array{Bool,1}(getfield.(reactions, :forwardable))
99-
if isa(reactions, Vector{Any})
100-
reactions = convert(Vector{ElementaryReaction}, reactions)
91+
function ReactiveInternalInterfaceConstantTPhi(domain1,domain2,reactions,T,A,phi=0.0)
92+
@assert domain1.T == domain2.T
93+
reactions = upgradekinetics(reactions,domain1,domain2)
94+
rxnarray = getinterfacereactioninds(domain1,domain2,reactions)
95+
M,Nrp1,Nrp2 = getstoichmatrix(domain1,domain2,reactions)
96+
Gpart = ArrayPartition(domain1.Gs,domain2.Gs)
97+
dGrxns = -M*Gpart
98+
kfs = getkf.(reactions,nothing,T,0.0,0.0,Ref([]),A,phi,dGrxns,0.0)
99+
Kc = getKcs(domain1.phase,domain2.phase,T,Nrp1,Nrp2,dGrxns)
100+
krevs = kfs./Kc
101+
reversibility = Array{Bool,1}(getfield.(reactions,:reversible))
102+
forwardability = Array{Bool,1}(getfield.(reactions,:forwardable))
103+
if isa(reactions,Vector{Any})
104+
reactions = convert(Vector{ElementaryReaction},reactions)
101105
end
102106
if isa(kfs, Vector{Any})
103107
kfs = convert(Vector{Float64}, kfs)

‎src/PhaseState.jl

Lines changed: 55 additions & 101 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@ using LinearAlgebra
44
using Tracker
55
using ReverseDiff
66
using RecursiveArrayTools
7+
using Logging
78

89
@inline function calcgibbs(ph::U,T::W) where {U<:IdealPhase,W<:Real}
910
return getGibbs.(getfield.(ph.species,:thermo),T)
@@ -50,15 +51,15 @@ end
5051

5152
export makespcsvector
5253

53-
@inline function getkf(rxn::ElementaryReaction,ph,T,P,C,ns,V,phi)
54+
@inline function getkf(rxn::ElementaryReaction,ph,T,P,C,ns,V,phi,dGrxn,d)
5455
if isdefined(rxn.kinetics,:efficiencies) && length(rxn.kinetics.efficiencies) > 0
5556
@views @inbounds @fastmath C += sum([ns[i]*val for (i,val) in rxn.kinetics.efficiencies])/V
5657
end
57-
return rxn.kinetics(T=T,P=P,C=C,phi=phi)
58+
return rxn.kinetics(T=T,P=P,C=C,phi=phi,dGrxn=dGrxn,d=d)
5859
end
5960
export getkf
6061

61-
@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}
62+
@inline function getkfs(ph::U,T::W1,P::W2,C::W3,ns::Q,V::W4,phi,dGrxns,d) where {U,W1,W2,W3,W4<:Real,Q<:AbstractArray}
6263
kfs = similar(ns,length(ph.reactions))
6364
i = 1
6465
oldind = 1
@@ -70,7 +71,7 @@ export getkf
7071
i += 1
7172
end
7273
@simd for i in ind+1:length(ph.reactions)
73-
@inbounds kfs[i] = getkf(ph.reactions[i],ph,T,P,C,ns,V,phi)
74+
@inbounds kfs[i] = getkf(ph.reactions[i],ph,T,P,C,ns,V,phi,dGrxns[i],d)
7475
end
7576
return kfs
7677
end
@@ -84,13 +85,13 @@ for 2 spc calculates using the Smolchowski equation
8485
for >2 spc calculates using the Generalized Smolchowski equation
8586
Equations from Flegg 2016
8687
"""
87-
@inline function getDiffusiveRate(spcs::Q,diffs::Array{W,1}) where {Q<:AbstractArray,W<:Real}
88+
@inline function getDiffusiveRate(spcs::Q,spcsinds::Q2,diffs::Array{W,1}) where {Q<:AbstractArray,W<:Real,Q2<:AbstractArray}
8889
if length(spcs) == 1
8990
return Inf
9091
elseif length(spcs) == 2
91-
@fastmath @inbounds kf = 4.0*Base.pi*(diffs[spcs[1].index]+diffs[spcs[2].index])*(spcs[1].radius+spcs[2].radius)*Na
92+
@fastmath @inbounds kf = 4.0*Base.pi*(diffs[spcsinds[1]]+diffs[spcsinds[2]])*(spcs[1].radius+spcs[2].radius)*Na
9293
else
93-
@views @inbounds diffusivity = diffs[getfield.(spcs,:index)]
94+
@views @inbounds diffusivity = diffs[spcsinds]
9495
N = length(spcs)
9596
@fastmath a = (3.0*length(spcs)-5.0)/2.0
9697
@fastmath Dinv = 1.0./diffusivity
@@ -103,76 +104,26 @@ Equations from Flegg 2016
103104
end
104105
export getDiffusiveRate
105106

106-
@inline function getKc(rxn::ElementaryReaction,ph::U,T::Z,Gs::Q,phi::V=0.0) where {U<:AbstractPhase,V,Q,Z<:Real}
107+
@inline function getKc(rxn::ElementaryReaction,ph::U,T::Z,dGrxn::Q,phi::V=0.0) where {U<:AbstractPhase,V,Q,Z<:Real}
107108
Nreact = length(rxn.reactantinds)
108109
Nprod = length(rxn.productinds)
109-
dGrxn = 0.0
110-
if Nreact == 1
111-
@fastmath @inbounds dGrxn -= Gs[rxn.reactantinds[1]]
112-
elseif Nreact == 2
113-
@fastmath @inbounds dGrxn -= Gs[rxn.reactantinds[1]]+Gs[rxn.reactantinds[2]]
114-
elseif Nreact == 3
115-
@fastmath @inbounds dGrxn -= Gs[rxn.reactantinds[1]]+Gs[rxn.reactantinds[2]]+Gs[rxn.reactantinds[3]]
116-
elseif Nreact == 4
117-
@fastmath @inbounds dGrxn -= Gs[rxn.reactantinds[1]]+Gs[rxn.reactantinds[2]]+Gs[rxn.reactantinds[3]]+Gs[rxn.reactantinds[4]]
118-
end
119-
if Nprod == 1
120-
@fastmath @inbounds dGrxn += Gs[rxn.productinds[1]]
121-
elseif Nprod == 2
122-
@fastmath @inbounds dGrxn += Gs[rxn.productinds[1]]+Gs[rxn.productinds[2]]
123-
elseif Nprod == 3
124-
@fastmath @inbounds dGrxn += Gs[rxn.productinds[1]]+Gs[rxn.productinds[2]]+Gs[rxn.productinds[3]]
125-
elseif Nprod == 4
126-
@fastmath @inbounds dGrxn += Gs[rxn.productinds[1]]+Gs[rxn.productinds[2]]+Gs[rxn.productinds[3]]+Gs[rxn.productinds[4]]
127-
end
128-
return @inbounds @fastmath exp(-(dGrxn+rxn.electronchange*phi)/(R*T))*(getC0(ph,T))^(Nprod-Nreact)
110+
return @inbounds @fastmath exp((dGrxn+rxn.electronchange*(phi*F))/(-R*T))*(getC0(ph,T))^(Nprod-Nreact)
129111
end
130112

131-
@inline function getKc(rxn::ElementaryReaction,phase1,phase2,Gs1,Gs2,T,phi=0.0) #for constant k interfaces
132-
dGrxn = 0.0
133-
dN1 = 0
134-
dN2 = 0
135-
for r in rxn.reactants
136-
isfirst = true
137-
ind = findfirst(isequal(r),phase1.species)
138-
if ind === nothing
139-
isfirst = false
140-
ind = findfirst(isequal(r),phase2.species)
141-
dGrxn -= Gs2[ind]
142-
dN2 -= 1
143-
else
144-
dGrxn -= Gs1[ind]
145-
dN1 -= 1
146-
end
147-
end
148-
for r in rxn.products
149-
isfirst = true
150-
ind = findfirst(isequal(r),phase1.species)
151-
if ind === nothing
152-
isfirst = false
153-
ind = findfirst(isequal(r),phase2.species)
154-
dGrxn += Gs2[ind]
155-
dN2 += 1
156-
else
157-
dGrxn += Gs1[ind]
158-
dN1 += 1
159-
end
160-
end
161-
return @inbounds @fastmath exp(-(dGrxn+rxn.electronchange*phi)/(R*T))*getC0(phase1,T)^dN1*getC0(phase2,T)^dN2
113+
@inline function getKcs(ph::U,T::Z,dGrxns::Q) where {U<:AbstractPhase,Q,Z<:Real}
114+
return @fastmath @inbounds exp.(dGrxns./(-R*T) .+ ph.Nrp.*log(getC0(ph,T)));
162115
end
163-
export getKc
164116

165-
@inline function getKcs(ph::U,T::Z,Gs::Q) where {U<:AbstractPhase,Q,Z<:Real}
166-
return @fastmath @inbounds exp.(ph.stoichmatrix*(Gs./(R*T)) .+ ph.Nrp.*log(getC0(ph,T)));
117+
@inline function getKcs(ph::U,T::Z,dGrxns::Q,phi::V) where {U<:AbstractPhase,Q,Z<:Real,V<:Real}
118+
return @fastmath @inbounds exp.(dGrxns./(-R*T) .+ ph.Nrp.*log(getC0(ph,T)));
167119
end
168120

169-
@inline function getKcs(ph::U,T::Z,Gs::Q,phi::V) where {U<:AbstractPhase,Q,Z<:Real,V<:Real}
170-
return @fastmath @inbounds exp.(ph.stoichmatrix*(Gs./(R*T)).+ph.electronchange.*(phi/(R*T)) .+ ph.Nrp.*log(getC0(ph,T)));
121+
@inline function getKcs(ph,T,Gs1,Gs2,dGrxns)
122+
return @fastmath @inbounds exp.(dGrxns./(-R*T) .+ ph.Nrp1.*log(getC0(ph.domain1.phase,T)) .+ ph.Nrp2.*log(getC0(ph.domain2.phase,T)));
171123
end
172124

173-
@inline function getKcs(ph,T,Gs1,Gs2)
174-
Gpart = ArrayPartition(Gs1,Gs2)
175-
return @fastmath @inbounds exp.(ph.stoichmatrix*(Gpart./(R*T)) .+ ph.Nrp1.*log(getC0(ph.domain1.phase,T)) .+ ph.Nrp2.*log(getC0(ph.domain2.phase,T)));
125+
@inline function getKcs(ph1,ph2,T,Nrp1,Nrp2,dGrxns)
126+
return @fastmath @inbounds exp.(dGrxns./(-R*T) .+ Nrp1.*log(getC0(ph1,T)) .+ Nrp2.*log(getC0(ph2,T)));
176127
end
177128

178129
export getKcs
@@ -181,30 +132,30 @@ export getKcs
181132
Calculates the forward and reverse rate coefficients for a given reaction, phase and state
182133
Maintains diffusion limitations if the phase has diffusionlimited=true
183134
"""
184-
@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}
185-
if signbit(kf)
135+
@inline function getkfkrev(rxn::ElementaryReaction,ph::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,dGrxn::Q2,diffs::Q3,V::W5,phi::W8,d;kf::W6=-1.0,f::W7=-1.0) where {U<:AbstractPhase,W8,W6,W7,W5,W4,W1,W2,W3<:Real,Q1,Q2,Q3<:AbstractArray}
136+
if signbit(kf)
186137
if signbit(f)
187-
kf = getkf(rxn,ph,T,P,C,ns,V,phi)
138+
kf = getkf(rxn,ph,T,P,C,ns,V,phi,dGrxn,d)
188139
else
189-
kf = getkf(rxn,ph,T,P,C,ns,V,phi)*f
140+
kf = getkf(rxn,ph,T,P,C,ns,V,phi,dGrxn,d)*f
190141
end
191142
end
192-
Kc = getKc(rxn,ph,T,Gs,phi)
143+
Kc = getKc(rxn,ph,T,dGrxn,phi)
193144
@fastmath krev = kf/Kc
194145
if ph.diffusionlimited
195146
if length(rxn.reactants) == 1
196147
if length(rxn.products) > 1
197-
krevdiff = getDiffusiveRate(rxn.products,diffs)
148+
krevdiff = getDiffusiveRate(rxn.products,rxn.productinds,diffs)
198149
@fastmath krev = krev*krevdiff/(krev+krevdiff)
199150
@fastmath kf = Kc*krev
200151
end
201152
elseif length(rxn.products) == 1
202-
kfdiff = getDiffusiveRate(rxn.reactants,diffs)
153+
kfdiff = getDiffusiveRate(rxn.reactants,rxn.reactantinds,diffs)
203154
@fastmath kf = kf*kfdiff/(kf+kfdiff)
204155
@fastmath krev = kf/Kc
205156
elseif length(rxn.products) == length(rxn.reactants)
206-
kfdiff = getDiffusiveRate(rxn.reactants,diffs)
207-
krevdiff = getDiffusiveRate(rxn.products,diffs)
157+
kfdiff = getDiffusiveRate(rxn.reactants,rxn.reactantinds,diffs)
158+
krevdiff = getDiffusiveRate(rxn.products,rxn.productinds,diffs)
208159
@fastmath kff = kf*kfdiff/(kf+kfdiff)
209160
@fastmath krevr = krev*krevdiff/(krev+krevdiff)
210161
@fastmath kfr = Kc*krevr
@@ -223,98 +174,101 @@ Maintains diffusion limitations if the phase has diffusionlimited=true
223174
end
224175
export getkfkrev
225176

226-
@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}
177+
@inline function getkfkrevs(phase::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Q2,diffs::Q3,V::W5,phi::W7,d;kfs::W6=nothing) where {U,W7,W6,W5<:Real,W1<:Real,W2<:Real,W3,W4,Q1<:AbstractArray,Q2,Q3<:AbstractArray}
178+
dGrxns = -(phase.stoichmatrix*Gs).+phase.electronchange.*(phi*F)
227179
if !phase.diffusionlimited && kfs === nothing
228-
kfs = getkfs(phase,T,P,C,ns,V,phi)
180+
kfs = getkfs(phase,T,P,C,ns,V,phi,dGrxns,d)
229181
if phi == 0.0
230-
krev = @fastmath kfs./getKcs(phase,T,Gs)
182+
krev = @fastmath kfs./getKcs(phase,T,dGrxns)
231183
else
232-
krev = @fastmath kfs./getKcs(phase,T,Gs,phi)
184+
krev = @fastmath kfs./getKcs(phase,T,dGrxns,phi)
233185
end
234186
elseif !phase.diffusionlimited && !(kfs === nothing)
235187
if phi == 0.0
236-
krev = @fastmath kfs./getKcs(phase,T,Gs)
188+
krev = @fastmath kfs./getKcs(phase,T,dGrxns)
237189
else
238-
krev = @fastmath kfs./getKcs(phase,T,Gs,phi)
190+
krev = @fastmath kfs./getKcs(phase,T,dGrxns,phi)
239191
end
240192
elseif phase.diffusionlimited && !(kfs === nothing)
241193
len = length(phase.reactions)
242194
krev = zeros(typeof(N),len)
243195
@simd for i = 1:len
244-
@fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi;kf=kfs[i])
196+
@fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,dGrxns[i],diffs,V,phi,d;kf=kfs[i])
245197
end
246198
else
247199
len = length(phase.reactions)
248200
kfs = zeros(typeof(N),len)
249201
krev = zeros(typeof(N),len)
250202
@simd for i = 1:len
251-
@fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi)
203+
@fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,dGrxns[i],diffs,V,phi,d)
252204
end
253205
end
254206
kfs .*= phase.forwardability
255207
krev .*= phase.reversibility
256208
return kfs,krev
257209
end
258210

259-
@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
211+
@inline function getkfkrevs(phase::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Q2,diffs::Q3,V::W5,phi::W7,d;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
212+
dGrxns = -(phase.stoichmatrix*Gs).+phase.electronchange.*(phi*F)
260213
if !phase.diffusionlimited && kfs === nothing
261-
kfs = getkfs(phase,T,P,C,ns,V,phi)
214+
kfs = getkfs(phase,T,P,C,ns,V,phi,dGrxns,d)
262215
if phi == 0.0
263-
krev = @fastmath kfs./getKcs(phase,T,Gs)
216+
krev = @fastmath kfs./getKcs(phase,T,dGrxns)
264217
else
265-
krev = @fastmath kfs./getKcs(phase,T,Gs,phi)
218+
krev = @fastmath kfs./getKcs(phase,T,dGrxns,phi)
266219
end
267220
elseif !phase.diffusionlimited && !(kfs === nothing)
268221
if phi == 0.0
269-
krev = @fastmath kfs./getKcs(phase,T,Gs)
222+
krev = @fastmath kfs./getKcs(phase,T,dGrxns)
270223
else
271-
krev = @fastmath kfs./getKcs(phase,T,Gs,phi)
224+
krev = @fastmath kfs./getKcs(phase,T,dGrxns,phi)
272225
end
273226
elseif phase.diffusionlimited && !(kfs === nothing)
274227
len = length(phase.reactions)
275228
krev = similar(kfs)
276229
kfsdiff = similar(kfs)
277230
@simd for i = 1:len
278-
@fastmath @inbounds kfsdiff[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi;kf=kfs[i])
231+
@fastmath @inbounds kfsdiff[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,dGrxns[i],diffs,V,phi,d;kf=kfs[i])
279232
end
280233
return kfsdiff, krev
281234
else
282235
len = length(phase.reactions)
283-
kfs = zeros(typeof(Gs[1]),len)
236+
kfs = zeros(typeof(Gs[1]),len)ss
284237
krev = zeros(typeof(Gs[1]),len)
285238
@simd for i = 1:len
286-
@fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi)
239+
@fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,dGrxns[i],diffs,V,phi,d)
287240
end
288241
end
289242
return kfs,krev
290243
end
291244

292-
@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
245+
@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,d;kfs::W6=nothing) where {U,W7,W6,W5<:Real,W1<:Real,W2<:Real,W3,W4,Q1<:AbstractArray,Q2<:ForwardDiff.Dual,Q3<:AbstractArray} #autodiff p
246+
dGrxns = -(phase.stoichmatrix*Gs).+phase.electronchange.*(phi*F)
293247
if !phase.diffusionlimited && kfs === nothing
294-
kfs = getkfs(phase,T,P,C,ns,V,phi)
248+
kfs = getkfs(phase,T,P,C,ns,V,phi,dGrxns,d)
295249
if phi == 0.0
296-
krev = @fastmath kfs./getKcs(phase,T,Gs)
250+
krev = @fastmath kfs./getKcs(phase,T,dGrxns)
297251
else
298-
krev = @fastmath kfs./getKcs(phase,T,Gs,phi)
252+
krev = @fastmath kfs./getKcs(phase,T,dGrxns,phi)
299253
end
300254
elseif !phase.diffusionlimited && !(kfs === nothing)
301255
if phi == 0.0
302-
krev = @fastmath kfs./getKcs(phase,T,Gs)
256+
krev = @fastmath kfs./getKcs(phase,T,dGrxns)
303257
else
304-
krev = @fastmath kfs./getKcs(phase,T,Gs,phi)
258+
krev = @fastmath kfs./getKcs(phase,T,dGrxns,phi)
305259
end
306260
elseif phase.diffusionlimited && !(kfs === nothing)
307261
len = length(phase.reactions)
308262
krev = similar(kfs)
309263
@simd for i = 1:len
310-
@fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi;kf=kfs[i])
264+
@fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,dGrxns[i],diffs,V,phi,d;kf=kfs[i])
311265
end
312266
else
313267
len = length(phase.reactions)
314268
kfs = zeros(typeof(Gs[1]),len)
315269
krev = zeros(typeof(Gs[1]),len)
316270
@simd for i = 1:len
317-
@fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi)
271+
@fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,dGrxns[i],diffs,V,phi,d)
318272
end
319273
end
320274
kfs .*= phase.forwardability

0 commit comments

Comments
 (0)