Add basic (NE)KMC code (#1)

src/NICE.jl

@@ -1,5 +1,180 @@
 module NICE
 
-# Write your package code here.
+
+export ReactionSystem, run_kmc, run_nekmc
+
+
+mutable struct ReactionSystem
+
+    # Constants
+    n_species::Int
+    n_reaction::Int
+    stoich::Matrix{Float64}
+    rev_rate_consts::Vector{Float64}
+    fwd_rate_consts::Vector{Float64}
+
+    # Mutables
+    n_iter::Int
+    time::Float64
+    concs::Vector{Float64}
+    rev_rates::Vector{Float64}
+    fwd_rates::Vector{Float64}
+    net_rates::Vector{Float64}
+end
+
+
+function ReactionSystem(stoich, concs, rev_rate_consts, fwd_rate_consts)
+    n_species = size(stoich, 1)
+    n_reaction = size(stoich, 2)
+    n_iter = 0
+    time = 0.
+    rev_rates = zeros(Float64, n_reaction)
+    fwd_rates = zeros(Float64, n_reaction)
+    net_rates = zeros(Float64, n_reaction)
+    return ReactionSystem(n_species, n_reaction, stoich, rev_rate_consts, fwd_rate_consts,
+        n_iter, time, concs, rev_rates, fwd_rates, net_rates)
+end
+
+
+function ReactionSystem(stoich, concs, keq_vals; φ::Real=1.0)
+    n_species = size(stoich, 1)
+    n_reaction = size(stoich, 2)
+    rev_rate_consts = φ ./ (keq_vals .+ 1)
+    fwd_rate_consts = φ .- rev_rate_consts
+    n_iter = 0
+    time = 0.
+    rev_rates = zeros(Float64, n_reaction)
+    fwd_rates = zeros(Float64, n_reaction)
+    net_rates = zeros(Float64, n_reaction)
+    return ReactionSystem(n_species, n_reaction, stoich, rev_rate_consts, fwd_rate_consts,
+        n_iter, time, concs, rev_rates, fwd_rates, net_rates)
+end
+
+
+function run_kmc(rxn_system, n_iter; ε=1.0e-4)
+    pvec = zeros(Float64, 2 * rxn_system.n_reaction)
+    for i in 1 : n_iter
+        update_rates(rxn_system)
+        i_rxn = kmc_select_reaction(rxn_system, pvec)
+        rxn_system.time += kmc_do_reaction(rxn_system, i_rxn, ε)
+        rxn_system.n_iter += 1
+    end
+end
+
+
+function run_nekmc(rxn_system, n_iter; ε=1.0e-4)
+    pvec = zeros(Float64, rxn_system.n_reaction)
+    for i in 1 : n_iter
+        update_rates(rxn_system)
+        i_rxn = nekmc_select_reaction(rxn_system, pvec)
+        rxn_system.time += nekmc_do_reaction(rxn_system, i_rxn, ε)
+        rxn_system.n_iter += 1
+    end
+end
+
+
+function update_rates(rxn_system)
+    for i in 1 : rxn_system.n_reaction
+        rev_rate = rxn_system.rev_rate_consts[i]
+        fwd_rate = rxn_system.fwd_rate_consts[i]
+        for j in 1 : rxn_system.n_species
+            s = rxn_system.stoich[j, i]
+            if s >= 0.
+                rev_rate *= rxn_system.concs[j] ^ s
+            else
+                fwd_rate *= rxn_system.concs[j] ^ abs(s)
+            end
+        end
+        rxn_system.rev_rates[i] = rev_rate
+        rxn_system.fwd_rates[i] = fwd_rate
+        rxn_system.net_rates[i] = fwd_rate - rev_rate
+    end
+end
+
+
+function kmc_select_reaction(rxn_system, pvec)
+    # Update probability vector
+    p = 0.
+    for i in 1 : rxn_system.n_reaction
+        p += rxn_system.fwd_rates[i]
+        pvec[i] = p
+    end
+    for i in 1 : rxn_system.n_reaction
+        p += rxn_system.rev_rates[i]
+        pvec[i + rxn_system.n_reaction] = p
+    end
+    # Select random reaction
+    p *= rand(Float64)
+    for i in 1 : 2 * rxn_system.n_reaction
+        if pvec[i] > p
+            return i
+        end
+    end
+    # sentinel return value
+    return 2 * rxn_system.n_reaction
+end
+
+
+function nekmc_select_reaction(rxn_system, pvec)
+    # Update probability vector
+    p = 0.
+    for i in 1 : rxn_system.n_reaction
+        p += abs(rxn_system.net_rates[i])
+        pvec[i] = p
+    end
+    # Select random reaction
+    p *= rand(Float64)
+    for i in 1 : rxn_system.n_reaction
+        if pvec[i] > p
+            return i
+        end
+    end
+    # sentinel return value
+    return rxn_system.n_reaction
+end
+
+
+function kmc_do_reaction(rxn_system, i_rxn, ε)
+    # Forward reaction (i_rxn <= n_rxn)
+    if i_rxn <= rxn_system.n_reaction
+        rate = rxn_system.fwd_rates[i_rxn]
+        if rate >= 0
+            for j = 1 : rxn_system.n_species
+                rxn_system.concs[j] += rxn_system.stoich[j, i_rxn] * ε
+            end
+            return ε / rate
+        end
+    # Reverse reaction (i_rxn > n_rxn)
+    else
+        i_rxn -= rxn_system.n_reaction
+        rate = rxn_system.rev_rates[i_rxn]
+        if rate >= 0
+            for j = 1 : rxn_system.n_species
+                rxn_system.concs[j] -= rxn_system.stoich[j, i_rxn] * ε
+            end
+            return ε / rate
+        end
+    end
+    # Sentinel return value
+    return 0.
+end
+
+
+function nekmc_do_reaction(rxn_system, i_rxn, ε)
+    # Net reaction
+    rate = rxn_system.net_rates[i_rxn]
+    if rate >= 0
+        for j = 1 : rxn_system.n_species
+            rxn_system.concs[j] += rxn_system.stoich[j, i_rxn] * ε
+        end
+        return ε / rate
+    else
+        for j = 1 : rxn_system.n_species
+            rxn_system.concs[j] -= rxn_system.stoich[j, i_rxn] * ε
+        end
+        return -ε / rate
+    end
+end
+
 
 end

test/runtests.jl

@@ -2,5 +2,43 @@ using NICE
 using Test
 
 @testset "NICE.jl" begin
-    # Write your tests here.
+
+    # Example 1
+
+    test_concs = [0.9261879203, 0.9623865752, 0.0926187920]
+
+    stoich = [[-0.5, 1.0, 0.0] [-0.5, -1.0, 1.0]]
+    keq_vals = [1.0, 0.1]
+
+    concs = [1.0, 0.2, 0.4]
+    rxn_system = ReactionSystem(stoich, concs, keq_vals, φ=1.0)
+
+    run_kmc(rxn_system, 1000000, ε=1.0e-5)
+    @test isapprox(rxn_system.concs, test_concs, atol=1.0e-5)
+
+    concs = [1.0, 0.2, 0.4]
+    rxn_system = ReactionSystem(stoich, concs, keq_vals, φ=1.0)
+
+    run_nekmc(rxn_system, 1000000, ε=1.0e-5)
+    @test isapprox(rxn_system.concs, test_concs, atol=1.0e-5)
+
+    # Example 2
+
+    test_concs = [9.0e-01, 1.10212630e-08, 1.00002433e-10, 9.98999891e-02, 0.0, 9.99999e-05]
+
+    stoich = [[-1, -1,  0,  1,  0,  0] [ 0, -1, -1,  0,  1,  0] [ 0,  0, -1, -1,  0,  1] [-1,  0,  0,  0, -1,  1]]
+    keq_vals = [1.0e7, 1.0e9, 1.0e7, 1.0e9]
+
+    concs = [1.0, 0.1, 1e-4, 0.0, 0.0, 0.0]
+    rxn_system = ReactionSystem(stoich, concs, keq_vals, φ=1.0)
+
+    run_kmc(rxn_system, 5000000, ε=1.0e-5)
+    @test isapprox(rxn_system.concs, test_concs, atol=1.0e-5)
+
+    concs = [1.0, 0.1, 1e-4, 0.0, 0.0, 0.0]
+    rxn_system = ReactionSystem(stoich, concs, keq_vals, φ=1.0)
+
+    run_nekmc(rxn_system, 1000000, ε=1.0e-5)
+    @test isapprox(rxn_system.concs, test_concs, atol=1.0e-5)
+
 end