Skip to content

Enabling modeling 5=5 reactions #201

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.

Sign up for GitHub

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

Open
wants to merge 4 commits into
base: main
Choose a base branch
from
Open

Enabling modeling 5=5 reactions #201

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
b85ccb5
Enabling five party reactions
xiaoruiDong May 20, 2022
0300330
enable 5 party reaction for edge analysis
xiaoruiDong May 20, 2022
4d29b0e
Clean up Reactor.jl format
xiaoruiDong Jun 9, 2022
fca2fff
BugFix: correct sim.p
xiaoruiDong Jun 10, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
198 changes: 129 additions & 69 deletions src/Domain.jl
View file
Open in desktop

Large diffs are not rendered by default.

26 changes: 13 additions & 13 deletions src/EdgeAnalysis.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -227,8 +227,8 @@ function getkeyselectioninds(coreedgedomains,coreedgeinters,domains,inters)
@inbounds rxnindexedge += length(coreedgedomains[i].phase.reactions)

@inbounds indend = length(domains[i].phase.reactions)
@inbounds reactantindices[:,ind:ind+indend-1] = domains[i].rxnarray[1:4,:]
@inbounds productindices[:,ind:ind+indend-1] = domains[i].rxnarray[5:8,:]
@inbounds reactantindices[:,ind:ind+indend-1] = domains[i].rxnarray[1:5,:]
@inbounds productindices[:,ind:ind+indend-1] = domains[i].rxnarray[6:10,:]
ind += indend
end

Expand All @@ -239,8 +239,8 @@ function getkeyselectioninds(coreedgedomains,coreedgeinters,domains,inters)
@inbounds index += length(coreedgeinters[i].phase.reactions)

@inbounds indend = length(inters[i].reactions)
@inbounds reactantindices[:,ind:ind+indend] = inters[i].rxnarray[1:4,:]
@inbounds productindices[:,ind:ind+indend] = inters[i].rxnarray[5:8,:]
@inbounds reactantindices[:,ind:ind+indend] = inters[i].rxnarray[1:5,:]
@inbounds productindices[:,ind:ind+indend] = inters[i].rxnarray[6:10,:]
ind += indend
end
end
Expand Down Expand Up @@ -270,8 +270,8 @@ function getkeyselectioninds(coreedgedomain::AbstractDomain,coreedgeinters,domai
end
end

@views @inbounds reactantindices = coreedgedomain.rxnarray[1:4,:]
@views @inbounds productindices = coreedgedomain.rxnarray[5:8,:]
@views @inbounds reactantindices = coreedgedomain.rxnarray[1:5,:]
@views @inbounds productindices = coreedgedomain.rxnarray[6:10,:]
coretoedgespcmap = Dict([i=>findfirst(isequal(spc),coreedgedomain.phase.species) for (i,spc) in enumerate(domain.phase.species)])
@simd for j = 3:length(domain.indexes)
@inbounds coretoedgespcmap[domain.indexes[j]] = coreedgedomain.indexes[j]
Expand Down Expand Up @@ -305,15 +305,15 @@ function processfluxes(sim::SystemSimulation,
if @inbounds any(d.rxnarray[:,i].>length(corespeciesconcentrations))
continue
end
for j = 1:4
for j = 1:5
if @inbounds d.rxnarray[j,i] != 0
@inbounds corespeciesconsumptionrates[d.rxnarray[j,i]] += frts[i+index]
@inbounds corespeciesproductionrates[d.rxnarray[j,i]] += rrts[i+index]
else
break
end
end
for j = 5:8
for j = 6:10
if @inbounds d.rxnarray[j,i] != 0
@inbounds corespeciesproductionrates[d.rxnarray[j,i]] += frts[i+index]
@inbounds corespeciesconsumptionrates[d.rxnarray[j,i]] += rrts[i+index]
Expand All @@ -330,15 +330,15 @@ function processfluxes(sim::SystemSimulation,
if @inbounds any(d.rxnarray[:,i].>length(corespeciesconcentrations))
continue
end
for j = 1:4
for j = 1:5
if @inbounds d.rxnarray[j,i] != 0
@inbounds corespeciesconsumptionrates[d.rxnarray[j,i]] += frts[i+index]
@inbounds corespeciesproductionrates[d.rxnarray[j,i]] += rrts[i+index]
else
break
end
end
for j = 5:8
for j = 6:10
if @inbounds d.rxnarray[j,i] != 0
@inbounds corespeciesproductionrates[d.rxnarray[j,i]] += frts[i+index]
@inbounds corespeciesconsumptionrates[d.rxnarray[j,i]] += rrts[i+index]
Expand Down Expand Up @@ -378,15 +378,15 @@ function processfluxes(sim::Simulation,corespcsinds,corerxninds,edgespcsinds,edg
if @inbounds any(d.rxnarray[:,i].>length(corespeciesconcentrations))
continue
end
for j = 1:4
for j = 1:5
if @inbounds d.rxnarray[j,i] != 0
@inbounds corespeciesconsumptionrates[d.rxnarray[j,i]] += frts[i]
@inbounds corespeciesproductionrates[d.rxnarray[j,i]] += rrts[i]
else
break
end
end
for j = 5:8
for j = 6:10
if @inbounds d.rxnarray[j,i] != 0
@inbounds corespeciesproductionrates[d.rxnarray[j,i]] += frts[i]
@inbounds corespeciesconsumptionrates[d.rxnarray[j,i]] += rrts[i]
Expand Down Expand Up @@ -449,7 +449,7 @@ function calcbranchingnumbers(sim,reactantinds,productinds,corespcsinds,corerxni
end
end
end
return branchingnums
return branchingnums
end

export calcbranchingnumbers
Expand Down
6 changes: 3 additions & 3 deletions src/Phase.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -147,7 +147,7 @@ function getstoichmatrix(rxnarray,spcs)
for (j,ind) in enumerate(rxnarray[:,i])
if ind == 0
continue
elseif j > 3
elseif j > 5
M[i,ind] -= 1.0
n += 1.0
else
Expand Down Expand Up @@ -257,7 +257,7 @@ Broadcast.broadcastable(p::T) where {T<:AbstractPhase} = Ref(p)
export broadcastable

function getreactionindices(spcs,rxns) where {Q<:AbstractPhase}
arr = zeros(Int64,(8,length(rxns)))
arr = zeros(Int64,(10,length(rxns)))
names = [spc.name for spc in spcs]
for (i,rxn) in enumerate(rxns)
inds = [findfirst(isequal(spc),spcs) for spc in rxn.reactants]
Expand All @@ -268,7 +268,7 @@ function getreactionindices(spcs,rxns) where {Q<:AbstractPhase}
end
for (j,spc) in enumerate(rxn.products)
ind = findfirst(isequal(spc),spcs)
arr[j+4,i] = ind
arr[j+5,i] = ind
rxn.productinds[j] = ind
end
if hasproperty(rxn.kinetics,:efficiencies) && length(rxn.kinetics.nameefficiencies) > 0
Expand Down
16 changes: 10 additions & 6 deletions src/PhaseState.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -110,20 +110,24 @@ export getDiffusiveRate
if Nreact == 1
@fastmath @inbounds dGrxn -= Gs[rxn.reactantinds[1]]
elseif Nreact == 2
@fastmath @inbounds dGrxn -= Gs[rxn.reactantinds[1]]+Gs[rxn.reactantinds[2]]
@fastmath @inbounds dGrxn -= Gs[rxn.reactantinds[1]] + Gs[rxn.reactantinds[2]]
elseif Nreact == 3
@fastmath @inbounds dGrxn -= Gs[rxn.reactantinds[1]]+Gs[rxn.reactantinds[2]]+Gs[rxn.reactantinds[3]]
@fastmath @inbounds dGrxn -= Gs[rxn.reactantinds[1]] + Gs[rxn.reactantinds[2]] + Gs[rxn.reactantinds[3]]
elseif Nreact == 4
@fastmath @inbounds dGrxn -= Gs[rxn.reactantinds[1]]+Gs[rxn.reactantinds[2]]+Gs[rxn.reactantinds[3]]+Gs[rxn.reactantinds[4]]
@fastmath @inbounds dGrxn -= Gs[rxn.reactantinds[1]] + Gs[rxn.reactantinds[2]] + Gs[rxn.reactantinds[3]] + Gs[rxn.reactantinds[4]]
elseif Nreact == 5
@fastmath @inbounds dGrxn -= Gs[rxn.reactantinds[1]] + Gs[rxn.reactantinds[2]] + Gs[rxn.reactantinds[3]] + Gs[rxn.reactantinds[4]] + Gs[rxn.reactantinds[5]]
end
if Nprod == 1
@fastmath @inbounds dGrxn += Gs[rxn.productinds[1]]
elseif Nprod == 2
@fastmath @inbounds dGrxn += Gs[rxn.productinds[1]]+Gs[rxn.productinds[2]]
@fastmath @inbounds dGrxn += Gs[rxn.productinds[1]] + Gs[rxn.productinds[2]]
elseif Nprod == 3
@fastmath @inbounds dGrxn += Gs[rxn.productinds[1]]+Gs[rxn.productinds[2]]+Gs[rxn.productinds[3]]
@fastmath @inbounds dGrxn += Gs[rxn.productinds[1]] + Gs[rxn.productinds[2]] + Gs[rxn.productinds[3]]
elseif Nprod == 4
@fastmath @inbounds dGrxn += Gs[rxn.productinds[1]]+Gs[rxn.productinds[2]]+Gs[rxn.productinds[3]]+Gs[rxn.productinds[4]]
@fastmath @inbounds dGrxn += Gs[rxn.productinds[1]] + Gs[rxn.productinds[2]] + Gs[rxn.productinds[3]] + Gs[rxn.productinds[4]]
elseif Nprod == 5
@fastmath @inbounds dGrxn += Gs[rxn.productinds[1]] + Gs[rxn.productinds[2]] + Gs[rxn.productinds[3]] + Gs[rxn.productinds[4]] + Gs[rxn.productinds[5]]
end
return @inbounds @fastmath exp(-(dGrxn+rxn.electronchange*phi)/(R*T))*(getC0(ph,T))^(Nprod-Nreact)
end
Expand Down
Loading